Hope this helps. spencer graves
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.
I then tried entering the formula directly:
PLN2 <- function(LL, Count, mu, sigma2) {
exp(-LL-(log(LL)-mu)^2/(2*sigma2))*LL^(Count-1)/(gamma(Count+1)*sqrt(2*pi*sigma2))
}
> integrate(PLN2, 0, 100, Count=73, mu=2, sigma2=2.8) 0.001287821 with absolute error < 2.6e-10 > integrate(PLN2, 0, 1000, Count=73, mu=2, sigma2=2.8) 0.001289726 with absolute error < 2.9e-08 > integrate(PLN2, 0, 10000, Count=73, mu=2, sigma2=2.8) 0.001289726 with absolute error < 9.7e-06 > integrate(PLN2, 0, 19100, Count=73, mu=2, sigma2=2.8) 1.160307e-08 with absolute error < 2.3e-08 > integrate(PLN2, 0, 19200, Count=73, mu=2, sigma2=2.8) Error in integrate(PLN2, 0, 19200, Count = 73, mu = 2, sigma2 = 2.8) : non-finite function value
And the same thing happens.
I assume that this is because for much of the range, the integral is basically zero.
Can anyone suggest a fix? Preferably one that will work with Count=320 and Count=0 (both of which I have in the data).
Bob
______________________________________________ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
