Hi there
I'm trying to construct a model of mortality risk in 2D space that
requires numerical integration of a hazard function, for which I'm using
the integrate function. I'm occasionally encountering parameter
combinations that cause integrate to terminate with error "Error in
integrate... the integral is probably divergent", which I'm not sure how
to interpret. The problem only crops up for a tiny part of the input
parameter space, but if I can't get this to work for the whole space I'm
a bit stuck! Plotting the integrand shows no sign of obvious problems
such as flatness, extreme spiking or non-finite area, and a tiny tweak
to the input is enough to get the integration to work, even though the
tweaked function looks essentially the same as the one that fails. Below
is some code that demonstrates the problem on my machine running R
2.9.0. Any suggestions on what might be causing this and what I can do
to avoid it would be very gratefully received.
Hoping for some insights
Marcus
##########################
#Problem example
#Works fine:
integrate(hazard, 0,Inf, v=0.1, gam=pi*232/200,
pos=list(r=5,th=atan(3/4)), a=10, b0=10, bt=0.1)
#Gives error "...integral probably divergent":
integrate(hazard, 0,Inf, v=0.1, gam=pi*231/200,
pos=list(r=5,th=atan(3/4)), a=10, b0=10, bt=0.1)
#Plot the integrands - doesn't look obviously problematic
h <-
hazard(seq(0,500,0.1),0.1,pi*231/200,pos=list(r=5,th=atan(3/4)),10,10,0.
1)
plot(seq(0,500,0.1),h,type="l",col=2)
h <-
hazard(seq(0,500,0.1),0.1,pi*232/200,pos=list(r=5,th=atan(3/4)),10,10,0.
1)
lines(seq(0,500,0.1),h)
#Functions used:
#hazard at a given point pos (a list of polar coordinates: distance r
and angle th); a, b0 and bt are model parameters
point.hazard <- function(pos,a,b0,bt) a * exp(-(pos$r^2/(2*b0^2))) *
exp(-(pos$th^2/(2*bt^2)))
#point.hazard for a point related to input point pos by time t and speed
v
hazard <- function(t, v, gam, pos, a, b0, bt)
{ pos2 <- zredef(pos,-t*v,gam)
pos2$th[pos2$th>pi] <- 2*pi-pos2$th[pos2$th>pi]
point.hazard(pos2,a,b0,bt)
}
#Returns a list of polar co-ordinates for a point defined by distance m
in direction gam from starting point pos
zredef <- function(pos, m, gam)
{ x <- pos$r*sin(pos$th) + m*sin(gam)
y <- pos$r*cos(pos$th) + m*cos(gam)
r <- (x^2 + y^2)^0.5
th <- atan(x/y)
th[x==0] <- 0
th[y<0] <- th[y<0] + pi
th[x<0 & y>0] <- th[x<0 & y>0] + 2*pi
list(r=r,th=th)
}
The Zoological Society of London is incorporated by Royal Charter
Principal Office England. Company Number RC000749
Registered address:
Regent's Park, London, England NW1 4RY
Registered Charity in England and Wales no. 208728
_________________________________________________________________________
This e-mail has been sent in confidence to the named add...{{dropped:20}}
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
