[Rd] dnbinom(,size<1,)=0 (PR#842)
maechler@stat.math.ethz.ch
maechler@stat.math.ethz.ch
Fri, 9 Feb 2001 09:08:58 +0100 (MET)
>>>>> "MM" == Martin Maechler <maechler@stat.math.ethz.ch> writes:
>>>>> "TL" == Thomas Lumley <tlumley@u.washington.edu> writes:
TL> This came up on r-help but indicates a bug.
TL> dnbinom(x,n,p) calls dbinom_raw(n-1,...)
TL> which returns 0 for n<1.
MM> Actually for all non-integer n.
wrong, Martin!
Thomas was right:
dbinom_raw(m, *)
seems fine for non-negative non-integers m.
MM> This is now fixed for "R-release" aka "R-patched" aka "R 1.2.2 to be" :
MM> --- src/nmath/dnbinom.c.~6~ Wed Jan 24 18:02:07 2001
MM> +++ src/nmath/dnbinom.c Fri Feb 9 08:37:14 2001
MM> @@ -47,6 +47,12 @@
MM> if (x < 0 || !R_FINITE(x)) return R_D__0;
MM> x = R_D_forceint(x);
MM> + if(R_D_nonint(n)) {
MM> + return R_D_exp(lfastchoose(x + n - 1, x)
MM> + + n * log(p) + x * log(1 - p));
MM> + }
MM> + else { /* n is integer */
MM> prob = dbinom_raw(n-1,n+x-1,p,1-p,give_log);
MM> return((give_log) ? log(p) + prob : p*prob);
MM> + }
MM> }
The above fix reverts to pre R 1.2.0 behaviour
for non-integer `n' (in C; `size' in R) which is `first order ok',
but far from optimal.
The real solution must be better,
and should be done along the code (see comments!) in dbeta.c.
I won't have time to look at this for the next time though.
MM> Martin Maechler <maechler@stat.math.ethz.ch> http://stat.ethz.ch/~maechler/
MM> Seminar fuer Statistik, ETH-Zentrum LEO D10 Leonhardstr. 27
MM> ETH (Federal Inst. Technology) 8092 Zurich SWITZERLAND
MM> phone: x-41-1-632-3408 fax: ...-1228 <><
TL> ---------- Forwarded message ----------
TL> Date: Thu, 08 Feb 2001 17:10:23 +0000
TL> From: Yudi Pawitan <yudi@stat.ucc.ie>
TL> To: Mark Myatt <mark@myatt.demon.co.uk>
TL> Cc: R-Help <r-help@stat.math.ethz.ch>
TL> Subject: Re: [R] Goodness of fit to Poisson / NegBinomial
TL> Sorry, Mark, the program worked for Rw1011, but this is
TL> what happened in Rw1021:
>>> dnbinom(0:5,.9,.4)
TL> [1] 0.4383833 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 #wrong!!
TL> while, the correct result from Rw1011 is:
>>> dnbinom(0:5,.9,.4)
TL> [1] 0.43838329 0.23672698 0.13493438 0.07826194 0.04578323 0.02692054
TL> I have added a corrected dnbinom() below, which makes the program
TL> works for Rw1021.
TL> -YP-
TL> x_ c(70, 38, 17, 10, 9, 6)
TL> kk_ 0:5
TL> n_ sum(x)
TL> dnbinom_ function(x,n,p){
TL> gamma(x+n)/gamma(n)/gamma(x+1)*p^n * (1-p)^x
TL> }
TL> fn_ function(p){
TL> r_ p[1]
TL> th_ p[2]
TL> ll _ -sum(x[1:max(kk)]*log(dnbinom(kk[-length(kk)],r,th))) -
TL> x[length(kk)]*log(1-pnbinom((max(kk)-1),r,th))
TL> return(ll)
TL> }
TL> a_ nlm(fn,p=c(.9,.5),hes=T)
TL> est_ a$est
TL> phat_ c(dnbinom(kk[-length(kk)],est[1],est[2]),
TL> 1-pnbinom((max(kk)-1),est[1],est[2]))
TL> E_ n*phat
TL> chi2_ sum((x-E)^2/E)
TL> cat('Chi2 goodness-of-fit = ',chi2,'\n')
TL> cat('P-value=',1-pchisq(chi2,df=length(kk)-1-2),'\n')
TL> print(cbind(O=x,E=round(E,1),err=round((x-E)^2/E,1)))
TL> ===========================================
TL> At 02:31 PM 2/8/01 +0000, Mark Wyatt wrote:
>>> Yudi,
>>>
>>> Thanks for this:
>>>
>>>> x_ c(70, 38, 17, 10, 9, 6)
>>>> kk_ 0:5
>>>> n_ sum(x)
>>>>
>>>> fn_ function(p){
>>>> r_ p[1]
>>>> th_ p[2]
>>>> ll _ -sum(x[1:max(kk)]*log(dnbinom(kk[-length(kk)],r,th))) -
>>>> x[length(kk)]*log(1-pnbinom((max(kk)-1),r,th))
>>>> return(ll)
>>>> }
>>>>
>>>> a_ nlm(fn,p=c(.9,.5),hes=T)
>>>> est_ a$est
>>>>
>>>> phat_ c(dnbinom(kk[-length(kk)],est[1],est[2]),
>>>> 1-pnbinom((max(kk)-1),est[1],est[2]))
>>>> E_ n*phat
>>>> chi2_ sum((x-E)^2/E)
>>>> cat('Chi2 goodness-of-fit = ',chi2,'\n')
>>>> cat('P-value=',1-pchisq(chi2,df=length(kk)-1-2),'\n')
>>>>
>>>> print(cbind(O=x,E=round(E,1),err=round((x-E)^2/E,1)))
>>>
>>> It is certainly more elegant than anything I would have come up with but
>>> I just cannot get it to work (Rwin 1.21). The call to nlm() returns:
>>>
>>> Warning messages:
>>> 1: NA/Inf replaced by maximum positive value
>>> 2: NA/Inf replaced by maximum positive value
>>> 3: NA/Inf replaced by maximum positive value
>>> 4: NA/Inf replaced by maximum positive value
>>> 5: NA/Inf replaced by maximum positive value
>>> 6: NA/Inf replaced by maximum positive value
>>> 7: NA/Inf replaced by maximum positive value
>>> 8: NA/Inf replaced by maximum positive value
>>>
>>> est holds the starting parameters:
>>>
>>> [1] 0.9 0.5
>>>
>>> E (from n*phat) holds:
>>>
>>> [1] 80.38301 0.00000 0.00000 0.00000 0.00000 3.90972
>>>
>>> which gives a very big chi-square.
>>>
>>> Sorry to be a nuisance.
>>>
>>> Mark
>>>
>>>
>>> --
>>> Mark Myatt
>>>
>>>
>>>
TL> Yudi Pawitan yudi@stat.ucc.ie
TL> Department of Statistics UCC
TL> Cork, Ireland
TL> Ph 353-21-490 2906
TL> Fax 353-21-427 1040
TL> -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
TL> r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
TL> Send "info", "help", or "[un]subscribe"
TL> (in the "body", not the subject !) To: r-help-request@stat.math.ethz.ch
TL> _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
TL> -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
TL> r-devel mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
TL> Send "info", "help", or "[un]subscribe"
TL> (in the "body", not the subject !) To: r-devel-request@stat.math.ethz.ch
TL> _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
MM> -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
MM> r-devel mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
MM> Send "info", "help", or "[un]subscribe"
MM> (in the "body", not the subject !) To: r-devel-request@stat.math.ethz.ch
MM> _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-devel mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !) To: r-devel-request@stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._