**From:** Anon. (*bob.ohara@helsinki.fi*)

**Date:** Tue 02 Mar 2004 - 19:55:14 EST

**Next message:**Jean-Noel: "[R] Import range of cells from Excel"**Previous message:**Prof Brian Ripley: "Re: [R] passing a string from .C()"**Next in thread:**Spencer Graves: "Re: [R] Problem with Integrate"**Reply:**Spencer Graves: "Re: [R] Problem with Integrate"**Reply:**Thomas Lumley: "Re: [R] Problem with Integrate"

Message-id: <40444BF2.9050900@helsinki.fi>

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'HaraDept. 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

______________________________________________ R-help@stat.math.ethz.ch 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

**Next message:**Jean-Noel: "[R] Import range of cells from Excel"**Previous message:**Prof Brian Ripley: "Re: [R] passing a string from .C()"**Next in thread:**Spencer Graves: "Re: [R] Problem with Integrate"**Reply:**Spencer Graves: "Re: [R] Problem with Integrate"**Reply:**Thomas Lumley: "Re: [R] Problem with Integrate"

*
This archive was generated by hypermail 2.1.3
: Fri 23 Apr 2004 - 18:36:27 EST
*