From maechler@stat.math.ethz.ch  Fri Feb  9 07:59:01 2001
From: maechler@stat.math.ethz.ch (maechler@stat.math.ethz.ch)
Date: Fri, 9 Feb 2001 08:59:01 +0100 (MET)
Subject: [Rd] dnbinom(,size<1,)=0 (PR#842)
Message-ID: <200102090759.IAA25635@pubhealth.ku.dk>

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


[Rd] dnbinom(,size<1,)=0 (PR#842)
maechler@stat.math.ethz.ch