From maechler at stat.math.ethz.ch  Fri Jul  4 11:51:02 2008
From: maechler at stat.math.ethz.ch (Martin Maechler)
Date: Fri, 4 Jul 2008 11:51:02 +0200
Subject: [Rd] update on dnbinom with large "size"
In-Reply-To: <alpine.LFD.1.10.0807040851240.21565@gannet.stats.ox.ac.uk>
References: <486D8C99.2050006@ufl.edu>
	<alpine.LFD.1.10.0807040851240.21565@gannet.stats.ox.ac.uk>
Message-ID: <18541.62086.316695.451171@stat.math.ethz.ch>

>>>>> "BDR" == Prof Brian Ripley <ripley at stats.ox.ac.uk>
>>>>>     on Fri, 4 Jul 2008 09:02:21 +0100 (BST) writes:

    BDR> On Thu, 3 Jul 2008, Ben Bolker wrote:
    >> -----BEGIN PGP SIGNED MESSAGE----- Hash: SHA1
    >> 
    >> 
    >> ~ turns out I don't need to look at the C code.
    >> 
    >> ~ if one uses the mu/size parameterization of the
    >> negative binomial, R computes size/(size+mu) to switch
    >> parameterizations.  If size>>mu this gets rounded to 1
    >> ... should be easy enough to test and return NA under
    >> these circumstances?

    BDR> It is all vectorized, so not so easy.  But why is NA
    BDR> appropriate, when a Poisson approximation seems more
    BDR> appropriate?

yes, definitely more appropriate.

    BDR> The same issue affects [pqr]nbinom, of course.

indeed.

    BDR> The short-term advice is of course 'don't do that' 
:-)

    BDR> but use a suitable approximation such as Poisson.  Maybe
    BDR> there are reasons to support values of size > 1e10, but
    BDR> I suspect only for completeness and maybe it is better
    BDR> to at least give a warning.

As you know I'm very much in favor of "completeness" here (and
have already spent an unreasonable/unpayable amount of hours to
achieve that in quite a few cases for R in the past).

I'd much prefer a solution without a warning that used the
appropriate approximation when needed.
I'm now between a workshop and a week of vacation, so I am not
promising anything, but I'd plan to get this fixed before the
end of the summer.

    >> - --------------------------
    >> ~ Check it out:
    >> 
    >> curve(dnbinom(1,mu=0.5,size=x),log="x",from=1,to=1e18)
    >> abline(h=dpois(1,lambda=0.5),col=2,lty=2)
    >> text(1,dpois(1,lambda=0.5)+0.02,"Poisson",col=2,pos=4)

a very nice example, thank you, Ben!

Regards, Martin

    >> 
    >> ~ I will take a look in the C code when I get a chance to
    >> see if I can offer a patch, but in the meantime wanted to
    >> alert people to this "feature" ...  (it looks like I will
    >> have to go through dnbinom and dbinom_raw to see where
    >> the problem is ...)
    >> 
    >> ~ Ben Bolker



[Rd] update on dnbinom with large "size"
Martin Maechler