Hi Marcus, I always use a smaller error tolerance in `integrate' than the default value. I generally use 1.e-07, whereas the default is only about 1.e-04. Sometimes you may also need to increase the number of subdivisions from its default value of 100.
Your problem disappears if you use a smaller tolerance for convergence: > integrate(hazard, 0, Inf, v=0.1, gam=pi*231/200, pos=list(r=5,th=atan(3/4)), rel.tol=1.e-07, a=10, b0=10, bt=0.1) 0.0002290673 with absolute error < 1.9e-11 > Ravi. ---------------------------------------------------------------------------- ------- Ravi Varadhan, Ph.D. Assistant Professor, The Center on Aging and Health Division of Geriatric Medicine and Gerontology Johns Hopkins University Ph: (410) 502-2619 Fax: (410) 614-9625 Email: rvarad...@jhmi.edu Webpage: http://www.jhsph.edu/agingandhealth/People/Faculty_personal_pages/Varadhan.h tml ---------------------------------------------------------------------------- -------- -----Original Message----- From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On Behalf Of Marcus Rowcliffe Sent: Wednesday, September 23, 2009 11:31 AM To: r-help@r-project.org Subject: [R] Numerical integration problem 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:8}} ______________________________________________ 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.
