Constructs a Wishart prior for the covariance matrix
420538
def construct_wishart(self,v,X):
        
        self.adjust_prior(list(range(int((len(self.latent_variables.z_list)-self.ylen-(self.ylen**2-self.ylen)/2)),
            int(len(self.latent_variables.z_list)))), fam.InverseWishart(v,X))
894685
def _setup_prior(self):
        
        if len(self._marg_prior) == 0:
            raise ValueError("A prior must be specified for the parameters "
                             "that you wish to marginalize the likelihood "
                             "over")
        marg_number = len([i for i in self._args if i != 0])
        if len(self._marg_prior) != marg_number:
            raise AttributeError("There is not a prior for each keyword "
                                 "argument")
        if self._margdist:
            bounds = self._marg_prior["distance"].bounds
            self._dist_array = numpy.linspace(bounds["distance"].min,
                                              bounds["distance"].max, 10**4)
            self._deltad = self._dist_array[1] - self._dist_array[0]
            self.dist_prior = numpy.array(
                                  [self._marg_prior["distance"].pdf(distance=i)
                                   for i in self._dist_array])
        if self._margtime:
            bounds = self._marg_prior["time"].bounds
            self._time_array = numpy.linspace(bounds["time"].min,
                                              bounds["time"].min, 10**4)
            self.time_prior = numpy.array(
                                  [self._marg_prior["time"].pdf(time=i) for
                                   i in self._time_array])
        if self._margphase:
            bounds = self._marg_prior["phase"].bounds
            self._phase_array = numpy.linspace(bounds["phase"].min,
                                               bounds["phase"].max, 10**4)
            self._deltap = self._phase_array[1] - self._phase_array[0]
            self.phase_prior = numpy.array(
                                   [self._marg_prior["phase"].pdf(phase=i) for
                                    i in self._phase_array])
109209
def posterior_covariance_between_points(self, X1, X2):
        
        return self.posterior.covariance_between_points(self.kern, self.X, X1, X2)
578443
def error_perturbation_var(C, S):
    

    """
    K = S.shape[0]
    cov = tmatrix_cov(C)
    for i in range(K):
        R = S[i, :, :]
        X[i] = (R[:, :, np.newaxis] * cov * R[:, np.newaxis, :]).sum()
969104
def _posterior(self, x):
        
        yedge = self._nuis_pdf.marginalization_bins()
        yc = 0.5 * (yedge[1:] + yedge[:-1])
        yw = yedge[1:] - yedge[:-1]

        like_array = self.like(x[:, np.newaxis], yc[np.newaxis, :]) * yw
        like_array /= like_array.sum()

        self._post = like_array.sum(1)
        self._post_interp = castro.Interpolator(x, self._post)
        return self._post
44949
def aligned_covariance(fit, type='noise'):
    
    cov = fit._covariance_matrix(type)
    # Rescale eigenvectors to sum to 1
    cov /= N.linalg.norm(cov)
    return dot(fit.axes,cov)
774414
def _build_likelihood(self):
        

        """

        # Get prior KL.
        KL = gauss_kl(self.q_mu, self.q_sqrt)

        # Get conditionals
        K = self.kern.K(self.X) + tf.eye(self.num_data, dtype=settings.float_type) * \
            settings.numerics.jitter_level
        L = tf.cholesky(K)

        fmean = tf.matmul(L, self.q_mu) + self.mean_function(self.X)  # NN,ND->ND

        q_sqrt_dnn = tf.matrix_band_part(self.q_sqrt, -1, 0)  # D x N x N

        L_tiled = tf.tile(tf.expand_dims(L, 0), tf.stack([self.num_latent, 1, 1]))

        LTA = tf.matmul(L_tiled, q_sqrt_dnn)  # D x N x N
        fvar = tf.reduce_sum(tf.square(LTA), 2)

        fvar = tf.transpose(fvar)

        # Get variational expectations.
        var_exp = self.likelihood.variational_expectations(fmean, fvar, self.Y)

        return tf.reduce_sum(var_exp) - KL
1117737
def covariance_between_points(self, kern, X, X1, X2):
        
        # ndim == 3 is a model for missing data
        if self.woodbury_chol.ndim != 2:
            raise RuntimeError("This method does not support posterior for missing data models")

        Kx1 = kern.K(X, X1)
        Kx2 = kern.K(X, X2)
        K12 = kern.K(X1, X2)

        tmp1 = dtrtrs(self.woodbury_chol, Kx1)[0]
        tmp2 = dtrtrs(self.woodbury_chol, Kx2)[0]
        var = K12 - tmp1.T.dot(tmp2)

        return var
1102802
def cov2corr(cov):
    
    assert cov.ndim == 2, 'covariance matrix should be 2D array'
    inv_sd = 1 / np.sqrt(np.diag(cov))
    corr = cov * inv_sd[None, :] * inv_sd[:, None]
    return corr
1153674
def _maximization_step(X, posteriors):
    

    ### Prior probabilities or class weights
    sum_post_proba = np.sum(posteriors, axis=0)
    prior_proba = sum_post_proba / (sum_post_proba.sum() + Epsilon)
    
    ### means
    means = np.dot(posteriors.T, X) / (sum_post_proba[:, np.newaxis] + Epsilon)

    ### covariance matrices
    n_components = posteriors.shape[1]
    n_features = X.shape[1]
    covars = np.empty(shape=(n_components, n_features, n_features), dtype=float)
  
    for i in range(n_components):
        post_i = posteriors[:, i]
        mean_i = means[i]
        diff_i = X - mean_i

        with np.errstate(under='ignore'):
            covar_i = np.dot(post_i * diff_i.T, diff_i) / (post_i.sum() + Epsilon)
        covars[i] = covar_i + Lambda * np.eye(n_features)


    _validate_params(prior_proba, means, covars)
    return(prior_proba, means, covars)
648627
def normal_prior(value, mean, sigma):
    
    return -0.5 * (2 * np.pi * sigma) - (value - mean) ** 2 / (2.0 * sigma)
564733
def show_corr(self):
        ""
        from .ndshow_gtk3 import view
        d = np.diag(self.covar) ** -0.5
        corr = self.covar * d[np.newaxis,:] * d[:,np.newaxis]
        view(corr, title='Correlation Matrix')
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
917828
def inference_minibatch(self, kern, X, Z, likelihood, Y, Zp):
        

        num_data, output_dim = Y.shape

        #see whether we've got a different noise variance for each datum
        beta = 1./np.fmax(likelihood.variance, 1e-6)
        het_noise = beta.size > 1
        # VVT_factor is a matrix such that tdot(VVT_factor) = VVT...this is for efficiency!
        #self.YYTfactor = beta*self.get_YYTfactor(Y)
        if self.Y_speedup and not het_noise:
            YYT_factor = self.get_YYTfactor(Y)
        else:
            YYT_factor = Y

        n_start = self.batch_pos
        batchsize = num_data if self.batchsize is None else self.batchsize
        n_end = min(batchsize+n_start, num_data)
        if n_end == num_data:
            isEnd = True
            self.batch_pos = 0
        else:
            isEnd = False
            self.batch_pos = n_end

        if batchsize == num_data:
            Y_slice = YYT_factor
            X_slice = X
        else:
            Y_slice = YYT_factor[n_start:n_end]
            X_slice = X[n_start:n_end]

        psi0 = kern._Kdiag(X_slice) #Kffdiag
        psi1 = kern.K(X_slice, Z) #Kfu
        betapsi1 = np.einsum('n,nm->nm', beta, psi1)

        X_slice = X_slice.values
        Z = Z.values

        Zp = Zp.gamma.values
        indX = np.int_(X_slice[:, -1])
        indZ = np.int_(Z[:, -1]) - Zp.shape[0]

        betaY = beta*Y_slice

        #======================================================================
        # Load Intermediate Results
        #======================================================================

        dL_dpsi2R = self.midRes['dL_dpsi2R']
        v = self.midRes['v']

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

        dL_dpsi0 = -.5*output_dim*(beta * Zp[indX, :]) #XxQ #TODO: Check this gradient

        dL_dpsi1 = np.dot(betaY, v.T)
        dL_dEZp  = psi1*dL_dpsi1
        dL_dpsi1 = Zp[np.ix_(indX, indZ)]*dL_dpsi1
        dL_dgamma = np.zeros(Zp.shape)
        for d in np.unique(indX):
            indd = indX == d
            betapsi1d = betapsi1[indd, :]
            psi1d = psi1[indd, :]
            Zpd = Zp[d, :]
            Zp2 = Zpd[:, None]*Zpd[None, :] - np.diag(np.power(Zpd, 2)) + np.diag(Zpd)
            dL_dpsi1[indd, :] += np.dot(betapsi1d, Zp2[np.ix_(indZ, indZ)] * dL_dpsi2R)*2.

            dL_EZp2 = dL_dpsi2R * (np.dot(psi1d.T, psi1d) * beta)*2.  # Zpd*Kufd*Kfud*beta
            #Gradient of Likelihood wrt gamma is calculated here
            EZ = Zp[d, indZ]
            for q in range(Zp.shape[1]):
                EZt = EZ.copy()
                indq = indZ == q
                EZt[indq] = .5
                dL_dgamma[d, q] = np.sum(dL_dEZp[np.ix_(indd, indq)]) + np.sum(dL_EZp2[:, indq]*EZt[:, None]) -\
                    .5*beta*(np.sum(psi0[indd, q]))

        #======================================================================
        # Compute dL_dthetaL
        #======================================================================
        if isEnd:
            dL_dthetaL = self.midRes['dL_dthetaL']
        else:
            dL_dthetaL = 0.

        grad_dict = {'dL_dKdiag': dL_dpsi0,
                     'dL_dKnm': dL_dpsi1,
                     'dL_dthetaL': dL_dthetaL,
                     'dL_dgamma': dL_dgamma}

        return isEnd, (n_start, n_end), grad_dict
138527
def log_prior(self):
        """"""
        for p, b in zip(self.parameter_vector, self.parameter_bounds):
            if b[0] is not None and p < b[0]:
                return -np.inf
            if b[1] is not None and p > b[1]:
                return -np.inf
        return 0.0
811708
def _create_model(self, X, Y):
        

        # --- define kernel
        self.input_dim = X.shape[1]
        if self.kernel is None:
            kern = GPy.kern.RBF(self.input_dim, variance=1.)
        else:
            kern = self.kernel
            self.kernel = None

        # --- define model
        noise_var = Y.var()*0.01 if self.noise_var is None else self.noise_var
        self.model = GPy.models.GPRegression(X, Y, kernel=kern, noise_var=noise_var)

        # --- Define prior on the hyper-parameters for the kernel (for integrated acquisitions)
        self.model.kern.set_prior(GPy.priors.Gamma.from_EV(2.,4.))
        self.model.likelihood.variance.set_prior(GPy.priors.Gamma.from_EV(2.,4.))

        # --- Restrict variance if exact evaluations of the objective
        if self.exact_feval:
            self.model.Gaussian_noise.constrain_fixed(1e-6, warning=False)
        else:
            self.model.Gaussian_noise.constrain_positive(warning=False)
155020
def fit(self, X):
        
        params_dict = _fit_gmm_params(X=X, n_mixtures=self.n_clusters, \
                        n_init=self.n_trials, init_method=self.init_method, \
                        n_iter=self.max_iter, tol=self.tol)
        self.priors_ = params_dict['priors']
        self.means_  = params_dict['means']
        self.covars_ = params_dict['covars']

        self.converged = True
        self.labels_ = self.predict(X)