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 ______________________________________________ [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