[R] Problem with Integrate
Peter Dalgaard
p.dalgaard at biostat.ku.dk
Tue Mar 2 17:37:38 CET 2004
"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
More information about the R-help
mailing list