log-prior for model parameters
488629
def lnprior(self, mass, age, feh,
                distance=None, AV=None, 
                use_local_fehprior=True):
        
        mass_prior = salpeter_prior(mass)
        if mass_prior==0:
            mass_lnprior = -np.inf
        else:
            mass_lnprior = np.log(mass_prior)

        if np.isnan(mass_lnprior):
            logging.warning('mass prior is nan at {}'.format(mass))

        age_lnprior = np.log(age * (2/(self.ic.maxage**2-self.ic.minage**2)))
        if np.isnan(age_lnprior):
            logging.warning('age prior is nan at {}'.format(age))


        if use_local_fehprior:
            fehdist = local_fehdist(feh)
        else:
            fehdist = 1/(self.ic.maxfeh - self.ic.minfeh)
        feh_lnprior = np.log(fehdist)
        if np.isnan(feh_lnprior):
            logging.warning('feh prior is nan at {}'.format(feh))


        if distance is not None:
            if distance <= 0:
                distance_lnprior = -np.inf
            else:
                distance_lnprior = np.log(3/self.max_distance**3 * distance**2)
        else:
            distance_lnprior = 0
        if np.isnan(distance_lnprior):
            logging.warning('distance prior is nan at {}'.format(distance))


        if AV is not None:
            AV_lnprior = np.log(1/self.maxAV)
        else:
            AV_lnprior = 0
        if np.isnan(AV_lnprior):
            logging.warning('AV prior is nan at {}'.format(AV))
            

        lnprior = (mass_lnprior + age_lnprior + feh_lnprior + 
                distance_lnprior + AV_lnprior)

        return lnprior
1254275
def grad_local_log_likelihood(self, x):
        
        T, D = self.T, self.D_latent
        assert x.shape == (T, D)
        gfun = grad(self.local_log_likelihood)

        g = np.zeros((T, D))
        for t in range(T):
            g[t] += gfun(x[t], self.data[t], self.inputs[t])
        return g
1151312
def sparse_hessian_log_joint(self, x):
        
        T, D = self.T, self.D_latent
        assert x.shape == (T, D)

        # Collect the Gaussian LDS prior terms
        J_diag, J_upper_diag = self.sparse_J_prior
        H_diag, H_upper_diag = -J_diag, -J_upper_diag

        # Collect the likelihood terms
        H_diag += self.hessian_local_log_likelihood(x)

        # Subtract a little bit to ensure negative definiteness
        H_diag -= 1e-8 * np.eye(D)

        return H_diag, H_upper_diag
419981
def inference_minibatch(self, kern, X, Z, likelihood, Y):
        

        num_data, output_dim = Y.shape

        if isinstance(X, VariationalPosterior):
            uncertain_inputs = True
        else:
            uncertain_inputs = False

        #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]

        if not uncertain_inputs:
            psi0 = kern.Kdiag(X_slice)
            psi1 = kern.K(X_slice, Z)
            psi2 = None
            betapsi1 = np.einsum('n,nm->nm',beta,psi1)
        elif het_noise:
            psi0 = kern.psi0(Z, X_slice)
            psi1 = kern.psi1(Z, X_slice)
            psi2 = kern.psi2(Z, X_slice)
            betapsi1 = np.einsum('n,nm->nm',beta,psi1)

        if het_noise:
            beta = beta[n_start] # assuming batchsize==1

        betaY = beta*Y_slice

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

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

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

        dL_dpsi0 = -output_dim * (beta * np.ones((n_end-n_start,)))/2.

        dL_dpsi1 = np.dot(betaY,v.T)

        if uncertain_inputs:
            dL_dpsi2 = beta* dL_dpsi2R
        else:
            dL_dpsi1 += np.dot(betapsi1,dL_dpsi2R)*2.
            dL_dpsi2 = None

        #======================================================================
        # Compute dL_dthetaL
        #======================================================================

        if het_noise:
            if uncertain_inputs:
                psiR = np.einsum('mo,mo->',dL_dpsi2R,psi2)
            else:
                psiR = np.einsum('nm,no,mo->',psi1,psi1,dL_dpsi2R)

            dL_dthetaL = ((np.square(betaY)).sum(axis=-1) + np.square(beta)*(output_dim*psi0)-output_dim*beta)/2. - np.square(beta)*psiR- (betaY*np.dot(betapsi1,v)).sum(axis=-1)
        else:
            if isEnd:
                dL_dthetaL = self.midRes['dL_dthetaL']
            else:
                dL_dthetaL = 0.

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

        return isEnd, (n_start,n_end), grad_dict
875462
def gaussian_tuple_prior_for_arguments(self, arguments):
        
        tuple_prior = TuplePrior()
        for prior_tuple in self.prior_tuples:
            setattr(tuple_prior, prior_tuple.name, arguments[prior_tuple.prior])
        return tuple_prior
156391
def map_forecast_estimate(self):
        
        assert self.forecasts is not None
        islog = self.pst.parameter_data.partrans == "log"
        par_map = self.map_parameter_estimate
        par_map.loc[islog,:] = np.log10(par_map.loc[islog,:])
        par_map = Matrix.from_dataframe(par_map.loc[:,["post_expt"]])
        posts,priors = [],[]
        post_expt = (self.predictions.T * par_map).to_dataframe()
        for fname in self.forecast_names:
            #fname = forecast.col_names[0]
            pr = self.pst.res.loc[fname,"modelled"]
            priors.append(pr)
            posts.append(pr + post_expt.loc[fname,"post_expt"])
        return pd.DataFrame(data=np.array([priors,posts]).transpose(),
                            columns=["prior_expt","post_expt"],
                            index=self.forecast_names)
851071
def data_log_likelihood(self, successes, trials, beta):
        ''''''
        return binom.logpmf(successes, trials, 1.0 / (1 + np.exp(-beta))).sum()
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
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
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)
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)
646451
def _build_likelihood(self):
        

        """
        K = self.kern.K(self.X) + tf.eye(tf.shape(self.X)[0], dtype=settings.float_type) * self.likelihood.variance
        L = tf.cholesky(K)
        m = self.mean_function(self.X)
        logpdf = multivariate_normal(self.Y, m, L)  # (R,) log-likelihoods for each independent dimension of Y

        return tf.reduce_sum(logpdf)
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)
1095810
def _margtimedist_loglr(self, mf_snr, opt_snr):
        
        logl = special.logsumexp(mf_snr, b=self._deltat)
        logl_marg = logl/self._dist_array
        opt_snr_marg = opt_snr/self._dist_array**2
        return special.logsumexp(logl_marg - 0.5*opt_snr_marg,
                                 b=self._deltad*self.dist_prior)
532389
def fgmc(log_fg_ratios, mu_log_vt, sigma_log_vt, Rf, maxfg):
    '''
    
    '''

    Lb = np.random.uniform(0., maxfg, len(Rf))
    pquit = 0

    while pquit < 0.1:
        # quit when the posterior on Lf is very close to its prior

        nsamp = len(Lb)
        Rf_sel = np.random.choice(Rf, nsamp)
        vt = np.random.lognormal(mu_log_vt, sigma_log_vt, len(Rf_sel))

        Lf = Rf_sel * vt

        log_Lf, log_Lb = log(Lf), log(Lb)

        plR = 0
        for lfr in log_fg_ratios:
            plR += np.logaddexp(lfr + log_Lf, log_Lb)

        plR -= (Lf + Lb)
        plRn = plR - max(plR)

        idx = np.exp(plRn) > np.random.random(len(plRn))

        pquit = ss.stats.ks_2samp(Lb, Lb[idx])[1]

        Lb = Lb[idx]

    return Rf_sel[idx], Lf[idx], Lb
1156744
def parameters_changed(self):
        
        # using the warped X to update
        self.X = self.transform_data(self.X_untransformed)
        super(InputWarpedGP, self).parameters_changed()
        # the gradient of log likelihood w.r.t. input AFTER warping is a product of dL_dK and dK_dX
        dL_dX = self.kern.gradients_X(self.grad_dict['dL_dK'], self.X)
        self.warping_function.update_grads(self.X_untransformed, dL_dX)
515429
private void findParameters(int iterations, double correctionConstant) {
        LOG.info("Performing max " + iterations + " iterations.");

        double prevLL = 0.0;
        for (int i = 1; i <= iterations; i++) {
            LOG.debug("Iteration " + i);
            double currLL = nextIteration(correctionConstant); // Core

            if (i > 1) {
                if (prevLL > currLL) {
                    LOG.error("Model Diverging: log likelihood decreased");
                    break;
                }
                if (currLL - prevLL < LLThreshold) {
                    break;
                }
            }
            prevLL = currLL;
        }

        // kill a bunch of these big objects now that we don't need them
        observedExpects = null;
        modelExpects = null;
        numTimesEventsSeen = null;
        trainingDataFeatNameIndices = null;
    }