我们从Python开源项目中，提取了以下50个代码示例，用于说明如何使用scipy.special.gammaln()。
def lda_e_step(doc_word_ids, doc_word_counts, alpha, beta, max_iter=100): gamma = np.ones(len(alpha)) expElogtheta = np.exp(dirichlet_expectation(gamma)) betad = beta[:, doc_word_ids] phinorm = np.dot(expElogtheta, betad) + 1e-100 counts = np.array(doc_word_counts) for _ in xrange(max_iter): lastgamma = gamma gamma = alpha + expElogtheta * np.dot(counts / phinorm, betad.T) Elogtheta = dirichlet_expectation(gamma) expElogtheta = np.exp(Elogtheta) phinorm = np.dot(expElogtheta, betad) + 1e-100 meanchange = np.mean(abs(gamma - lastgamma)) if (meanchange < meanchangethresh): break likelihood = np.sum(counts * np.log(phinorm)) likelihood += np.sum((alpha - gamma) * Elogtheta) likelihood += np.sum(sp.gammaln(gamma) - sp.gammaln(alpha)) likelihood += sp.gammaln(np.sum(alpha)) - sp.gammaln(np.sum(gamma)) return (likelihood, gamma)
def test_value(self): with self.test_session(use_gpu=True): def _test_value(given, temperature, logits): given = np.array(given, np.float32) logits = np.array(logits, np.float32) n = logits.shape[-1] t = temperature target_log_p = special.gammaln(n) + (n - 1) * np.log(t) + \ (logits - t * given).sum(axis=-1) - \ n * np.log(np.exp(logits - t * given).sum(axis=-1)) con = ExpConcrete(temperature, logits=logits) log_p = con.log_prob(given) self.assertAllClose(log_p.eval(), target_log_p) p = con.prob(given) self.assertAllClose(p.eval(), np.exp(target_log_p)) _test_value([np.log(0.25), np.log(0.25), np.log(0.5)], 0.1, [1., 1., 1.2]) _test_value([[np.log(0.25), np.log(0.25), np.log(0.5)], [np.log(0.1), np.log(0.5), np.log(0.4)]], 0.5, [[1., 1., 1.], [.5, .5, .4]])
def test_value(self): with self.test_session(use_gpu=True): def _test_value(given, temperature, logits): given = np.array(given, np.float32) logits = np.array(logits, np.float32) n = logits.shape[-1] t = temperature target_log_p = special.gammaln(n) + (n - 1) * np.log(t) + \ (logits - (t + 1) * np.log(given)).sum(axis=-1) - \ n * np.log(np.exp(logits - t * np.log(given)).sum(axis=-1)) con = Concrete(temperature, logits=logits) log_p = con.log_prob(given) self.assertAllClose(log_p.eval(), target_log_p) p = con.prob(given) self.assertAllClose(p.eval(), np.exp(target_log_p)) _test_value([0.25, 0.25, 0.5], 0.1, [1., 1., 1.2]) _test_value([[0.25, 0.25, 0.5], [0.1, 0.5, 0.4]], 0.5, [[1., 1., 1.], [.5, .5, .4]])
def neg_loglikelihood(y, mean, scale, shape, skewness): """ Negative loglikelihood function Parameters ---------- y : np.ndarray univariate time series mean : np.ndarray array of location parameters for the Poisson distribution scale : float scale parameter for the Poisson distribution shape : float tail thickness parameter for the Poisson distribution skewness : float skewness parameter for the Poisson distribution Returns ---------- - Negative loglikelihood of the Poisson family """ return -np.sum(-mean + np.log(mean)*y - sp.gammaln(y + 1))
def sample_gamma_posterior (data, shape0, scale0): S = numpy.sum(data) Slog = numpy.sum(map(numpy.log, data)) N = len(data) logP = 1.0+Slog q = 1.0 + S r = 1.0 + N s = 1.0 + N # lnf_shape = lambda al: (al-1)*logP + special.gammaln(s*al+1) - (1-al*s)*numpy.log(q) -r*special.gammaln(al) lnf_shape = lambda al: (al-1)*logP - q/scale0 - r*special.gammaln(al) - (al*s)*numpy.log(scale0) shape_new = sampling.slice (lnf_shape, shape0, thin=5, N=1, lower=0, upper=numpy.inf)[0] # print "Sfoo", shape0, lnf_shape(shape0) # print "....", shape_new, lnf_shape(shape_new) # lnf_rate = lambda beta: -beta*q + s*shape_new*numpy.log(beta) # rate_new = sampling.slice (lnf_rate, 1.0/scale0, thin=5, N=1, lower=0, upper=numpy.inf)[0] lnf_scale = lambda beta: (shape_new-1)*logP - q/beta - r*special.gammaln(shape_new) - (shape_new*s)*numpy.log(beta) scale_new = sampling.slice (lnf_scale, scale0, thin=5, N=1, lower=0, upper=numpy.inf)[0] # print "Rfoo", 1.0/scale0, lnf_scale(scale0) # print "...", scale_new, lnf_scale(scale_new) return [shape_new, scale_new]
def gamma_transition_kernel (data, shape0, scale0, shape1, scale1): S = numpy.sum(data) Slog = numpy.sum(map(numpy.log, data)) N = len(data) logP = 1.0+Slog q = 1.0 + S r = 1.0 + N s = 1.0 + N lnf_shape = lambda al: (al-1)*logP - q/scale0 - r*special.gammaln(al) - (al*s)*numpy.log(scale0) p_shape = lnf_shape(shape1) lnf_scale = lambda beta: (shape1-1)*logP - q/beta - r*special.gammaln(shape1) - (shape1*s)*numpy.log(beta) p_scale = lnf_scale (scale1) return p_shape + p_scale # proposal
def log_marg_for_copy(self, copy_components): """Return log marginal of data and component assignments: p(X, z)""" # Log probability of component assignment P(z|alpha) # Equation (10) in Wood and Black, 2008 # Use \Gamma(n) = (n - 1)! facts_ = gammaln(copy_components.counts[:copy_components.K]) facts_[copy_components.counts[:copy_components.K] == 0] = 0 # definition of log(0!) log_prob_z = ( (copy_components.K - 1)*math.log(self.alpha) + gammaln(self.alpha) - gammaln(np.sum(copy_components.counts[:copy_components.K]) + self.alpha) + np.sum(facts_) ) log_prob_X_given_z = copy_components.log_marg() return log_prob_z + log_prob_X_given_z
def log_marg(self): """Return log marginal of data and component assignments: p(X, z)""" # Log probability of component assignment P(z|alpha) # Equation (10) in Wood and Black, 2008 # Use \Gamma(n) = (n - 1)! facts_ = gammaln(self.components.counts[:self.components.K]) facts_[self.components.counts[:self.components.K] == 0] = 0 # definition of log(0!) log_prob_z = ( (self.components.K - 1)*math.log(self.alpha) + gammaln(self.alpha) - gammaln(np.sum(self.components.counts[:self.components.K]) + self.alpha) + np.sum(facts_) ) log_prob_X_given_z = self.components.log_marg() return log_prob_z + log_prob_X_given_z
def log_marg(self): """Return log marginal of data and component assignments: p(X, z)""" # Log probability of component assignment, (24.24) in Murphy, p. 842 log_prob_z = ( gammaln(self.kalpha) - gammaln(self.kalpha + np.sum(self.components.counts)) + np.sum( gammaln( self.components.counts + float(self.kalpha)/self.components.K_max ) - gammaln(self.kalpha/self.components.K_max) ) ) # Log probability of data in each component log_prob_X_given_z = self.components.log_marg() return log_prob_z + log_prob_X_given_z
def _ll(self, x): result = 0.0 # P(w|z) result += self.K * special.gammaln(self.beta * self.K) result += -np.sum(special.gammaln(np.sum(self.n_z_w, axis=1))) result += np.sum(special.gammaln(self.n_z_w)) result += -self.V * special.gammaln(self.beta) # P(z|Lambda) alpha = self.get_alpha(x) result += np.sum(special.gammaln(np.sum(alpha, axis=1))) result += -np.sum(special.gammaln( np.sum(self.n_m_z+alpha, axis=1))) result += np.sum(special.gammaln(self.n_m_z+alpha)) result += -np.sum(special.gammaln(alpha)) # P(Lambda) result += -self.K / 2.0 * np.log(2.0 * np.pi * (self.sigma ** 2)) result += -1.0 / (2.0 * (self.sigma ** 2)) * np.sum(x ** 2) result = -result return result
def visit_fn(self, temperature): factor1 = np.exp(np.log(temperature) / (self.qv - 1.0)) factor2 = np.exp((4.0 - self.qv) * np.log(self.qv - 1.0)) factor3 = np.exp((2.0 - self.qv) * np.log(2.0) / (self.qv - 1.0)) factor4 = np.sqrt(np.pi) * factor1 * factor2 / (factor3 * ( 3.0 - self.qv)) factor5 = 1.0 / (self.qv - 1.0) - 0.5 d1 = 2.0 - factor5 factor6 = np.pi * (1.0 - factor5) / np.sin( np.pi * (1.0 - factor5)) / np.exp(gammaln(d1)) sigmax = np.exp(-(self.qv - 1.0) * np.log( factor6 / factor4) / (3.0 - self.qv)) x = sigmax * self.gaussian_fn(1) y = self.gaussian_fn(0) den = np.exp( (self.qv - 1.0) * np.log((np.fabs(y))) / (3.0 - self.qv)) return x / den
def prob_jk(self, j, k): # -1 because table of current sample topic jk, is not conditioned on njdotk = self.count_k_by_j[j, k] if njdotk == 1: return np.ones(1) possible_ms = np.arange(1, njdotk) # +1-1 log_alpha_beta_k = self.get_log_alpha_beta(k) alpha_beta_k = np.exp(log_alpha_beta_k) normalizer = gammaln(alpha_beta_k) - gammaln(alpha_beta_k + njdotk) log_stir = self.stirling_mat(njdotk, possible_ms) params = normalizer + log_stir + possible_ms*log_alpha_beta_k return lognormalize(params)
def init_with_data(cls, lda_data, nhaps=3, haplotypes=None, use_ll_offset=False): V = len(lda_data.alphabets) status = cls(lda_data) zero = 1e-30 # value for zero count status.nhaps = nhaps status.ll_offset = gammaln(nhaps + 1) if use_ll_offset else 0. status.alpha = np.array([1. for _ in xrange(nhaps)]) # Uniform distribution #status.alpha = np.array([1. / nhaps for _ in xrange(nhaps)]) if haplotypes: status.log_beta = np.log(np.array([[[1 if a == haplotypes[n][m] else zero for n in xrange(nhaps)] for a in lda_data.alphabets] for m in xrange(lda_data.nsites)])) else: status.log_beta = np.array([[[np.random.random() for _ in xrange(nhaps)] for _ in lda_data.alphabets] for _ in xrange(lda_data.nsites)]) status.log_beta = log_normalize(status.log_beta, axis=1) # adhoc initialization status.log_psi = status.log_beta[np.newaxis, :, :, :] status.gamma = status.alpha[np.newaxis, :] + np.exp(lda_data.log_d[:, :, :, np.newaxis] + status.log_psi).sum(axis=1).sum(axis=1) status.ll = status.calc_lower_ll() #(status.gamma, status.log_psi) = _calc_gamma_psi(log_d, status.alpha, status.log_beta, gamma, log_psi) return status
def log_p_kmer_table(kmer_counts, p_kmers): """log_p_kmer_table calculates the probabilty of an observed kmer table given as a multinomial given the multinomial probabilities. The kmer table is preferably collapsed, i.e. if a given sequence is present, its reverse complement is not. Args: kmer_counts: the kmer table p_kmers: the multinomial probabilities Returns: the total log probability as a float """ total = 0 log_p = 0 total_p = sum(p_kmers.values()) for kmer in p_kmers.keys(): p_kmers[kmer] /= total_p for kmer, count in kmer_counts.iteritems(): log_p += count*np.log(p_kmers[kmer]) - special.gammaln(count+1) total += count log_p += special.gammaln(total+1) return log_p
def log_p_multinomial(X, p): """log_p_multinomial returns the probability of drawing a vector of counts from a multinomial parameterized by p. Args: X: np.array of counts representing a multinomial draw. p: np.array of probabilities, corresponding to X Returns: the log probability of the multinomial draw. """ if sum(X) == 0: return 0.0 # check that input is valid assert len(X) == len(p) eps = 0.0001 assert abs(sum(p)- 1.0) < eps # calculate log prob. log_n_choices = special.gammaln(sum(X)+1) - sum([special.gammaln(x_i+1) for x_i in X]) log_p_items = sum(x_i*np.log(p_i) for (x_i, p_i) in zip(X, p)) return log_n_choices + log_p_items
def log_prob_given_alpha(self): """log_prob_given_alpha returns the probability of the mixing proportions, gamma, under the CRP prior Returns: the log probability as a float. """ K = len(self.component_counts.keys()) # number of seen components N = sum(self.component_counts.values()) log_numerator = K*np.log(self.alpha_crp) + sum( special.gammaln(N_k) for N_k in self.component_counts.values() ) log_denom = (special.gammaln(N + self.alpha_crp) - special.gammaln(self.alpha_crp)) return log_numerator - log_denom
def makePlot_pdf_Phi( nu=0, tau=0, phi_grid=None, ngrid=1000, min_phi=-100, max_phi=100): label = 'nu=%7.2f' % (nu) cPrior = - gammaln(nu) + gammaln(nu-tau) + gammaln(tau) if phi_grid is None: phi_grid = np.linspace(min_phi, max_phi, ngrid) logpdf_grid = tau * phi_grid - nu * c_Phi(phi_grid) - cPrior pdf_grid = np.exp(logpdf_grid) IntegralVal = np.trapz(pdf_grid, phi_grid) mu_grid = phi2mu(phi_grid) ExpectedPhiVal = np.trapz(pdf_grid * phi_grid, phi_grid) ExpectedMuVal = np.trapz(pdf_grid * mu_grid, phi_grid) print '%s Integral=%.4f E[phi]=%6.3f E[mu]=%.4f' % ( label, IntegralVal, ExpectedPhiVal, ExpectedMuVal) pylab.plot(phi_grid, pdf_grid, '-', label=label) pylab.xlabel('phi (log odds ratio)') pylab.ylabel('density p(phi)')
def makePlot_pdf_Mu( nu=0, tau=0, phi_grid=None, ngrid=1000, min_phi=-100, max_phi=100): label = 'nu=%7.2f' % (nu,) cPrior = - gammaln(nu) + gammaln(nu-tau) + gammaln(tau) mu_grid = np.linspace(1e-15, 1-1e-15, ngrid) phi_grid = mu2phi(mu_grid) logpdf_grid = tau * phi_grid - nu * c_Phi(phi_grid) - cPrior logJacobian_grid = logJacobian_mu(mu_grid) pdf_grid = np.exp(logpdf_grid + logJacobian_grid) IntegralVal = np.trapz(pdf_grid, mu_grid) ExpectedMuVal = np.trapz(pdf_grid * mu_grid, mu_grid) ExpectedPhiVal = np.trapz(pdf_grid * phi_grid, mu_grid) print '%s Integral=%.4f E[phi]=%6.3f E[mu]=%.4f' % ( label, IntegralVal, ExpectedPhiVal, ExpectedMuVal) pylab.plot(mu_grid, pdf_grid, '-', label=label) pylab.xlabel('mu') pylab.ylabel('density p(mu)')
def L_hdp(beta, omega, Tvec, alpha): ''' Compute top-tier of hdp bound. ''' K = omega.size assert K == beta.size assert K == Tvec.size rho = OROB.beta2rho(beta, K) eta1 = rho * omega eta0 = (1-rho) * omega digamma_omega = digamma(omega) digamma_eta1 = digamma(eta1) digamma_eta0 = digamma(eta0) ElogU = digamma_eta1 - digamma_omega Elog1mU = digamma_eta0 - digamma_omega Lscore = \ np.sum(gammaln(eta1) + gammaln(eta0) - gammaln(omega)) \ + np.inner(nDoc + 1 - eta1, ElogU) \ + np.inner(nDoc * OROB.kvec(K) + gamma - eta0, Elog1mU) \ + alpha * np.inner(beta, Tvec) return Lscore
def c_Func(nu, logdetB, M, logdetV): ''' Evaluate cumulant function c at given parameters. c is the cumulant of the MatrixNormal-Wishart, using common params. Returns -------- c : scalar real value of cumulant function at provided args ''' if logdetB.ndim >= 2: logdetB = np.log(np.linalg.det(logdetB)) if logdetV.ndim >= 2: logdetV = np.log(np.linalg.det(logdetV)) D, E = M.shape dvec = np.arange(1, D + 1, dtype=np.float) return - 0.25 * D * (D - 1) * LOGPI \ - 0.5 * D * LOGTWO * nu \ - np.sum(gammaln(0.5 * (nu + 1 - dvec))) \ + 0.5 * nu * logdetB \ - 0.5 * D * E * LOGTWOPI \ + 0.5 * E * logdetV
def c_Diff(nu1, logdetB1, D, nu2, logdetB2): ''' Evaluate difference of cumulant functions c(params1) - c(params2) May be more numerically stable than directly using c_Func to find the difference. Returns ------- diff : scalar real value of the difference in cumulant functions ''' if logdetB1.ndim >= 2: assert D == logdetB1.shape[-1] logdetB1 = np.log(np.linalg.det(logdetB1)) logdetB2 = np.log(np.linalg.det(logdetB2)) dvec = np.arange(1, D + 1, dtype=np.float) return - 0.5 * D * LOGTWO * (nu1 - nu2) \ - np.sum(gammaln(0.5 * (nu1 + 1 - dvec))) \ + np.sum(gammaln(0.5 * (nu2 + 1 - dvec))) \ + 0.5 * (nu1 * logdetB1 - nu2 * logdetB2)
def _calcPredProbVec_Fast(self, SS, x): p = np.zeros(SS.K) nu, beta, m, kappa = self.calcPostParams(SS) kbeta = beta kbeta *= ((kappa + 1) / kappa)[:, np.newaxis] base = np.square(x - m) base /= kbeta base += 1 # logp : 2D array, size K x D logp = (-0.5 * (nu + 1))[:, np.newaxis] * np.log(base) logp += (gammaln(0.5 * (nu + 1)) - gammaln(0.5 * nu))[:, np.newaxis] logp -= 0.5 * np.log(kbeta) # p : 1D array, size K p = np.sum(logp, axis=1) p -= np.max(p) np.exp(p, out=p) return p
def c_Diff(nu1, beta1, m1, kappa1, nu2, beta2, m2, kappa2): ''' Evaluate difference of cumulant functions c(params1) - c(params2) May be more numerically stable than directly using c_Func to find the difference. Returns ------- diff : scalar real value of the difference in cumulant functions ''' cDiff = - 0.5 * LOGTWO * (nu1 - nu2) \ - gammaln(0.5 * nu1) + gammaln(0.5 * nu2) \ + 0.5 * (np.log(kappa1) - np.log(kappa2)) \ + 0.5 * (nu1 * np.log(beta1) - nu2 * np.log(beta2)) return np.sum(cDiff)
def _calcPredProbVec_Fast(self, SS, x): nu, B, m, kappa = self.calcPostParams(SS) kB = B kB *= ((kappa + 1) / kappa)[:, np.newaxis, np.newaxis] logp = np.zeros(SS.K) p = logp # Rename so its not confusing what we're returning for k in xrange(SS.K): cholKB = scipy.linalg.cholesky(kB[k], lower=1) logdetKB = 2 * np.sum(np.log(np.diag(cholKB))) mVec = np.linalg.solve(cholKB, x - m[k]) mDist_k = np.inner(mVec, mVec) logp[k] = -0.5 * logdetKB - 0.5 * \ (nu[k] + 1) * np.log(1.0 + mDist_k) logp += gammaln(0.5 * (nu + 1)) - gammaln(0.5 * (nu + 1 - self.D)) logp -= np.max(logp) np.exp(logp, out=p) return p
def c_Func(nu, logdetB, m, kappa): ''' Evaluate cumulant function at given params. Returns -------- c : scalar real value of cumulant function at provided args ''' if logdetB.ndim >= 2: logdetB = np.log(np.linalg.det(logdetB)) D = m.size dvec = np.arange(1, D + 1, dtype=np.float) return - 0.5 * D * LOGTWOPI \ - 0.25 * D * (D - 1) * LOGPI \ - 0.5 * D * LOGTWO * nu \ - np.sum(gammaln(0.5 * (nu + 1 - dvec))) \ + 0.5 * D * np.log(kappa) \ + 0.5 * nu * logdetB
def c_Diff(nu1, logdetB1, m1, kappa1, nu2, logdetB2, m2, kappa2): ''' Evaluate difference of cumulant functions c(params1) - c(params2) May be more numerically stable than directly using c_Func to find the difference. Returns ------- diff : scalar real value of the difference in cumulant functions ''' if logdetB1.ndim >= 2: logdetB1 = np.log(np.linalg.det(logdetB1)) if logdetB2.ndim >= 2: logdetB2 = np.log(np.linalg.det(logdetB2)) D = m1.size dvec = np.arange(1, D + 1, dtype=np.float) return - 0.5 * D * LOGTWO * (nu1 - nu2) \ - np.sum(gammaln(0.5 * (nu1 + 1 - dvec))) \ + np.sum(gammaln(0.5 * (nu2 + 1 - dvec))) \ + 0.5 * D * (np.log(kappa1) - np.log(kappa2)) \ + 0.5 * (nu1 * logdetB1 - nu2 * logdetB2)
def c_Beta(g1, g0): ''' Calculate cumulant function of the Beta distribution Input can be vectors, in which case we compute sum over several cumulant functions of the independent distributions: \prod_k Beta(g1[k], g0[k]) Args ---- g1 : 1D array, size K first parameter of a Beta distribution g0 : 1D array, size K second parameter of a Beta distribution Returns ------- c : scalar sum of the cumulants defined by provided parameters ''' return np.sum(gammaln(g1 + g0) - gammaln(g1) - gammaln(g0))
def calcELBO_SingleDocFromSparseResp( spResp_d, ElogLik_d, wc_d, alphaEbeta): ''' Calculate single document contribution to the ELBO objective. This isolates all ELBO terms that depend on local parameters of this doc. Returns ------- L : scalar float value of ELBO objective, up to additive constant. This constant is independent of any local parameter attached to doc d. ''' DocTopicCount_d = wc_d * spResp_d # Sparse multiply theta_d = DocTopicCount_d + alphaEbeta R_d = spResp_d.toarray() wR_d = wc_d[:, np.newaxis] * R_d L_alloc = np.sum(gammaln(theta_d)) L_data = np.sum(wR_d * ElogLik_d) RlogR = np.sum(wR_d * np.log(R_d + 1e-100)) return L_alloc + L_data - RlogR
def log_prob_z(self): """ Return the log marginal probability of component assignment P(z). See (24.24) in Murphy, p. 842. """ log_prob_z = ( gammaln(self.alpha) - gammaln(self.alpha + np.sum(self.components.counts)) + np.sum( gammaln( self.components.counts + float(self.alpha)/self.components.K_max ) - gammaln(self.alpha/self.components.K_max) ) ) return log_prob_z
def log_marg(self): """Return log marginal of data and component assignments: p(X, z)""" log_prob_z = self.log_prob_z() log_prob_X_given_z = self.log_prob_X_given_z() # # Log probability of component assignment, (24.24) in Murphy, p. 842 # log_prob_z = ( # gammaln(self.alpha) # - gammaln(self.alpha + np.sum(self.components.counts)) # + np.sum( # gammaln( # self.components.counts # + float(self.alpha)/self.components.K_max # ) # - gammaln(self.alpha/self.components.K_max) # ) # ) # # Log probability of data in each component # log_prob_X_given_z = self.components.log_marg() return log_prob_z + log_prob_X_given_z # @profile
def _vlb(self): vlb = 0. vlb += sum(l.get_vlb() for l in self.labels_list) vlb += self.weights.get_vlb() vlb += sum(c.get_vlb() for c in self.components) for l in self.labels_list: vlb += np.sum([r.dot(c.expected_log_likelihood(l.data)) for c,r in zip(self.components, l.r.T)]) # add in symmetry factor (if we're actually symmetric) if len(set(type(c) for c in self.components)) == 1: vlb += special.gammaln(len(self.components)+1) return vlb ### SVI
def get_vlb(self): # E[ln p(mu) / q(mu)] part h_0, J_0, h_mf, J_mf = self.h_0, self.J_0, self.h_mf, self.J_mf Emu, Emu2 = self._E_mu p_mu_avgengy = -1./2*J_0*Emu2 + h_0*Emu \ - 1./2*(h_0**2/J_0) + 1./2*np.log(J_0) - 1./2*np.log(2*np.pi) q_mu_entropy = 1./2*np.log(2*np.pi*np.e/J_mf) # E[ln p(sigmasq) / q(sigmasq)] part alpha_0, beta_0, alpha_mf, beta_mf = \ self.alpha_0, self.beta_0, self.alpha_mf, self.beta_mf (Esigmasqinv, Elnsigmasq) = self._E_sigmasq p_sigmasq_avgengy = (-alpha_0-1)*Elnsigmasq + (-beta_0)*Esigmasqinv \ - (special.gammaln(alpha_0) - alpha_0*np.log(beta_0)) q_sigmasq_entropy = alpha_mf + np.log(beta_mf) + special.gammaln(alpha_mf) \ - (1+alpha_mf)*special.digamma(alpha_mf) return p_mu_avgengy + q_mu_entropy + p_sigmasq_avgengy + q_sigmasq_entropy
def get_vlb(self): # E[ln p(mu) / q(mu)] part h_0, J_0, h_mf, J_mf = self.h_0, self.J_0, self.h_mf, self.J_mf Emu, Emu2 = self._E_mu p_mu_avgengy = -1./2*J_0*Emu2 + h_0*Emu \ - 1./2*(h_0**2/J_0) + 1./2*np.log(J_0) - 1./2*np.log(2*np.pi) q_mu_entropy = 1./2*np.log(2*np.pi*np.e/J_mf) # E[ln p(sigmasq) / q(sigmasq)] part alpha_0, beta_0, alpha_mf, beta_mf = \ self.alpha_0, self.beta_0, self.alpha_mf, self.beta_mf (Esigmasqinv, Elnsigmasq) = self._E_sigmasq p_sigmasq_avgengy = (-alpha_0-1)*Elnsigmasq + (-beta_0)*Esigmasqinv \ - (special.gammaln(alpha_0) - alpha_0*np.log(beta_0)) q_sigmasq_entropy = alpha_mf + np.log(beta_mf) + special.gammaln(alpha_mf) \ - (1+alpha_mf)*special.digamma(alpha_mf) return np.sum(p_mu_avgengy + q_mu_entropy + p_sigmasq_avgengy + q_sigmasq_entropy)
def log_likelihood(self,x,r=None,p=None): r = r if r is not None else self.r p = p if p is not None else self.p x = np.array(x,ndmin=1) if self.p > 0: xnn = x[x >= 0] raw = np.empty(x.shape) raw[x>=0] = special.gammaln(r + xnn) - special.gammaln(r) \ - special.gammaln(xnn+1) + r*np.log(1-p) + xnn*np.log(p) raw[x<0] = -np.inf return raw if isinstance(x,np.ndarray) else raw[0] else: raw = np.log(np.zeros(x.shape)) raw[x == 0] = 0. return raw if isinstance(x,np.ndarray) else raw[0]
def _get_statistics(self,data): # NOTE: since this isn't really in exponential family, this method needs # to look at hyperparameters. form posterior hyperparameters for the p # parameters here so we can integrate them out and get the r statistics n, tot = super(NegativeBinomialIntegerR,self)._get_statistics(data) if n > 0: alpha_n, betas_n = self.alpha_0 + tot, self.beta_0 + self.r_support*n data = flattendata(data) log_marg_likelihoods = \ special.betaln(alpha_n, betas_n) \ - special.betaln(self.alpha_0, self.beta_0) \ + (special.gammaln(data[:,na]+self.r_support) - special.gammaln(data[:,na]+1) \ - special.gammaln(self.r_support)).sum(0) else: log_marg_likelihoods = np.zeros_like(self.r_support) return n, tot, log_marg_likelihoods
def log_prob(self, x): def _compute(df, loc, scale_tril, x): k = scale_tril.shape[-1] ildj = np.sum(np.log(np.abs(np.diag(scale_tril))), axis=-1) logz = ildj + k * (0.5 * np.log(df) + 0.5 * np.log(np.pi) + special.gammaln(0.5 * df) - special.gammaln(0.5 * (df + 1.))) y = linalg.solve_triangular(scale_tril, np.matrix(x - loc).T, lower=True, overwrite_b=True) logs = -0.5 * (df + 1.) * np.sum(np.log1p(y**2. / df), axis=-2) return logs - logz if not self._df.shape: return _compute(self._df, self._loc, self._scale_tril, x) return np.concatenate([ [_compute(self._df[i], self._loc[i], self._scale_tril[i], x[:, i, :])] for i in range(len(self._df))]).T
def logprob_dc(counts, prior, axis=None): """Non-normalized log probability of a Dirichlet-Categorical distribution. See https://en.wikipedia.org/wiki/Dirichlet-multinomial_distribution """ # Note that this excludes the factorial(counts) term, since we explicitly # track permutations in assignments. return gammaln(np.add(counts, prior, dtype=np.float32)).sum(axis)
def neg_bin_log_pmf(k, mu, phi): """ Log(PMF) of negative binomial distribution with mean mu and dispersion phi, conveniently parameterized. Args: k (int) - NB random variable mu (float) - mean phi (float) - dispersion Returns: The log of the pmf at k. """ r = 1.0/phi return gammaln(r+k) - (gammaln(r) + gammaln(k+1)) + \ k * np.log(mu/(r+mu)) + \ r * np.log(r/(r+mu))
def trunc_logLik(data, mu, alpha): log_1_plus_a_mu = np.log(1 + alpha*mu) log_1_minus_prob_zero = np.log(1.0 - np.exp(-np.log(1.0+alpha*mu)/alpha)) alpha_inv = 1.0/alpha lim = int(np.max(data.keys())) holding_val=0.0 log_L=0.0 for i in range(1, lim+1): holding_val += np.log(1+alpha*(i-1)) log_L += data[i]* (holding_val - special.gammaln(i) + i*np.log(mu)-(i+alpha_inv)*log_1_plus_a_mu - log_1_minus_prob_zero) return log_L
def log_multinomial_coefficient(list_of_primitive_rules): """ Count copies of each primitive and calculate a log multinomial coefficient. ~ """ list_of_tuples = [p.as_tuple() for p in list_of_primitive_rules] counts = dict.fromkeys(set(list_of_tuples),0) for t in list_of_tuples: counts[t] += 1 N = len(list_of_primitive_rules) v = np.array(counts.values()) return gammaln(N+1)-np.sum(gammaln(v+1))
def log_multinomial_coefficient(list_of_tuples): counts = dict.fromkeys(set(list_of_tuples),0) for t in list_of_tuples: counts[t] += 1 N = len(list_of_tuples) v = np.array(counts.values()) return gammaln(N+1)-np.sum(gammaln(v+1)) #### # Parameters (this should be made user-configurable later)
def test_serialize(self): from scipy.special import gammaln x = range(1, 5) expected = list(map(gammaln, x)) observed = self.sc.parallelize(x).map(gammaln).collect() self.assertEqual(expected, observed)
def estimate_with_fixed_beta(self, data, beta): mu = median(data) v = mean((data - mu) ** 2) alpha = sqrt(v * exp(gammaln(1. / beta) - gammaln(3. / beta))) return mu, alpha
def log_prob(self, x): mu = self.mu alpha = self.alpha beta = self.beta return log(beta / (2.0 * alpha)) - gammaln(1. / beta) - power(fabs(x - mu) / alpha, beta)
def log_prob(self, x): a, b = self['alpha'], self['beta'] return a * log(b) - gammaln(clip(a, 1e-308, 1e308)) + \ (a - 1) * log(clip(x, 1e-308, 1e308)) - b * x
