Extract priors from a likelihood object.
517658
def get_priors(like):
    """"""

    npar = len(like.params())

    vals = np.ones(npar)
    errs = np.ones(npar)
    has_prior = np.array([False] * npar)

    for i, p in enumerate(like.params()):

        prior = like[i].log_prior()

        if prior is None:
            continue

        par_names = pyLike.StringVector()
        prior.getParamNames(par_names)

        if not 'Mean' in par_names:
            raise Exception('Failed to find Mean in prior parameters.')

        if not 'Sigma' in par_names:
            raise Exception('Failed to find Sigma in prior parameters.')

        for t in par_names:

            if t == 'Mean':
                vals[i] = prior.parameter(t).getValue()

            if t == 'Sigma':
                errs[i] = prior.parameter(t).getValue()

        has_prior[i] = True

    return vals, errs, has_prior
559142
def posterior_to_xarray(self):
        """"""
        posterior = self.posterior
        posterior_model = self.posterior_model
        # filter posterior_predictive and log_likelihood
        posterior_predictive = self.posterior_predictive
        if posterior_predictive is None:
            posterior_predictive = []
        elif isinstance(posterior_predictive, str):
            posterior_predictive = [posterior_predictive]
        log_likelihood = self.log_likelihood
        if not isinstance(log_likelihood, str):
            log_likelihood = []
        else:
            log_likelihood = [log_likelihood]

        ignore = posterior_predictive + log_likelihood

        data = get_draws_stan3(posterior, model=posterior_model, ignore=ignore)

        return dict_to_dataset(data, library=self.stan, coords=self.coords, dims=self.dims)
28703
def sample_prior(batch_size):
    cat, _ = get_distributions(DIST_PRIOR_PARAM[:NUM_CLASS], DIST_PRIOR_PARAM[NUM_CLASS:])
    sample_cat = tf.one_hot(cat.sample(batch_size), NUM_CLASS)

    
    sample_uni = tf.random_uniform([batch_size, NUM_UNIFORM], -1, 1)
    samples = tf.concat([sample_cat, sample_uni], axis=1)
    return samples
370801
def _expectation(self, X):
        """"""

        # Calculate probabilities of X belonging to the different clusters
        weighted_likelihoods = self._get_likelihoods(X) * self.priors
        sum_likelihoods = np.expand_dims(np.sum(weighted_likelihoods, axis=1), axis=1)
        # Determine responsibility as P(X|y)*P(y)/P(X)
        self.responsibility = weighted_likelihoods / sum_likelihoods
        # Assign samples to cluster that has largest probability
        self.sample_assignments = self.responsibility.argmax(axis=1)
        # Save value for convergence check
        self.responsibilities.append(np.max(self.responsibility, axis=1))
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
619129
def inference(self, kern, X, likelihood, Y, Y_metadata=None):

        
        N = X.shape[0] #number of training points
        D = X.shape[1] #number of dimensions

        Kds = np.zeros(D, dtype=object) #vector for holding covariance per dimension
        Qs = np.zeros(D, dtype=object) #vector for holding eigenvectors of covariance per dimension
        QTs = np.zeros(D, dtype=object) #vector for holding transposed eigenvectors of covariance per dimension
        V_kron = 1 # kronecker product of eigenvalues

        # retrieve the one-dimensional variation of the designated kernel
        oneDkernel = kern.get_one_dimensional_kernel(D)

        for d in range(D):
            xg = list(set(X[:,d])) #extract unique values for a dimension
            xg = np.reshape(xg, (len(xg), 1))
            oneDkernel.lengthscale = kern.lengthscale[d]
            Kds[d] = oneDkernel.K(xg)
            [V, Q] = np.linalg.eig(Kds[d])
            V_kron = np.kron(V_kron, V)
            Qs[d] = Q
            QTs[d] = Q.T

        noise = likelihood.variance + 1e-8

        alpha_kron = self.kron_mvprod(QTs, Y)
        V_kron = V_kron.reshape(-1, 1)
        alpha_kron = alpha_kron / (V_kron + noise)
        alpha_kron = self.kron_mvprod(Qs, alpha_kron)

        log_likelihood = -0.5 * (np.dot(Y.T, alpha_kron) + np.sum((np.log(V_kron + noise))) + N*log_2_pi)

        # compute derivatives wrt parameters Thete
        derivs = np.zeros(D+2, dtype='object')
        for t in range(len(derivs)):
            dKd_dTheta = np.zeros(D, dtype='object')
            gamma = np.zeros(D, dtype='object')
            gam = 1
            for d in range(D):
                xg = list(set(X[:,d]))
                xg = np.reshape(xg, (len(xg), 1))
                oneDkernel.lengthscale = kern.lengthscale[d]
                if t < D:
                    dKd_dTheta[d] = oneDkernel.dKd_dLen(xg, (t==d), lengthscale=kern.lengthscale[t]) #derivative wrt lengthscale
                elif (t == D):
                    dKd_dTheta[d] = oneDkernel.dKd_dVar(xg) #derivative wrt variance
                else:
                    dKd_dTheta[d] = np.identity(len(xg)) #derivative wrt noise
                gamma[d] = np.diag(np.dot(np.dot(QTs[d], dKd_dTheta[d].T), Qs[d]))
                gam = np.kron(gam, gamma[d])
            
            gam = gam.reshape(-1,1)
            kappa = self.kron_mvprod(dKd_dTheta, alpha_kron)
            derivs[t] = 0.5*np.dot(alpha_kron.T,kappa) - 0.5*np.sum(gam / (V_kron + noise))

        # separate derivatives
        dL_dLen = derivs[:D]
        dL_dVar = derivs[D]
        dL_dThetaL = derivs[D+1]

        return GridPosterior(alpha_kron=alpha_kron, QTs=QTs, Qs=Qs, V_kron=V_kron), \
                log_likelihood, {'dL_dLen':dL_dLen, 'dL_dVar':dL_dVar, 'dL_dthetaL':dL_dThetaL}
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
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)
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)
562419
def posterior_to_xarray(self):
        """"""
        data = self.posterior
        if not isinstance(data, dict):
            raise TypeError("DictConverter.posterior is not a dictionary")

        if "log_likelihood" in data:
            warnings.warn(
                "log_likelihood found in posterior."
                " For stats functions log_likelihood needs to be in sample_stats.",
                SyntaxWarning,
            )

        return dict_to_dataset(data, library=None, coords=self.coords, dims=self.dims)
955951
def map_parameter_estimate(self):
        
        res = self.pst.res
        assert res is not None
        # build the prior expectation parameter vector
        prior_expt = self.pst.parameter_data.loc[:,["parval1"]].copy()
        islog = self.pst.parameter_data.partrans == "log"
        prior_expt.loc[islog] = prior_expt.loc[islog].apply(np.log10)
        prior_expt = Matrix.from_dataframe(prior_expt)
        prior_expt.col_names = ["prior_expt"]
        # build the residual vector
        res_vec = Matrix.from_dataframe(res.loc[:,["residual"]])

        # form the terms of Schur's complement
        b = self.parcov * self.jco.T
        c = ((self.jco * self.parcov * self.jco.T) + self.obscov).inv
        bc = Matrix((b * c).x, row_names=b.row_names, col_names=c.col_names)

        # calc posterior expectation
        upgrade = bc * res_vec
        upgrade.col_names = ["prior_expt"]
        post_expt = prior_expt + upgrade

        # post processing - back log transform
        post_expt = pd.DataFrame(data=post_expt.x,index=post_expt.row_names,
                                 columns=["post_expt"])
        post_expt.loc[:,"prior_expt"] = prior_expt.x.flatten()
        post_expt.loc[islog,:] = 10.0**post_expt.loc[islog,:]
        # unc_sum = self.get_parameter_summary()
        # post_expt.loc[:,"standard_deviation"] = unc_sum.post_var.apply(np.sqrt)
        post_expt.sort_index(inplace=True)
        return post_expt
1267629
def log_uniform_prior(value, umin=0, umax=None):
    
    if value > 0 and value >= umin:
        if umax is not None:
            if value <= umax:
                return 1 / value
            else:
                return -np.inf
        else:
            return 1 / value
    else:
        return -np.inf
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])
1308294
def posterior(self, x, sigma=1.):
        """"""
        pr0 = 1. / self.scale**2  # prior precision
        prd = x.size / sigma**2  # data precision
        varp = 1. / (pr0 + prd)  # posterior variance
        mu = varp * (pr0 * self.loc + prd * x.mean())
        return Normal(loc=mu, scale=np.sqrt(varp))
422359
def prior_predictive_to_xarray(self):
        """"""
        prior = self.prior
        prior_predictive = self.prior_predictive
        data = get_draws(prior, variables=prior_predictive)
        return dict_to_dataset(data, library=self.pystan, coords=self.coords, dims=self.dims)
765754
def param_labels(self):
        """"""

        paramnames_labels = []
        prior_class_dict = self.variable.prior_class_dict
        prior_prior_model_dict = self.variable.prior_prior_model_dict

        for prior_name, prior in self.variable.prior_tuples_ordered_by_id:
            param_string = self.label_config.label(prior_name)
            prior_model = prior_prior_model_dict[prior]
            cls = prior_class_dict[prior]
            cls_string = "{}{}".format(self.label_config.subscript(cls), prior_model.component_number + 1)
            param_label = "{}_{{\\mathrm{{{}}}}}".format(param_string, cls_string)
            paramnames_labels.append(param_label)

        return paramnames_labels
554548
def _build_likelihood(self):
        

        """
        K = self.kern.K(self.X)
        L = tf.cholesky(
            K + tf.eye(tf.shape(self.X)[0], dtype=settings.float_type) * settings.numerics.jitter_level)
        F = tf.matmul(L, self.V) + self.mean_function(self.X)

        return tf.reduce_sum(self.likelihood.logp(F, self.Y))