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
--
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
Have you done a search of "www.r-project.org" -> search ->
"R site
search" for "hermite quadrature"? I just got 11 hits on this,
the
second of which referred to "gauss.quad {statmod}". This said,
"See
also ... gauss.quad.prob {statmod}". I haven't tried it, but it looks
like what you want.
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
>
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