[Numpy-discussion] Multivariate hypergeometric distribution?
josef.pktd at gmail.com
josef.pktd at gmail.com
Mon Jul 2 23:27:56 EDT 2012
On Mon, Jul 2, 2012 at 10:53 PM, Fernando Perez <fperez.net at gmail.com> wrote:
> On Mon, Jul 2, 2012 at 7:49 PM, Fernando Perez <fperez.net at 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 at scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion
More information about the NumPy-Discussion
mailing list