[R] Problem with Integrate

[Thomas Lumley]

On Tue, 2 Mar 2004, Anon. wrote:

> The background: I'm trying to fit a Poisson-lognormal distrbutuion to
> some data.  This is a way of modelling species abundances:
> N ~ Pois(lam)
> log(lam) ~ N(mu, sigma2)
> The number of individuals are Poisson distributed with an abundance
> drawn from a log-normal distrbution.
>
> To fit this to data, I need to integrate out lam.  In principle, I can
> do it this way:
>
> PLN1 <- function(lam, Count, mu, sigma2) {
>    dpois(Count, exp(lam), log=F)*dnorm(LL, mu, sqrt(sigma2))
> }
>
> and integrate between -Inf and Inf.  For example, with mu=2, and
> sigma2=2.8 (which are roughly right for the data), and Count=73, I get this:
>
>  > integrate(PLN1, -10, 10, Count=73, mu=2, sigma2=2.8)
> 0.001289726 with absolute error < 2.5e-11
>  > integrate(PLN1, -20, 20, Count=73, mu=2, sigma2=2.8)
> 0.001289726 with absolute error < 2.5e-11
>  > integrate(PLN1, -100, 100, Count=73, mu=2, sigma2=2.8)
> 2.724483e-10 with absolute error < 5.3e-10
>  > integrate(PLN1, -500, 500, Count=73, mu=2, sigma2=2.8)
> 1.831093e-73 with absolute error < 3.6e-73
>  > integrate(PLN1, -1000, 1000, Count=73, mu=2, sigma2=2.8)
> Error in integrate(PLN1, -1000, 1000, Count = 73, mu = 2, sigma2 = 2.8):
>          non-finite function value
> In addition: Warning message:
> NaNs produced in: dpois(x, lambda, log)
>
> So, the integral gets smaller, and then gives an error.
<snip>
> I assume that this is because for much of the range, the integral is
> basically zero.

The help page for integrate() says

     When integrating over infinite intervals do so explicitly, rather
     than just using a large number as the endpoint.  This increases
     the chance of a correct answer - any function whose integral over
     an infinite interval is finite must be near zero for most of that
     interval.

That is, if you want an integral from 0 to Inf, do that.

	-thomas