On Mon, Jul 2, 2012 at 10:53 PM, Fernando Perez <fperez....@gmail.com> wrote: > On Mon, Jul 2, 2012 at 7:49 PM, Fernando Perez <fperez....@gmail.com> wrote: >> It appears you're right!
nice idea: https://github.com/pymc-devs/pymc/blob/master/pymc/distributions.py#L1670 >> >> http://pymc-devs.github.com/pymc/distributions.html?highlight=hypergeometric#pymc.distributions.multivariate_hypergeometric_like > > Furthermore, the code actually calls a sampler implemented in Fortran: > > http://pymc-devs.github.com/pymc/_modules/pymc/distributions.html#multivariate_hypergeometric_like > > which calls this: > > https://github.com/pymc-devs/pymc/blob/master/pymc/flib.f#L4379 > > Thanks again to both of you for the help. It's always productive to > ask around here ;) the loglikelihood function is just some calls to scipy.special.gammaln I changed my version to get better numerical precision #taken from scipy.misc but returns log of comb def log_comb(n, k): from scipy import special from scipy.special import gammaln k = np.asarray(k) n = np.asarray(n) cond = (k <= n) & (n >= 0) & (k >= 0) sv = special.errprint(0) vals = gammaln(n + 1) - gammaln(n - k + 1) - gammaln(k + 1) sv = special.errprint(sv) return np.where(cond, vals, -np.inf) def log_pmf(x, n_balls): x = np.asarray(x) n_balls = np.asarray(n_balls) ret = np.sum(log_comb(n_balls, x)) - log_comb(n_balls.sum(), x.sum()) return ret def pmf(x, n_balls): x = np.asarray(x) n_balls = np.asarray(n_balls) ret = np.sum(log_comb(n_balls, x)) - log_comb(n_balls.sum(), x.sum()) ret = np.exp(ret) return ret proof by picture https://picasaweb.google.com/106983885143680349926/Joepy#5760777856741667746 https://picasaweb.google.com/106983885143680349926/Joepy#5760777861891130962 Cheers, Josef > > Cheers, > > f > _______________________________________________ > NumPy-Discussion mailing list > NumPy-Discussion@scipy.org > http://mail.scipy.org/mailman/listinfo/numpy-discussion _______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
