[R] Problem with Integrate

[Peter Dalgaard]

"Anon." <bob.ohara at helsinki.fi> writes:

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

It does seem to help a bit if you modify the integrand to

> PLN1
function(lam, Count, mu, sigma2) {
   t1 <- dpois(Count, exp(lam), log=F)
   t2 <- dnorm(lam, mu, sqrt(sigma2))
   ifelse(t1==0|t2==0,0,t1*t2)
}

> integrate(PLN1, -Inf, Inf, Count=73, mu=2, sigma2=2.8)
0.001289726 with absolute error < 6.5e-05
Warning messages:
1: NaNs produced in: dpois(x, lambda, log)
2: NaNs produced in: dpois(x, lambda, log)


For the generic problem, you might need to expand around the maximum,
e.g.

# -log(PLN1) -- better behaved for optimizing
pln2 <-function(lam, Count, mu, sigma2) {
   t1 <- dpois(Count, exp(lam), log=T)
   t2 <- dnorm(lam, mu, sqrt(sigma2), log=T)
   ifelse(t1==0|t2==0,0,-t1-t2)
}

xmax <- nlm(pln2,0, Count=73, mu=2, sigma2=2.8)$estimate # or optimize()

r1 <- uniroot(function(x,...)(PLN1(x,...)-1e-30),c(-100,xmax), 
              Count=73, mu=2, sigma2=2.8)$root
r2 <- uniroot(function(x,...)(PLN1(x,...)-1e-30),c(xmax,100), 
              Count=73, mu=2, sigma2=2.8)$root

integrate(PLN1, r1, r2, Count=73, mu=2, sigma2=2.8)

-- 
   O__  ---- Peter Dalgaard             Blegdamsvej 3  
  c/ /'_ --- Dept. of Biostatistics     2200 Cph. N   
 (*) \(*) -- University of Copenhagen   Denmark      Ph: (+45) 35327918
~~~~~~~~~~ - (p.dalgaard at biostat.ku.dk)             FAX: (+45) 35327907