Example of regressing with a student t likelihood using Laplace
500057
def student_t_approx(optimize=True, plot=True):
    
    real_std = 0.1
    #Start a function, any function
    X = np.linspace(0.0, np.pi*2, 100)[:, None]
    Y = np.sin(X) + np.random.randn(*X.shape)*real_std
    Y = Y/Y.max()
    Yc = Y.copy()

    X_full = np.linspace(0.0, np.pi*2, 500)[:, None]
    Y_full = np.sin(X_full)
    Y_full = Y_full/Y_full.max()

    #Slightly noisy data
    Yc[75:80] += 1

    #Very noisy data
    #Yc[10] += 100
    #Yc[25] += 10
    #Yc[23] += 10
    #Yc[26] += 1000
    #Yc[24] += 10
    #Yc = Yc/Yc.max()

    #Add student t random noise to datapoints
    deg_free = 1
    print("Real noise: ", real_std)
    initial_var_guess = 0.5
    edited_real_sd = initial_var_guess

    # Kernel object
    kernel1 = GPy.kern.RBF(X.shape[1]) + GPy.kern.White(X.shape[1])
    kernel2 = GPy.kern.RBF(X.shape[1]) + GPy.kern.White(X.shape[1])
    kernel3 = GPy.kern.RBF(X.shape[1]) + GPy.kern.White(X.shape[1])
    kernel4 = GPy.kern.RBF(X.shape[1]) + GPy.kern.White(X.shape[1])

    #Gaussian GP model on clean data
    m1 = GPy.models.GPRegression(X, Y.copy(), kernel=kernel1)
    # optimize
    m1['.*white'].constrain_fixed(1e-5)
    m1.randomize()

    #Gaussian GP model on corrupt data
    m2 = GPy.models.GPRegression(X, Yc.copy(), kernel=kernel2)
    m2['.*white'].constrain_fixed(1e-5)
    m2.randomize()

    #Student t GP model on clean data
    t_distribution = GPy.likelihoods.StudentT(deg_free=deg_free, sigma2=edited_real_sd)
    laplace_inf = GPy.inference.latent_function_inference.Laplace()
    m3 = GPy.core.GP(X, Y.copy(), kernel3, likelihood=t_distribution, inference_method=laplace_inf)
    m3['.*t_scale2'].constrain_bounded(1e-6, 10.)
    m3['.*white'].constrain_fixed(1e-5)
    m3.randomize()

    #Student t GP model on corrupt data
    t_distribution = GPy.likelihoods.StudentT(deg_free=deg_free, sigma2=edited_real_sd)
    laplace_inf = GPy.inference.latent_function_inference.Laplace()
    m4 = GPy.core.GP(X, Yc.copy(), kernel4, likelihood=t_distribution, inference_method=laplace_inf)
    m4['.*t_scale2'].constrain_bounded(1e-6, 10.)
    m4['.*white'].constrain_fixed(1e-5)
    m4.randomize()
    print(m4)
    debug=True
    if debug:
        m4.optimize(messages=1)
        from matplotlib import pyplot as pb
        pb.plot(m4.X, m4.inference_method.f_hat)
        pb.plot(m4.X, m4.Y, 'rx')
        m4.plot()
        print(m4)
        return m4

    if optimize:
        optimizer='scg'
        print("Clean Gaussian")
        m1.optimize(optimizer, messages=1)
        print("Corrupt Gaussian")
        m2.optimize(optimizer, messages=1)
        print("Clean student t")
        m3.optimize(optimizer, messages=1)
        print("Corrupt student t")
        m4.optimize(optimizer, messages=1)

    if plot:
        plt.figure(1)
        plt.suptitle('Gaussian likelihood')
        ax = plt.subplot(211)
        m1.plot(ax=ax)
        plt.plot(X_full, Y_full)
        plt.ylim(-1.5, 1.5)
        plt.title('Gaussian clean')

        ax = plt.subplot(212)
        m2.plot(ax=ax)
        plt.plot(X_full, Y_full)
        plt.ylim(-1.5, 1.5)
        plt.title('Gaussian corrupt')

        plt.figure(2)
        plt.suptitle('Student-t likelihood')
        ax = plt.subplot(211)
        m3.plot(ax=ax)
        plt.plot(X_full, Y_full)
        plt.ylim(-1.5, 1.5)
        plt.title('Student-t rasm clean')

        ax = plt.subplot(212)
        m4.plot(ax=ax)
        plt.plot(X_full, Y_full)
        plt.ylim(-1.5, 1.5)
        plt.title('Student-t rasm corrupt')

    return m1, m2, m3, m4
913300
def log_likelihood(C, T):
    
    C = C.tocsr()
    T = T.tocsr()
    ind = scipy.nonzero(C)
    relT = np.array(T[ind])[0, :]
    relT = np.log(relT)
    relC = np.array(C[ind])[0, :]

    return relT.dot(relC)
518588
def _backward_step(self, deltaY, deltaY_sigma2inv_rho_weightT,
                       sigma2_e, weight, mu, mu_Gamma_inv, Gamma_inv,
                       Lambda_0, Lambda_1, H):
        
        n_T = len(Gamma_inv)
        # All the terms with hat before are parameters of posterior
        # distributions of X conditioned on data from all time points,
        # whereas the ones without hat calculated by _forward_step
        # are mean and covariance of posterior of X conditioned on
        # data up to the time point.
        Gamma_inv_hat = [None] * n_T
        mu_Gamma_inv_hat = [None] * n_T
        mu_hat = [None] * n_T

        mu_hat[-1] = mu[-1].copy()
        mu_Gamma_inv_hat[-1] = mu_Gamma_inv[-1].copy()
        Gamma_inv_hat[-1] = Gamma_inv[-1].copy()

        for t in np.arange(n_T - 2, -1, -1):
            tmp = np.linalg.solve(Gamma_inv_hat[t + 1] - Gamma_inv[t + 1]
                                  + Lambda_1, H)
            Gamma_inv_hat[t] = Gamma_inv[t] + Lambda_0 - np.dot(H.T, tmp)
            mu_Gamma_inv_hat[t] = mu_Gamma_inv[t] \
                - deltaY_sigma2inv_rho_weightT[t, :] + np.dot(
                    mu_Gamma_inv_hat[t + 1] - mu_Gamma_inv[t + 1]
                    + np.dot(deltaY[t, :] / sigma2_e, weight.T), tmp)
            mu_hat[t] = np.linalg.solve(Gamma_inv_hat[t],
                                        mu_Gamma_inv_hat[t])
        return mu_hat, mu_Gamma_inv_hat, Gamma_inv_hat
420853
def _laplace_fit(self,obj_type):
        

        # Get Mode and Inverse Hessian information
        y = self.fit(method='PML',printer=False)

        if y.ihessian is None:
            raise Exception("No Hessian information - Laplace approximation cannot be performed")
        else:

            self.latent_variables.estimation_method = 'Laplace'

            theta, Y, scores, states, states_var, X_names = self._categorize_model_output(self.latent_variables.get_z_values())

            # Change this in future
            try:
                latent_variables_store = self.latent_variables.copy()
            except:
                latent_variables_store = self.latent_variables

            return LaplaceResults(data_name=self.data_name,X_names=X_names,model_name=self.model_name,
                model_type=self.model_type, latent_variables=latent_variables_store,data=Y,index=self.index,
                multivariate_model=self.multivariate_model,objective_object=obj_type, 
                method='Laplace',ihessian=y.ihessian,signal=theta,scores=scores,
                z_hide=self._z_hide,max_lag=self.max_lag,states=states,states_var=states_var)
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)
892160
def fit_mle(self, data, b=None):
        

        data = np.array(data)
        length = len(data)

        if not b:
            b = np.sum(data)

        return _trunc_logser_solver(length, b), b
383895
def samples(self, gp, Y_metadata=None):
        
        orig_shape = gp.shape
        gp = gp.flatten()
        #FIXME: Very slow as we are computing a new random variable per input!
        #Can't get it to sample all at the same time
        #student_t_samples = np.array([stats.t.rvs(self.v, self.gp_link.transf(gpj),scale=np.sqrt(self.sigma2), size=1) for gpj in gp])
        dfs = np.ones_like(gp)*self.v
        scales = np.ones_like(gp)*np.sqrt(self.sigma2)
        student_t_samples = stats.t.rvs(dfs, loc=self.gp_link.transf(gp),
                                        scale=scales)
        return student_t_samples.reshape(orig_shape)
786755
def fit_deriv(x, amplitude, mean, stddev):
        
        import operator
        return list(map(
            operator.neg, Gaussian1D.fit_deriv(x, amplitude, mean, stddev)))
762343
def estimate(self, maxiter=250, convergence=1e-7):

        

        self.loglik = np.zeros(maxiter)

        iter = 0

        while iter < maxiter:

            self.loglik[iter] = self.E_step()
            if np.isnan(self.loglik[iter]):
                    print("undefined log-likelihood")
                    break
            self.M_step()

            if self.loglik[iter] - self.loglik[iter - 1] < 0 and iter > 0:
                    print("log-likelihood decreased by %f at iteration %d"
                          % (self.loglik[iter] - self.loglik[iter - 1],
                             iter))
            elif self.loglik[iter] - self.loglik[iter - 1] < convergence \
                    and iter > 0:
                    print("convergence at iteration %d, loglik = %f" %
                           (iter, self.loglik[iter]))
                    self.loglik = self.loglik[self.loglik < 0]
                    break

            iter += 1
1236690
def _EM_gradient(A,Q,H,R,m_init,P_init,measurements, M, P, G, dA, dQ, dH, dR, dm_init, dP_init):
        
        import pdb; pdb.set_trace();

        param_number = dA.shape[-1]
        d_log_likelihood_update = np.empty((param_number,1))

        sample_no = measurements.shape[0]
        P_1 = P[1:,:,:] # remove 0-th step
        P_2 = P[0:-1,:,:] # remove 0-th step

        M_1 = M[1:,:] # remove 0-th step
        M_2 = M[0:-1,:] # remove the last step

        Sigma = np.mean(P_1,axis=0) + np.dot(M_1.T, M_1) / sample_no #
        Phi =   np.mean(P_2,axis=0) + np.dot(M_2.T, M_2) / sample_no #

        B = np.dot( measurements.T, M_1 )/ sample_no
        C =   (sp.einsum( 'ijk,ikl', P_1, G) + np.dot(M_1.T, M_2)) / sample_no #

#        C1 = np.zeros( (P_1.shape[1],P_1.shape[1]) )
#        for k in range(P_1.shape[0]):
#            C1 += np.dot(P_1[k,:,:],G[k,:,:]) + sp.outer( M_1[k,:], M_2[k,:] )
#        C1 = C1 / sample_no

        D = np.dot( measurements.T, measurements ) / sample_no

        try:
            P_init_inv = sp.linalg.inv(P_init)

            if np.max( np.abs(P_init_inv)) > 10e13:
                compute_P_init_terms = False
            else:
                compute_P_init_terms = True
        except np.linalg.LinAlgError:
            compute_P_init_terms = False

        try:
            Q_inv = sp.linalg.inv(Q)

            if np.max( np.abs(Q_inv)) > 10e13:
                compute_Q_terms = False
            else:
                compute_Q_terms = True
        except np.linalg.LinAlgError:
            compute_Q_terms = False

        try:
            R_inv = sp.linalg.inv(R)

            if np.max( np.abs(R_inv)) > 10e13:
                compute_R_terms = False
            else:
                compute_R_terms = True
        except np.linalg.LinAlgError:
            compute_R_terms = False


        d_log_likelihood_update = np.zeros((param_number,1))
        for j in range(param_number):
            if compute_P_init_terms:
                d_log_likelihood_update[j,:] -= 0.5 * np.sum(P_init_inv* dP_init[:,:,j].T ) #p #m

                M0_smoothed = M[0]; M0_smoothed.shape = (M0_smoothed.shape[0],1)
                tmp1 = np.dot( dP_init[:,:,j], np.dot( P_init_inv, (P[0,:,:] + sp.outer( (M0_smoothed - m_init), (M0_smoothed - m_init) )) )  ) #p #m
                d_log_likelihood_update[j,:] += 0.5 * np.sum(P_init_inv* tmp1.T )

                tmp2 = sp.outer( dm_init[:,j], M0_smoothed )
                tmp2 += tmp2.T
                d_log_likelihood_update[j,:] += 0.5 * np.sum(P_init_inv* tmp2.T )

            if compute_Q_terms:

                d_log_likelihood_update[j,:] -=  sample_no/2.0 * np.sum(Q_inv* dQ[:,:,j].T ) #m

                tmp1 = np.dot(C,A.T); tmp1 += tmp1.T; tmp1 = Sigma - tmp1 + np.dot(A, np.dot(Phi,A.T)) #m
                tmp1 = np.dot( dQ[:,:,j], np.dot( Q_inv, tmp1) )
                d_log_likelihood_update[j,:] += sample_no/2.0 * np.sum(Q_inv * tmp1.T)

                tmp2 = np.dot( dA[:,:,j], C.T); tmp2 += tmp2.T;
                tmp3 = np.dot(dA[:,:,j], np.dot(Phi,A.T)); tmp3 += tmp3.T
                d_log_likelihood_update[j,:] -= sample_no/2.0 * np.sum(Q_inv.T * (tmp3 - tmp2) )

            if compute_R_terms:
                d_log_likelihood_update[j,:] -=  sample_no/2.0 * np.sum(R_inv* dR[:,:,j].T )

                tmp1 = np.dot(B,H.T); tmp1 += tmp1.T; tmp1 = D - tmp1 + np.dot(H, np.dot(Sigma,H.T))
                tmp1 = np.dot( dR[:,:,j], np.dot( R_inv, tmp1) )
                d_log_likelihood_update[j,:] += sample_no/2.0 * np.sum(R_inv * tmp1.T)

                tmp2 = np.dot( dH[:,:,j], B.T); tmp2 += tmp2.T;
                tmp3 = np.dot(dH[:,:,j], np.dot(Sigma,H.T)); tmp3 += tmp3.T
                d_log_likelihood_update[j,:] -= sample_no/2.0 * np.sum(R_inv.T * (tmp3 - tmp2) )

        return d_log_likelihood_update
1309185
def regress_poly(degree, data, remove_mean=True, axis=-1):
    
    timepoints = data.shape[0]
    # Generate design matrix
    X = np.ones((timepoints, 1))  # quick way to calc degree 0
    for i in range(degree):
        polynomial_func = Legendre.basis(i + 1)
        value_array = np.linspace(-1, 1, timepoints)
        X = np.hstack((X, polynomial_func(value_array)[:, np.newaxis]))
    non_constant_regressors = X[:, :-1] if X.shape[1] > 1 else np.array([])
    betas = np.linalg.pinv(X).dot(data)
    if remove_mean:
        datahat = X.dot(betas)
    else:  # disregard the first layer of X, which is degree 0
        datahat = X[:, 1:].dot(betas[1:, ...])
    regressed_data = data - datahat
    return regressed_data, non_constant_regressors
715730
def _distance_covariance_sqr_naive(x, y, exponent=1):
    
    a = _distance_matrix(x, exponent=exponent)
    b = _distance_matrix(y, exponent=exponent)

    return mean_product(a, b)
459067
def _Aterm(N, h, m, k, dW):
    """"""
    sqrt2h = np.sqrt(2.0/h)
    Xk = np.random.normal(0.0, 1.0, (N, m, 1))
    Yk = np.random.normal(0.0, 1.0, (N, m, 1))
    term1 = _dot(Xk, _t(Yk + sqrt2h*dW))
    term2 = _dot(Yk + sqrt2h*dW, _t(Xk))
    return (term1 - term2)/k
1366201
def _twofilter_smoothing_ON2(self, t, ti, info, phi, lwinfo):
        
        sp, sw = 0., 0.
        upb = lwinfo.max() + self.wgt[t].lw.max()
        if hasattr(self.model, 'upper_bound_trans'):
            upb += self.model.upper_bound_trans(t + 1)
        # Loop over n, to avoid having in memory a NxN matrix
        for n in range(self.N):
            omegan = np.exp(lwinfo + self.wgt[t].lw[n] - upb
                            + self.model.logpt(t + 1, self.X[t][n],
                                               info.hist.X[ti]))
            sp += np.sum(omegan * phi(self.X[t][n], info.hist.X[ti]))
            sw += np.sum(omegan)
        return sp / sw
314745
def fit(self,Params0=None,grad_threshold=1e-2):
        

        # GPVD initialization
        lik = limix.CLikNormalNULL()
        # Initial Params
        if Params0==None:
            n_params = self.C1.getNumberParams()
            n_params+= self.C2.getNumberParams()
            Params0 = SP.rand(n_params)
        # MultiGP Framework
        covar = []
        gp    = []
        mean  = []
        for i in range(self.N):
            covar.append(limix.CLinCombCF())
            covar[i].addCovariance(self.C1)
            covar[i].addCovariance(self.C2)
            coeff = SP.array([self.eigen[i],1])
            covar[i].setCoeff(coeff)
            mean.append(limix.CLinearMean(self.Yt[:,i],SP.eye(self.P)))
        gpVD = limix.CGPvarDecomp(covar[0],lik,mean[0],SP.ones(self.N),self.P,self.Yt,Params0)
        for i in range(self.N):
            gp.append(limix.CGPbase(covar[i],lik,mean[i]))
            gp[i].setY(self.Yt[:,i])
            gpVD.addGP(gp[i])
                
        # Optimization
        gpVD.initGPs()
        gpopt = limix.CGPopt(gpVD)
        LML0=-1.0*gpVD.LML()
        start_time = time.time()
        conv = gpopt.opt()
        time_train = time.time() - start_time
        LML=-1.0*gpVD.LML()
        LMLgrad = SP.linalg.norm(gpVD.LMLgrad()['covar'])
        Params = gpVD.getParams()['covar']
    
        # Check whether limix::CVarianceDecomposition.train() has converged
        if conv!=True or LMLgrad>grad_threshold or Params.max()>10:
            print('limix::CVarianceDecomposition::train has not converged')
            res=None
        else:
            res = {
                'Params0':          Params0,
                'Params':           Params,
                'LML':              SP.array([LML]),
                'LML0':             SP.array([LML0]),
                'LMLgrad':          SP.array([LMLgrad]),
                'time_train':       SP.array([time_train]),
                }
        return res
        pass
956087
def gradient_log_joint(self, x):
        
        T, D = self.T, self.D_latent
        assert x.shape == (T, D)

        # Collect the Gaussian LDS prior terms
        _, h_init, _ = self.info_init_params
        _, _, _, h1, h2, _ = self.info_dynamics_params
        H_diag, H_upper_diag = self.sparse_J_prior

        # Compute the gradient from the prior
        g = -1 * symm_block_tridiag_matmul(H_diag, H_upper_diag, x)
        g[0] += h_init
        g[:-1] += h1
        g[1:] += h2

        # Compute gradient from the likelihood terms
        g += self.grad_local_log_likelihood(x)

        return g
1088259
def model1(inps, lims, logger=None):
    
    a_arr, c_ = inps
    ilim, jlim = lims

    i_bs = sympy.symbols('i_lb i_ub', integer=True)
    i = sympy.Idx('i', i_bs)

    j_bs = sympy.symbols('j_lb j_ub', integer=True)
    j = sympy.Idx('j', j_bs)

    a = sympy.IndexedBase('a')

    c = sympy.Symbol('c', real=True)

    x = sympy.IndexedBase('x')
    y = sympy.IndexedBase('y')

    eqs = [
        sympy.Eq(x[i], (a[i]/3-1)**i+c),
        sympy.Eq(y[j], a[j]-j),
    ]

    ex_code = ExampleCode(eqs, (a[i], c), (i, j),
                          logger=logger, save_temp=True)
    x_, y_ = ex_code(inps, bounds=(ilim, jlim))
    x_ref = (a_arr/3-1)**np.arange(ilim[0], ilim[1]) + c_
    y_ref = a_arr[jlim[0]:jlim[1]] - np.arange(jlim[0], jlim[1])
    assert np.allclose(x_, x_ref)
    assert np.allclose(y_, y_ref)