[SciPy-user] looking for a negative binomial distribution

[Robert Kern]

On Sat, Aug 23, 2008 at 02:25, Chris Fonnesbeck <listservs at mac.com> wrote:
> Chris Fonnesbeck <listservs <at> mac.com> writes:
>> Robert Kern <robert.kern <at> gmail.com> writes:
>> > Fixed in numpy SVN. This will be part of the upcoming 1.2.0 release
>> > and probably 1.1.2 if we do one.
>> >
>>
>> Thanks Robert,
>>
>> I would have poked around the source code myself, but I'm on a deadline. I will
>> update and test the code now though.
>>
>> The easiest way to generate NB random variables is to sample from a gamma,
>> then use that value as the mean for a sample from a poisson.

Yes, that's exactly what I'm doing.

long rk_negative_binomial(rk_state *state, double n, double p)
{
    double Y;

    Y = rk_gamma(state, n, (1-p)/p);
    return rk_poisson(state, Y);
}

> The fix appears to work, though the samples appear to be biased high relative to
> thosegenerated using another method, based on a couple tests. I'll make sure
> I'm not making an error first, before complaining again. The first parameter does
> accept non-integer arguments now, however.

I cannot see any particular deviation from the CDF or PDF (via Q-Q
plots and histograms). The sample means appear to match the
theoretical mean over a fairly wide range of parameters (see below).
If you are still having trouble getting the distribution to match your
other method, please let me know with details about your other method
and the results that numpy generates for you.

In [132]: def meandiff(r, mean):
   .....:     p = float(r) / (r+mean)
   .....:     return np.random.negative_binomial(r,p,size=10000).mean() - mean
   .....:

-- 
Robert Kern

"I have come to believe that the whole world is an enigma, a harmless
enigma that is made terrible by our own mad attempt to interpret it as
though it had an underlying truth."
 -- Umberto Eco