[Rd] dnbinom(,size<1,)=0 (PR#842)
maechler@stat.math.ethz.ch
maechler@stat.math.ethz.ch
Fri, 9 Feb 2001 08:59:01 +0100 (MET)
>>>>> "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.
Actually for all non-integer n.
This is now fixed for "R-release" aka "R-patched" aka "R 1.2.2 to be" :
--- src/nmath/dnbinom.c.~6~ Wed Jan 24 18:02:07 2001
+++ src/nmath/dnbinom.c Fri Feb 9 08:37:14 2001
@@ -47,6 +47,12 @@
if (x < 0 || !R_FINITE(x)) return R_D__0;
x = R_D_forceint(x);
+ if(R_D_nonint(n)) {
+ return R_D_exp(lfastchoose(x + n - 1, x)
+ + n * log(p) + x * log(1 - p));
+ }
+ else { /* n is integer */
prob = dbinom_raw(n-1,n+x-1,p,1-p,give_log);
return((give_log) ? log(p) + prob : p*prob);
+ }
}
Martin Maechler <maechler@stat.math.ethz.ch> http://stat.ethz.ch/~maechler/
Seminar fuer Statistik, ETH-Zentrum LEO D10 Leonhardstr. 27
ETH (Federal Inst. Technology) 8092 Zurich SWITZERLAND
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> _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._