[R] Problem with Integrate

[Prof Brian Ripley]

I think what you actually want is an expectation with respect to a normal 
random variable.  If so, you should be integrating wrt to a standard 
normal (after a linear change) and not wrt to Lebesgue measure.  That is, 
use something like Gauss-Hermite integration, although I expect 
integrate() might work well enough over the range (-4, 4).

On Tue, 2 Mar 2004, Anon. wrote:

> Thomas Lumley wrote:
> > 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.
> > 
> Sorry, I forgot to put that in my message.  It gives the same error as a 
> large value.  I assume it's all a result of the NaN's being returned.
> 
> Bob
> 
> 

-- 
Brian D. Ripley,                  ripley at stats.ox.ac.uk
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford,             Tel:  +44 1865 272861 (self)
1 South Parks Road,                     +44 1865 272866 (PA)
Oxford OX1 3TG, UK                Fax:  +44 1865 272595