[R] Problem with Integrate

Anon. bob.ohara at helsinki.fi
Tue Mar 2 16:23:23 CET 2004


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

-- 
Bob O'Hara

Dept. of Mathematics and Statistics
P.O. Box 4 (Yliopistonkatu 5)
FIN-00014 University of Helsinki
Finland
Telephone: +358-9-191 23743
Mobile: +358 50 599 0540
Fax:  +358-9-191 22 779
WWW:  http://www.RNI.Helsinki.FI/~boh/
Journal of Negative Results - EEB: http://www.jnr-eeb.org




More information about the R-help mailing list