[Numpy-discussion] Multivariate hypergeometric distribution?

josef.pktd@gmai... josef.pktd@gmai...
Mon Jul 2 22:27:56 CDT 2012


On Mon, Jul 2, 2012 at 10:53 PM, Fernando Perez <fperez.net@gmail.com> wrote:
> On Mon, Jul 2, 2012 at 7:49 PM, Fernando Perez <fperez.net@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


More information about the NumPy-Discussion mailing list