Model is X_1,...,X_n ~ N(theta, s^2), theta~self, s fixed
326373
def posterior(self, x, s=1.):
        """"""
        pr0 = 1. / self.sigma**2  # prior precision
        prd = x.size / s**2  # data precision
        varp = 1. / (pr0 + prd)  # posterior variance
        mu = varp * (pr0 * self.mu + prd * x.mean())
        return TruncNormal(mu=mu, sigma=np.sqrt(varp), a=self.a, b=self.b)
414054
def set_model(self, f0, A, model='trigonometric'):
        
        self.f0 = f0
        self.y_clean, vary = trigonometric_model(self.t, f0, A)
        # y_clean is normalized to have unit standard deviation
        self.y_clean = self.y_clean/np.sqrt(vary)
635337
def pol2cart(theta, rho):
    """"""
    x = rho * np.cos(theta)
    y = rho * np.sin(theta)
    return x, y
463897
def gradF_t(X, Y, S, M_E, E, m0, rho):
    ''' 
    '''
    n, r = X.shape
    m, r = Y.shape

    XS = np.dot(X, S)
    YS = np.dot(Y, S.T)
    XSY = np.dot(XS, Y.T)

    Qx = np.dot(np.dot(X.T, ((M_E - XSY) * E)), YS) / n
    Qy = np.dot(np.dot(Y.T, ((M_E - XSY) * E).T), XS) / m

    W = np.dot((XSY - M_E) * E, YS) + np.dot(X, Qx) + rho * Gp(X, m0, r)
    Z = np.dot(((XSY - M_E) * E).T, XS) + np.dot(Y, Qy) + rho * Gp(Y, m0, r)

    return W, Z
755854
def _fixedpoint(D, tol=1e-7, maxiter=None):
    ''''''
    N, K = D.shape
    logp = log(D).mean(axis=0)
    a0 = _init_a(D)

    # Start updating
    if maxiter is None:
        maxiter = MAXINT
    for i in xrange(maxiter):
        a1 = _ipsi(psi(a0.sum()) + logp)
        # if norm(a1-a0) < tol:
        if abs(loglikelihood(D, a1)-loglikelihood(D, a0)) < tol: # much faster
            return a1
        a0 = a1
    raise Exception('Failed to converge after {} iterations, values are {}.'
                    .format(maxiter, a1))
884947
def multi_normal(X, t):
    
    m, L = update_params(X, t)
    z = np.random.standard_normal(np.shape(m)) # generate i.i.d N(0,1)
    return  la.solve(L.T,z)+m
1247906
def xstep(self):
        
        """

        ngsit = 0
        gsrrs = np.inf
        while gsrrs > self.opt['GSTol'] and ngsit < self.opt['MaxGSIter']:
            self.X = self.GaussSeidelStep(self.S, self.X,
                                          self.cnst_AT(self.Y-self.U),
                                          self.rho, self.lcw, self.Wdf2)
            gsrrs = sl.rrs(
                self.rho*self.cnst_AT(self.cnst_A(self.X)) +
                self.Wdf2*self.X, self.Wdf2*self.S +
                self.rho*self.cnst_AT(self.Y - self.U))
            ngsit += 1

        self.xs = (ngsit, gsrrs)
137140
def cbpdnmd_ustep(k):
    

    mp_Z_U0[k] += mp_DX[k] - mp_Z_Y0[k] - mp_S[k]
    mp_Z_U1[k] += mp_Z_X[k] - mp_Z_Y1[k]
855572
function Model(inputMatrix, rowLabels, colLabels) {
	this.rowLabels = rowLabels;	// labels for the rows
	this.colLabels = colLabels; // labels for the columns
	this.input = inputMatrix;	// input data
	
	// estimated data, initialized to all zeros
	this.estimated = sylvester.Matrix.Zeros(this.input.rows(),this.input.cols());
}
492780
def as_2d_array(theta):
    
    v = theta.view(np.float)
    N = theta.shape[0]
    v.shape = (N, - 1)
    # raise an error if v cannot be reshaped without creating a copy
    return v
294129
function mpn_shr1(r, n, l) {
	var c = 0;
	while (l > 0) {
		var w = n[--l] | c;
		c = w << 31;
		r[l] = w >>> 1;
	}
	return c;
}
350379
def K(self, X, X2=None):        
    # model :   -a d^2y/dx^2  + b dy/dt + c * y = U
    # kernel Kyy rbf spatiol temporal
    # vyt Y temporal variance  vyx Y spatiol variance   lyt Y temporal lengthscale   lyx Y spatiol lengthscale
    # kernel Kuu doper( doper(Kyy))
    # a   b    c    lyt   lyx    vyx*vyt
        """"""        
        X,slices = X[:,:-1],index_to_slices(X[:,-1])
        if X2 is None:
            X2,slices2 = X,slices
            K = np.zeros((X.shape[0], X.shape[0]))
        else:
            X2,slices2 = X2[:,:-1],index_to_slices(X2[:,-1])
            K = np.zeros((X.shape[0], X2.shape[0]))


        tdist = (X[:,0][:,None] - X2[:,0][None,:])**2
        xdist = (X[:,1][:,None] - X2[:,1][None,:])**2

        ttdist = (X[:,0][:,None] - X2[:,0][None,:])
        #rdist = [tdist,xdist]
        #dist = np.abs(X - X2.T)
        vyt = self.variance_Yt
        vyx = self.variance_Yx
        
        lyt=1/(2*self.lengthscale_Yt)
        lyx=1/(2*self.lengthscale_Yx)

        a = self.a ## -a is used in the model, negtive diffusion
        b = self.b
        c = self.c

        kyy = lambda tdist,xdist: np.exp(-lyt*(tdist) -lyx*(xdist))

        k1 = lambda tdist: (2*lyt - 4*lyt**2 * (tdist) )

        k2 = lambda xdist: ( 4*lyx**2 * (xdist)  - 2*lyx )

        k3 = lambda xdist: ( 3*4*lyx**2 - 6*8*xdist*lyx**3 + 16*xdist**2*lyx**4 )

        k4 = lambda ttdist: 2*lyt*(ttdist)

        for i, s1 in enumerate(slices):
            for j, s2 in enumerate(slices2):
                for ss1 in s1:
                    for ss2 in s2:
                        if i==0 and j==0:
                            K[ss1,ss2] = vyt*vyx*kyy(tdist[ss1,ss2],xdist[ss1,ss2])
                        elif i==0 and j==1:
                            K[ss1,ss2] = (-a*k2(xdist[ss1,ss2]) + b*k4(ttdist[ss1,ss2]) + c)*vyt*vyx*kyy(tdist[ss1,ss2],xdist[ss1,ss2])
                            #K[ss1,ss2]=  np.where(  rdist[ss1,ss2]>0 , kuyp(np.abs(rdist[ss1,ss2])), kuyn(np.abs(rdist[ss1,ss2]) )   )
                            #K[ss1,ss2]=  np.where(  rdist[ss1,ss2]>0 , kuyp(rdist[ss1,ss2]), kuyn(rdist[ss1,ss2] )   )
                        elif i==1 and j==1:
                            K[ss1,ss2] = ( b**2*k1(tdist[ss1,ss2]) - 2*a*c*k2(xdist[ss1,ss2]) + a**2*k3(xdist[ss1,ss2]) + c**2 )* vyt*vyx* kyy(tdist[ss1,ss2],xdist[ss1,ss2])
                        else:
                            K[ss1,ss2] = (-a*k2(xdist[ss1,ss2]) - b*k4(ttdist[ss1,ss2]) + c)*vyt*vyx*kyy(tdist[ss1,ss2],xdist[ss1,ss2])
                            #K[ss1,ss2]= np.where(  rdist[ss1,ss2]>0 , kyup(np.abs(rdist[ss1,ss2])), kyun(np.abs(rdist[ss1,ss2]) )   )
                            #K[ss1,ss2] = np.where(  rdist[ss1,ss2]>0 , kyup(rdist[ss1,ss2]), kyun(rdist[ss1,ss2] )   )
        
        #stop
        return K
315240
def cv_compute(self, b, A, B, C, mK, f, m1, m2):
        '''
        

        '''

        A = np.sum([l * a for l, a in zip(self.lam[b], A)
                    if l is not None], axis=0)
        B = np.sum([l * b for l, b in zip(self.lam[b], B)
                    if l is not None], axis=0)
        W = np.linalg.solve(mK + A + C, f)
        if self.transit_model is None:
            model = np.dot(B, W)
        else:
            w_pld = np.concatenate([l * np.dot(self.X(n, m2).T, W)
                                    for n, l in enumerate(self.lam[b])
                                    if l is not None])
            model = np.dot(np.hstack(
                [self.X(n, m1) for n, l in enumerate(self.lam[b])
                 if l is not None]), w_pld)
        model -= np.nanmedian(model)

        return model
604224
def B_(self):
        '''
        (n_, n) = self.H.shape
        ls = self.linear_system
        if self._B_ is None or self._B_.shape[1] < n_:
            # compute B_
            if ls.self_adjoint:
                self._B_ = self.C.T.conj()
                if n_ > n:
                    self._B_ = numpy.r_[self._B_,
                                        utils.inner(self.V[:, [-1]],
                                                    self.projection.AU,
                                                    ip_B=ls.ip_B)]
            else:
                self._B_ = utils.inner(self.V, self.projection.AU,
                                       ip_B=ls.ip_B)
        return self._B_
844847
def hilbertrot(n, x, y, rx, ry):
    
    if ry == 0:
        if rx == 1:
            x = n - 1 - x
            y = n - 1 - y
        return y, x
    return x, y
167683
def fit(self):
        
        # solve for estimator of B
        n,p          = np.shape(self._X)
        self._df     = float(n - p)
        self._Cx     = np.linalg.pinv(np.dot(self._X.T,self._X))
        self._Bhat   = np.dot(np.dot(self._Cx, self._X.T), self._A)
        self._Y_rec  = self._compute_prediction(self._X)
302153
def funTransSectY(theta, rho):
    
    return np.matrix([[1, rho * theta], [0, 1]], dtype=np.double)