[Rd] dnbinom(,size<1,)=0 (PR#842)

[tlumley@u.washington.edu]

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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._