Thu 01 Aug 2013 02:02:57 PM UTC, comment #8:
Hi,
why not to change the relevant part in nbinpdf.m to
according to the definition given in comment #2?
About bincoeff.m, I do not understand the formula given on line 92. Where does it come from?
Marco

Thu 01 Aug 2013 10:54:13 AM UTC, comment #7:
Thanks for your help with this issue. I believe the R statistical software might be helpful.
http://stat.ethz.ch/Rmanual/Rpatched/library/stats/html/NegBinomial.html
Following the source code, I see:
http://svn.rproject.org/R/trunk/src/nmath/pnbinom.c
http://svn.rproject.org/R/trunk/src/nmath/pbeta.c
http://svn.rproject.org/R/trunk/src/nmath/toms708.c
(bratio function for incomplete beta)
Maybe there is something in there that could help you patch the discontinuities.

Thu 01 Aug 2013 03:18:17 AM UTC, comment #6:
Ok, I see the point you are making. There are similar discontinuities at other integer values of N +/ a few epsilon. Visually,

Thu 01 Aug 2013 01:48:15 AM UTC, comment #5:
We obviously can't just compute the binomial coefficient using the gamma function for all arguments, since it will fail for everything but the tiniest numbers. If you look at the actual implementation, you'll see that we're using bincoeff to compute this, which in turn computes the binomial coefficients using the usual loggamma approach.
The bug, if there is one, is in the bincoeff function then. And again, the problem is that we are essentially computing the value of the beta function around a removable singularity, using sums of loggamma functions. You'll note that on line 92 here,
http://hg.savannah.gnu.org/hgweb/octave/file/26589abbc78d/scripts/miscellaneous/bincoeff.m#l91
since we've determined that n is not an integer, then we can take the lgamma of it. But then we're actually very close to a singularity of the gamma function. The whole thing still sort of manages to cancel and not give a completely ridiculous answer, but I don't know if there is any sort of code that we can implement here to make it better.
Do you have a suggestion?

Thu 01 Aug 2013 12:47:47 AM UTC, comment #4:
I understand about floating point issues and whether something is an integer or not, but if I want to nbinpdf for real arguments, then I would expect continuity.
octave:2> nbinpdf(0, 3, .5)
ans = 0.12500
octave:3> nbinpdf(0, 2.9999999999999, .5)
ans = 0.12495
octave:4> nbinpdf(0, 1/(12/3), .5)
ans = 0.18110 # discontinuous by bad floating pt luck
octave:5> nbinpdf(0, 3.0000000000001, .5)
ans = 0.12475
By definition in terms of the gamma function,
nbinpdf(x,n,p) = gamma(x+n) / (gamma(x+1)gamma(n)) p^n*(1p)^x
octave:10> x=0; n=1/(12/3); p=.5; gamma(x+n) / (gamma(x+1)gamma(n)) p^n*(1p)^x
ans = 0.12500
octave:11> x=0; n=1/(12/3); p=.5; gamma(x+n) / (gamma(x+1)gamma(n)) p^n*(1p)^x
ans = 0.12500 # gives the right answer

Wed 31 Jul 2013 11:53:02 PM UTC, comment #3:
No, if that's your reasoning then there is no bug here. The problem is that you think your junk variable is still exactly equal to 3, which it is not.
The help for nbinpdf clearly states that it will accept values for N that are either integer or noninteger. It is up to the caller to ensure that the value is actually an integer if that's what you wanted to call it with.
If I misunderstood, you can respond and clarify what you think the bug is, but I am closing this as invalid.

Wed 31 Jul 2013 11:42:12 PM UTC, comment #2:
The PMF for negative binomial should be for my example
nbinpdf(0,3,.5)
((0+31) choose 0) * (.5)^3 * (.5)^0 = .125
see http://en.wikipedia.org/wiki/Negative_binomial
Since 1/(12/3) = 1/(1/3) = 3, this should also give the same answer, but it gives .18110 instead.
nbinpdf(0,1/(12/3),.5)
Perhaps the implementation needs to round the second argument to an exact integer?

Wed 31 Jul 2013 07:40:20 PM UTC, comment #1:
Is this a bug? Ultimately it boils down to discontinuities of the beta function at negative integers, where the gamma function has its poles.

Wed 31 Jul 2013 07:10:55 PM UTC, original submission:
> nbinpdf(0,3,.5)
ans = 0.12500 # this is correct
> junk=1/(12/3) # equals 3 in exact arithmetic
junk = 3.0000
> nbinpdf(0,junk,.5)
ans = 0.18110 # this is incorrect
