PDF for Poisson prior

        Parameters
        ----------
        mu : float
            Latent variable for which the prior is being formed over

        Returns
        ----------
        - p(mu)
168404
def pdf(self, mu):
        
        if self.transform is not None:
            mu = self.transform(mu)                
        return ss.poisson.pmf(mu, self.lmd0)
900282
def zip_ll_row(params, data_row):
    
    l = params[0]
    pi = params[1]
    d0 = (data_row==0)
    likelihood = d0*pi + (1-pi)*poisson.pmf(data_row, l)
    return -np.log(likelihood+eps).sum()
965549
def toy_poisson_rbf_1d_laplace(optimize=True, plot=True):
    """"""
    optimizer='scg'
    x_len = 100
    X = np.linspace(0, 10, x_len)[:, None]
    f_true = np.random.multivariate_normal(np.zeros(x_len), GPy.kern.RBF(1).K(X))
    Y = np.array([np.random.poisson(np.exp(f)) for f in f_true])[:,None]

    kern = GPy.kern.RBF(1)
    poisson_lik = GPy.likelihoods.Poisson()
    laplace_inf = GPy.inference.latent_function_inference.Laplace()

    # create simple GP Model
    m = GPy.core.GP(X, Y, kernel=kern, likelihood=poisson_lik, inference_method=laplace_inf)

    if optimize:
        m.optimize(optimizer)
    if plot:
        m.plot()
        # plot the real underlying rate function
        pb.plot(X, np.exp(f_true), '--k', linewidth=2)

    return m
648559
def gaussian_prior_model_for_arguments(self, arguments):
        
        return CollectionPriorModel(
            {
                key: value.gaussian_prior_model_for_arguments(arguments)
                if isinstance(value, AbstractPriorModel)
                else value
                for key, value in self.__dict__.items() if key not in ('component_number', 'item_number', 'id')
            }
        )
1281840
def pdf(self, mu):
        
        if self.transform is not None:
            mu = self.transform(mu)    
        return self.pdf_internal(mu, df=self.df0, loc=self.loc0, scale=self.scale0, gamma=self.gamma0)
715307
def guass(self, mu: float, sigma: float) -> float:
        
        return float(
            lib.TCOD_random_get_gaussian_double(self.random_c, mu, sigma)
        )
636119
def variable_relative(self, r: float) -> mm.ModelMapper:
        
        return self.previous_variable.mapper_from_gaussian_tuples(self.gaussian_tuples, r=r)
1167254
def vp(self):
        
        return np.sqrt(self.mu * (2 / (self.rp) - 1 / self.kep.a))
1067128
def probability_density(self, X):
        
        self.check_fit()

        # make cov positive semi-definite
        covariance = self.covariance * np.identity(self.covariance.shape[0])
        return stats.multivariate_normal.pdf(X, cov=covariance)
403101
def pdf(self, mu):
        
        if self.transform is not None:
            mu = self.transform(mu)    
        if mu < self.lower and self.lower is not None:
            return 0.0
        elif mu > self.upper and self.upper is not None:
            return 0.0       
        else:
            return (1/float(self.sigma0))*np.exp(-(0.5*(mu-self.mu0)**2)/float(self.sigma0**2))
154391
def pdf(self, phi):
        

        """
        if self._n is None:
            pdf = np.empty_like(phi)
            pdf.fill(.5/np.pi)
            return pdf

        phi = np.asarray(phi)

        pdf = self._pdf(phi)
        pdf /= 2.*np.pi

        return pdf
739783
def draw_latent_variables(self, nsims=5000):
        
        if self.latent_variables.estimation_method is None:
            raise Exception("No latent variables estimated!")
        elif self.latent_variables.estimation_method == 'BBVI':
            return np.array([i.q.draw_variable_local(size=nsims) for i in self.latent_variables.z_list])
        elif self.latent_variables.estimation_method == "M-H":
            chain = np.array([self.latent_variables.z_list[i].sample for i in range(len(self.latent_variables.z_list))])
            return chain[:,np.random.choice(chain.shape[1], nsims)]
        else:
            raise Exception("No latent variables estimated through Bayesian inference")
496377
def rnegative_binomial(mu, alpha, size=None):
    
    # Using gamma-poisson mixture rather than numpy directly
    # because numpy apparently rounds
    mu = np.asarray(mu, dtype=float)
    pois_mu = np.random.gamma(alpha, mu / alpha, size)
    return np.random.poisson(pois_mu, size)
1056822
def Poisson(mu: vertex_constructor_param_types, label: Optional[str]=None) -> Vertex:
    
    return Integer(context.jvm_view().PoissonVertex, label, cast_to_double_vertex(mu))
111337
def add_poisson(image, exp_time):
    
    """
    adds a poison (or Gaussian) distributed noise with mean given by surface brightness
    """
    if isinstance(exp_time, int) or isinstance(exp_time, float):
        if exp_time <= 0:
            exp_time = 1
    else:
        mean_exp_time = np.mean(exp_time)
        exp_time[exp_time < mean_exp_time/10] = mean_exp_time/10
    sigma = np.sqrt(np.abs(image)/exp_time) # Gaussian approximation for Poisson distribution, normalized to exposure time
    nx, ny = np.shape(image)
    poisson = np.random.randn(nx, ny) * sigma
    return poisson
1026286
def rmv_normal_cov(mu, C, size=1):
    
    mu_size = np.shape(mu)
    if size == 1:
        return np.random.multivariate_normal(mu, C, size).reshape(mu_size)
    else:
        return np.random.multivariate_normal(
            mu, C, size).reshape((size,) + mu_size)
518376
def inference_likelihood(self, kern, X, Z, likelihood, Y):
        

        num_data, output_dim = Y.shape
        input_dim = Z.shape[0]
        if self.mpi_comm != None:
            from mpi4py import MPI
            num_data_all = np.array(num_data,dtype=np.int32)
            self.mpi_comm.Allreduce([np.int32(num_data), MPI.INT], [num_data_all, MPI.INT])
            num_data = num_data_all

        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
        if het_noise:
            self.batchsize = 1

        psi0_full, psi1Y_full, psi2_full, YRY_full = self.gatherPsiStat(kern, X, Z, Y, beta, uncertain_inputs)

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

        Kmm = kern.K(Z).copy()
        diag.add(Kmm, self.const_jitter)
        if not np.isfinite(Kmm).all():
            print(Kmm)
        Lm = jitchol(Kmm)
        LmInv = dtrtri(Lm)

        LmInvPsi2LmInvT = LmInv.dot(psi2_full.dot(LmInv.T))
        Lambda = np.eye(Kmm.shape[0])+LmInvPsi2LmInvT
        LL = jitchol(Lambda)
        LLInv = dtrtri(LL)
        logdet_L = 2.*np.sum(np.log(np.diag(LL)))
        LmLLInv = LLInv.dot(LmInv)
        
        b  = psi1Y_full.dot(LmLLInv.T)
        bbt = np.square(b).sum()
        v = b.dot(LmLLInv).T
        LLinvPsi1TYYTPsi1LLinvT = tdot(b.T)
        
        tmp = -LLInv.T.dot(LLinvPsi1TYYTPsi1LLinvT+output_dim*np.eye(input_dim)).dot(LLInv)
        dL_dpsi2R = LmInv.T.dot(tmp+output_dim*np.eye(input_dim)).dot(LmInv)/2.
        
        # Cache intermediate results
        self.midRes['dL_dpsi2R'] = dL_dpsi2R
        self.midRes['v'] = v

        #======================================================================
        # Compute log-likelihood
        #======================================================================
        if het_noise:
            logL_R = -np.log(beta).sum()
        else:
            logL_R = -num_data*np.log(beta)
        logL = -(output_dim*(num_data*log_2_pi+logL_R+psi0_full-np.trace(LmInvPsi2LmInvT))+YRY_full-bbt)/2.-output_dim*logdet_L/2.

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

        dL_dKmm =  dL_dpsi2R - output_dim*LmInv.T.dot(LmInvPsi2LmInvT).dot(LmInv)/2.

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

        if not self.Y_speedup or het_noise:
            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=v, K=Kmm, mean=None, cov=None, K_chol=Lm)
        else:
            post = None

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

        if not het_noise:
            dL_dthetaL = (YRY_full*beta + beta*output_dim*psi0_full - num_data*output_dim*beta)/2. - beta*(dL_dpsi2R*psi2_full).sum() - beta*(v.T*psi1Y_full).sum()
            self.midRes['dL_dthetaL'] = dL_dthetaL

        return logL, dL_dKmm, post