Creates likelihood markov blanket of the model

        Parameters
        ----------
        beta : np.array
            Contains untransformed starting values for latent variables

        Returns
        ----------
        - Negative loglikelihood
542664
def likelihood_markov_blanket(self, beta):
             
        states = beta[self.z_no:self.z_no+self.data_length] # the local level (untransformed)
        parm = np.array([self.latent_variables.z_list[k].prior.transform(beta[k]) for k in range(self.z_no)]) # transformed distribution parameters
        scale, shape, skewness = self._get_scale_and_shape(parm)
        return self.family.markov_blanket(self.data, self.link(states), scale, shape, skewness)
275293
def inference(self, kern, X, Z, likelihood, Y, indexD, output_dim, Y_metadata=None, Lm=None, dL_dKmm=None, Kuu_sigma=None):
        

        input_dim = Z.shape[0]

        uncertain_inputs = isinstance(X, VariationalPosterior)

        beta = 1./likelihood.variance
        if len(beta)==1:
            beta = np.zeros(output_dim)+beta

        beta_exp = np.zeros(indexD.shape[0])
        for d in range(output_dim):
            beta_exp[indexD==d] = beta[d]

        psi0, psi1, psi2 = self.gatherPsiStat(kern, X, Z, Y, beta, uncertain_inputs)

        psi2_sum = (beta_exp[:,None,None]*psi2).sum(0)/output_dim

        #======================================================================
        # Compute Common Components
        #======================================================================

        Kmm = kern.K(Z).copy()
        if Kuu_sigma is not None:
            diag.add(Kmm, Kuu_sigma)
        else:
            diag.add(Kmm, self.const_jitter)
        Lm = jitchol(Kmm)

        logL = 0.
        dL_dthetaL = np.zeros(output_dim)
        dL_dKmm = np.zeros_like(Kmm)
        dL_dpsi0 = np.zeros_like(psi0)
        dL_dpsi1 = np.zeros_like(psi1)
        dL_dpsi2 = np.zeros_like(psi2)
        wv = np.empty((Kmm.shape[0],output_dim))

        for d in range(output_dim):
            idx_d = indexD==d
            Y_d = Y[idx_d]
            N_d = Y_d.shape[0]
            beta_d = beta[d]

            psi2_d = psi2[idx_d].sum(0)*beta_d
            psi1Y = Y_d.T.dot(psi1[idx_d])*beta_d
            psi0_d = psi0[idx_d].sum()*beta_d
            YRY_d = np.square(Y_d).sum()*beta_d

            LmInvPsi2LmInvT = backsub_both_sides(Lm, psi2_d, 'right')

            Lambda = np.eye(Kmm.shape[0])+LmInvPsi2LmInvT
            LL = jitchol(Lambda)
            LmLL = Lm.dot(LL)

            b  = dtrtrs(LmLL, psi1Y.T)[0].T
            bbt = np.square(b).sum()
            v = dtrtrs(LmLL, b.T, trans=1)[0].T
            LLinvPsi1TYYTPsi1LLinvT = tdot(b.T)

            tmp = -backsub_both_sides(LL, LLinvPsi1TYYTPsi1LLinvT)
            dL_dpsi2R = backsub_both_sides(Lm, tmp+np.eye(input_dim))/2

            logL_R = -N_d*np.log(beta_d)
            logL += -((N_d*log_2_pi+logL_R+psi0_d-np.trace(LmInvPsi2LmInvT))+YRY_d- bbt)/2.

            dL_dKmm +=  dL_dpsi2R - backsub_both_sides(Lm, LmInvPsi2LmInvT)/2

            dL_dthetaL[d:d+1] = (YRY_d*beta_d + beta_d*psi0_d - N_d*beta_d)/2. - beta_d*(dL_dpsi2R*psi2_d).sum() - beta_d*np.trace(LLinvPsi1TYYTPsi1LLinvT)

            dL_dpsi0[idx_d] = -beta_d/2.
            dL_dpsi1[idx_d] = beta_d*np.dot(Y_d,v)
            dL_dpsi2[idx_d] = beta_d*dL_dpsi2R
            wv[:,d] = v

        LmInvPsi2LmInvT = backsub_both_sides(Lm, psi2_sum, 'right')

        Lambda = np.eye(Kmm.shape[0])+LmInvPsi2LmInvT
        LL = jitchol(Lambda)
        LmLL = Lm.dot(LL)
        logdet_L = 2.*np.sum(np.log(np.diag(LL)))
        dL_dpsi2R_common = dpotri(LmLL)[0]/-2.
        dL_dpsi2 += dL_dpsi2R_common[None,:,:]*beta_exp[:,None,None]

        for d in range(output_dim):
            dL_dthetaL[d] += (dL_dpsi2R_common*psi2[indexD==d].sum(0)).sum()*-beta[d]*beta[d]

        dL_dKmm += dL_dpsi2R_common*output_dim

        logL += -output_dim*logdet_L/2.

        #======================================================================
        # Compute dL_dKmm
        #======================================================================

        # dL_dKmm =  dL_dpsi2R - output_dim* backsub_both_sides(Lm, LmInvPsi2LmInvT)/2 #LmInv.T.dot(LmInvPsi2LmInvT).dot(LmInv)/2.

        #======================================================================
        # Compute the Posterior distribution of inducing points p(u|Y)
        #======================================================================

        LLInvLmT = dtrtrs(LL, Lm.T)[0]
        cov = tdot(LLInvLmT.T)

        wd_inv = backsub_both_sides(Lm, np.eye(input_dim)- backsub_both_sides(LL, np.identity(input_dim), transpose='left'), transpose='left')
        post = Posterior(woodbury_inv=wd_inv, woodbury_vector=wv, K=Kmm, mean=None, cov=cov, K_chol=Lm)

        #======================================================================
        # Compute dL_dthetaL for uncertian input and non-heter noise
        #======================================================================

        # for d in range(output_dim):
        #     dL_dthetaL[d:d+1] += - beta[d]*beta[d]*(dL_dpsi2R[None,:,:] * psi2[indexD==d]/output_dim).sum()
        # dL_dthetaL += - (dL_dpsi2R[None,:,:] * psi2_sum*D beta*(dL_dpsi2R*psi2).sum()

        #======================================================================
        # Compute dL_dpsi
        #======================================================================

        if not uncertain_inputs:
            dL_dpsi1 += (psi1[:,None,:]*dL_dpsi2).sum(2)*2.

        if uncertain_inputs:
            grad_dict = {'dL_dKmm': dL_dKmm,
                         'dL_dpsi0':dL_dpsi0,
                         'dL_dpsi1':dL_dpsi1,
                         'dL_dpsi2':dL_dpsi2,
                         'dL_dthetaL':dL_dthetaL}
        else:
            grad_dict = {'dL_dKmm': dL_dKmm,
                         'dL_dKdiag':dL_dpsi0,
                         'dL_dKnm':dL_dpsi1,
                         'dL_dthetaL':dL_dthetaL}

        return post, logL, grad_dict
818263
def graph_lasso(X, num_folds):
    
    print("GraphLasso (sklearn)")
    model = GraphLassoCV(cv=num_folds)
    model.fit(X)
    print("   lam_: {}".format(model.alpha_))
    return model.covariance_, model.precision_, model.alpha_
985675
def _init_parameters_random(self, X_bin):
        

        _, n_features = X_bin.shape

        # The parameter class_log_prior_ has shape (2,). The values represent
        # 'match' and 'non-match'.
        rand_vals = np.random.rand(2)
        class_prior = rand_vals / np.sum(rand_vals)

        # make empty array of feature log probs
        # dimensions 2xn_features
        feature_prob = np.zeros((2, n_features))

        feat_i = 0

        for i, bin in enumerate(self._binarizers):

            bin_len = bin.classes_.shape[0]

            rand_vals_0 = np.random.rand(bin_len)
            feature_prob[0, feat_i:feat_i + bin_len] = \
                rand_vals_0 / np.sum(rand_vals_0)

            rand_vals_1 = np.random.rand(bin_len)
            feature_prob[1, feat_i:feat_i + bin_len] = \
                rand_vals_1 / np.sum(rand_vals_1)

            feat_i += bin_len

        return np.log(class_prior), np.log(feature_prob)
578920
def predict_proba(self, L):
        
        n = L.shape[0]
        Y_p = np.random.rand(n, self.k)
        Y_p /= Y_p.sum(axis=1).reshape(-1, 1)
        return Y_p
1190905
def lnprior(pars):
    
    # Limit norm and B to be positive
    logprob = (
        naima.uniform_prior(pars[0], 0.0, np.inf)
        + naima.uniform_prior(pars[1], -1, 5)
        + naima.uniform_prior(pars[3], 0, np.inf)
    )

    return logprob
697712
def _create_latent_variables(self):
        

        for p_term in range(self.p):
            self.latent_variables.add_z('p(' + str(p_term+1) + ')', fam.Normal(0,0.5,transform='logit'), fam.Normal(0,3))
            if p_term == 0:
                self.latent_variables.z_list[-1].start = 3.00
            else:
                self.latent_variables.z_list[-1].start = -4.00

        for q_term in range(self.q):
            self.latent_variables.add_z('q(' + str(q_term+1) + ')', fam.Normal(0,0.5,transform='logit'), fam.Normal(0,3))
            if q_term == 0:
                self.latent_variables.z_list[-1].start = -1.50  
            else: 
                self.latent_variables.z_list[-1].start = -4.00  

        self.latent_variables.add_z('v', fam.Flat(transform='exp'), fam.Normal(0,3))
        self.latent_variables.add_z('GARCH-M', fam.Normal(0, 3, transform=None),fam.Normal(0,3))

        for parm in range(len(self.X_names)):
            self.latent_variables.add_z('Vol Beta ' + self.X_names[parm], fam.Normal(0,10,transform=None), fam.Normal(0,3))

        for parm in range(len(self.X_names)):
            self.latent_variables.add_z('Returns Beta ' + self.X_names[parm], fam.Normal(0,10,transform=None), fam.Normal(0,3))

        # Starting values        

        for i in range(self.p+self.q, self.z_no):
            self.latent_variables.z_list[i].start = 0.0

        self.latent_variables.z_list[self.p+self.q].start = 2.0
204166
def rweibull(alpha, beta, size=None):
    
    tmp = -np.log(runiform(0, 1, size))
    return beta * (tmp ** (1. / alpha))
305490
def grad_local_log_likelihood(self, x):
        
        C, D, u, y = self.C, self.D, self.inputs, self.data
        psi = x.dot(C.T) + u.dot(D.T)
        p = 1. / (1 + np.exp(-psi))
        return (y - p).dot(C)
1293184
def _create_latent_variables(self):
        

        # Create latent variables
        for no, i in enumerate(self.kernel.build_latent_variables()):
            self.latent_variables.add_z(i[0],i[1],i[2])
            self.latent_variables.z_list[no].start = i[3]

        self.z_no = len(self.kernel.build_latent_variables())

        # Use an ARIMA model to find starting point for the initial noise latent variable
        arma_start = arma.ARIMA(self.data, ar=self.ar, ma=0, integ=self.integ)
        x = arma_start.fit()
        arma_starting_values = arma_start.latent_variables.get_z_values()
        self.latent_variables.z_list[0].start = np.log(np.exp(np.power(arma_starting_values[-1],2)))
791801
def log_likelihood_pairwise(data, params):
    """"""
    loglik = 0
    for winner, loser in data:
        loglik -= np.logaddexp(0, -(params[winner] - params[loser]))
    return loglik
1212095
def expected_hmm_logprob(pi_0, trans_matrix, stats):
    
    E_z, sum_E_ztztp1T, _ = stats
    T, K = E_z.shape
    assert sum_E_ztztp1T.shape == (K, K)

    out = 0
    out += np.dot(E_z[0], np.log(pi_0))
    out += np.sum(sum_E_ztztp1T * np.log(trans_matrix))
    return out
848441
def loglike(self):
        '''
        
        '''
        sum = logp_of_set(self.children)
        if self.verbose > 2:
            print_('\t' + self._id + ' Current log-likelihood ', sum)
        return sum
100501
def start_naive_bayes(automated_run, session, path):
    
    module = functions.import_string_code_as_module(automated_run.source)
    random_state = 8 if not hasattr(module, 'random_state') else module.random_state
    assert module.metric_to_optimize in automated_run.base_learner_origin.metric_generators

    # get non-searchable parameters
    base_estimator = automated_run.base_learner_origin.return_estimator()
    base_estimator.set_params(**module.default_params)
    default_params = functions.make_serializable(base_estimator.get_params())
    non_searchable_params = dict((key, val) for key, val in iteritems(default_params)
                                 if key not in module.pbounds)

    # get already calculated base learners in search space
    existing_base_learners = []
    for base_learner in automated_run.base_learner_origin.base_learners:
        if not base_learner.job_status == 'finished':
            continue
        in_search_space = True
        for key, val in iteritems(non_searchable_params):
            if base_learner.hyperparameters[key] != val:
                in_search_space = False
                break  # If no match, move on to the next base learner
        if in_search_space:
            existing_base_learners.append(base_learner)

    # build initialize dictionary
    target = []
    initialization_dict = dict((key, list()) for key in module.pbounds.keys())
    for base_learner in existing_base_learners:
        # check if base learner's searchable hyperparameters are all numerical
        all_numerical = True
        for key in module.pbounds.keys():
            if not isinstance(base_learner.hyperparameters[key], numbers.Number):
                all_numerical = False
                break
        if not all_numerical:
            continue  # if there is a non-numerical hyperparameter, skip this.

        for key in module.pbounds.keys():
            initialization_dict[key].append(base_learner.hyperparameters[key])
        target.append(base_learner.individual_score[module.metric_to_optimize])
    initialization_dict['target'] = target if not module.invert_metric \
        else list(map(lambda x: -x, target))
    print('{} existing in initialization dictionary'.
          format(len(initialization_dict['target'])))

    # Create function to be optimized
    func_to_optimize = return_func_to_optimize(
        path, session, automated_run.base_learner_origin, module.default_params,
        module.metric_to_optimize, module.invert_metric, set(module.integers)
    )

    # Create Bayes object
    bo = BayesianOptimization(func_to_optimize, module.pbounds)

    bo.initialize(initialization_dict)

    np.random.seed(random_state)

    bo.maximize(**module.maximize_config)
557383
def _build_likelihood(self):
        
        pX = DiagonalGaussian(self.X_mean, self.X_var)

        num_inducing = len(self.feature)
        psi0 = tf.reduce_sum(expectation(pX, self.kern))
        psi1 = expectation(pX, (self.kern, self.feature))
        psi2 = tf.reduce_sum(expectation(pX, (self.kern, self.feature), (self.kern, self.feature)), axis=0)
        Kuu = features.Kuu(self.feature, self.kern, jitter=settings.jitter)
        L = tf.cholesky(Kuu)
        sigma2 = self.likelihood.variance
        sigma = tf.sqrt(sigma2)

        # Compute intermediate matrices
        A = tf.matrix_triangular_solve(L, tf.transpose(psi1), lower=True) / sigma
        tmp = tf.matrix_triangular_solve(L, psi2, lower=True)
        AAT = tf.matrix_triangular_solve(L, tf.transpose(tmp), lower=True) / sigma2
        B = AAT + tf.eye(num_inducing, dtype=settings.float_type)
        LB = tf.cholesky(B)
        log_det_B = 2. * tf.reduce_sum(tf.log(tf.matrix_diag_part(LB)))
        c = tf.matrix_triangular_solve(LB, tf.matmul(A, self.Y), lower=True) / sigma

        # KL[q(x) || p(x)]
        dX_var = self.X_var if len(self.X_var.get_shape()) == 2 else tf.matrix_diag_part(self.X_var)
        NQ = tf.cast(tf.size(self.X_mean), settings.float_type)
        D = tf.cast(tf.shape(self.Y)[1], settings.float_type)
        KL = -0.5 * tf.reduce_sum(tf.log(dX_var)) \
             + 0.5 * tf.reduce_sum(tf.log(self.X_prior_var)) \
             - 0.5 * NQ \
             + 0.5 * tf.reduce_sum((tf.square(self.X_mean - self.X_prior_mean) + dX_var) / self.X_prior_var)

        # compute log marginal bound
        ND = tf.cast(tf.size(self.Y), settings.float_type)
        bound = -0.5 * ND * tf.log(2 * np.pi * sigma2)
        bound += -0.5 * D * log_det_B
        bound += -0.5 * tf.reduce_sum(tf.square(self.Y)) / sigma2
        bound += 0.5 * tf.reduce_sum(tf.square(c))
        bound += -0.5 * D * (tf.reduce_sum(psi0) / sigma2 -
                             tf.reduce_sum(tf.matrix_diag_part(AAT)))
        bound -= KL
        return bound
599227
def _create_likelihood(self, srcmdl=None):
        """"""

        self._like = SummedLikelihood()
        for c in self.components:
            c._create_binned_analysis(srcmdl)
            self._like.addComponent(c.like)

        self.like.model = self.like.components[0].model
        self._fitcache = None
        self._init_roi_model()
640624
def create_gamma_model(alignment, missing_data=None, ncat=4):
    
    model = alignment.parameters.partitions.model
    freqs = alignment.parameters.partitions.frequencies
    alpha = alignment.parameters.partitions.alpha
    if model == 'LG':
        subs_model = LG(freqs)
    elif model == 'WAG':
        subs_model = WAG(freqs)
    elif model == 'GTR':
        rates = alignment.parameters.partitions.rates
        subs_model = GTR(rates, freqs, True)
    else:
        raise ValueError("Can't handle this model: {}".format(model))
    tm = TransitionMatrix(subs_model)
    gamma = GammaMixture(alpha, ncat)
    gamma.init_models(tm, alignment_to_partials(alignment, missing_data))
    return gamma