# [Rd] update on dnbinom with large "size"

Martin Maechler maechler at stat.math.ethz.ch
Fri Jul 4 15:23:18 CEST 2008

```>>>>> "BB" == Ben Bolker <bolker at ufl.edu>
>>>>>     on Fri, 04 Jul 2008 08:43:21 -0400 writes:

BB> -----BEGIN PGP SIGNED MESSAGE-----
BB> Hash: SHA1

BB> ~  thanks!

BB> ~  Poisson+warning seems very sensible.

better than now, in the extreme cases,
but far from optimal (in those cases where it is only approximate).
I'm currently working on a much better solution, at least for
[dpr]nbinom();  qnbinom() needs a bit more work, and may have to
wait a bit longer.

Martin

BB> ~  Looks like things start to get wonky around size=5e14
BB> (this is for mu=1)
BB> but thinking a little harder about the numerical analysis
BB> will probably give a more precise criterion.

BB> curve(dnbinom(1,mu=1,size=x),from=5e14,to=1e16,log="x",ylim=c(0.32,0.37))

BB> ~  Ben Bolker

BB> Prof Brian Ripley wrote:
BB> | On Fri, 4 Jul 2008, Martin Maechler wrote:
BB> |
BB> |>>>>>> "BDR" == Prof Brian Ripley <ripley at stats.ox.ac.uk>
BB> |>>>>>>     on Fri, 4 Jul 2008 09:02:21 +0100 (BST) writes:
BB> |>
BB> |>    BDR> On Thu, 3 Jul 2008, Ben Bolker wrote:
BB> |>    >> -----BEGIN PGP SIGNED MESSAGE----- Hash: SHA1
BB> |>    >>
BB> |>    >>
BB> |>    >> ~ turns out I don't need to look at the C code.
BB> |>    >>
BB> |>    >> ~ if one uses the mu/size parameterization of the
BB> |>    >> negative binomial, R computes size/(size+mu) to switch
BB> |>    >> parameterizations.  If size>>mu this gets rounded to 1
BB> |>    >> ... should be easy enough to test and return NA under
BB> |>    >> these circumstances?
BB> |>
BB> |>    BDR> It is all vectorized, so not so easy.  But why is NA
BB> |>    BDR> appropriate, when a Poisson approximation seems more
BB> |>    BDR> appropriate?
BB> |>
BB> |> yes, definitely more appropriate.
BB> |>
BB> |>    BDR> The same issue affects [pqr]nbinom, of course.
BB> |>
BB> |> indeed.
BB> |>
BB> |>    BDR> The short-term advice is of course 'don't do that'
BB> |> :-)
BB> |>
BB> |>    BDR> but use a suitable approximation such as Poisson.  Maybe
BB> |>    BDR> there are reasons to support values of size > 1e10, but
BB> |>    BDR> I suspect only for completeness and maybe it is better
BB> |>    BDR> to at least give a warning.
BB> |>
BB> |> As you know I'm very much in favor of "completeness" here (and
BB> |> have already spent an unreasonable/unpayable amount of hours to
BB> |> achieve that in quite a few cases for R in the past).
BB> |>
BB> |> I'd much prefer a solution without a warning that used the
BB> |> appropriate approximation when needed.
BB> |
BB> | The argument for warning is that one possible reason for doing this is
BB> | to explore the difference between negative binomial and Poisson, and so
BB> | it would be helpful to be warned that a Poisson had been used (as we
BB> | warn in some other cases where accuracy is lost, if not entirely
BB> | consistently).
BB> |
BB> |> I'm now between a workshop and a week of vacation, so I am not
BB> |> promising anything, but I'd plan to get this fixed before the
BB> |> end of the summer.
BB> |>
BB> |>    >> - --------------------------
BB> |>    >> ~ Check it out:
BB> |>    >>
BB> |>    >> curve(dnbinom(1,mu=0.5,size=x),log="x",from=1,to=1e18)
BB> |>    >> abline(h=dpois(1,lambda=0.5),col=2,lty=2)
BB> |>    >> text(1,dpois(1,lambda=0.5)+0.02,"Poisson",col=2,pos=4)
BB> |>
BB> |> a very nice example, thank you, Ben!
BB> |>
BB> |> Regards, Martin
BB> |>
BB> |>    >>
BB> |>    >> ~ I will take a look in the C code when I get a chance to
BB> |>    >> see if I can offer a patch, but in the meantime wanted to
BB> |>    >> alert people to this "feature" ...  (it looks like I will
BB> |>    >> have to go through dnbinom and dbinom_raw to see where
BB> |>    >> the problem is ...)
BB> |>    >>
BB> |>    >> ~ Ben Bolker
BB> |>
BB> |

BB> -----BEGIN PGP SIGNATURE-----
BB> Version: GnuPG v1.4.6 (GNU/Linux)
BB> Comment: Using GnuPG with Mozilla - http://enigmail.mozdev.org

BB> iD8DBQFIbhrpc5UpGjwzenMRAmCsAJ91HwBzzl29NiJ+LyqpMEr9ROC2DgCeNoMC
BB> n+2CjBTgBn7uPXfr2lo5ydc=
BB> =Kcj0
BB> -----END PGP SIGNATURE-----

```