This came up on r-help but indicates a bug. dnbinom(x,n,p) calls dbinom_raw(n-1,...) which returns 0 for n<1. -thomas ---------- Forwarded message ---------- Date: Thu, 08 Feb 2001 17:10:23 +0000 From: Yudi Pawitan <yudi@stat.ucc.ie> To: Mark Myatt <mark@myatt.demon.co.uk> Cc: R-Help <r-help@stat.math.ethz.ch> Subject: Re: [R] Goodness of fit to Poisson / NegBinomial Sorry, Mark, the program worked for Rw1011, but this is what happened in Rw1021:> dnbinom(0:5,.9,.4)[1] 0.4383833 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 #wrong!! while, the correct result from Rw1011 is:> dnbinom(0:5,.9,.4)[1] 0.43838329 0.23672698 0.13493438 0.07826194 0.04578323 0.02692054 I have added a corrected dnbinom() below, which makes the program works for Rw1021. -YP- x_ c(70, 38, 17, 10, 9, 6) kk_ 0:5 n_ sum(x) dnbinom_ function(x,n,p){ gamma(x+n)/gamma(n)/gamma(x+1)*p^n * (1-p)^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))) ==========================================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 > > >Yudi Pawitan yudi@stat.ucc.ie Department of Statistics UCC Cork, Ireland Ph 353-21-490 2906 Fax 353-21-427 1040 -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.- r-help 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-help-request@stat.math.ethz.ch _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._ -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.- 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 _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
>>>>> "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 _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
>>>>> "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 _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._