[R] Discontinuities in a simple graph (machine precision?)

[Peter Dalgaard]

"Simon Cullen" <cullens at tcd.ie> writes:

> Hi,
> 
> I've got an ugly but fairly simple function:
> 
> mdevstdev <- function(a){
> 	l <- dnorm(a)/(1-pnorm(a))
> 	integrand  <- function(z)(abs(z-l)*dnorm(z))
> 	inted <- integrate(integrand, a, Inf)
> 	inted[[1]]/((1- pnorm(a))*sqrt((1 + a*l - l^2)))
> }
> 
> I wanted to quickly produce a graph of this over the range [-3,3] so I
> used:
> 
> plotit <-function(x=seq(-3,3,0.01),...){
> 	y<-sapply(x,mdevstdev)
> 	plot(x,y,...)
> }
> 
> > plotit()
> 
> This produces the graph, but some discontinuities appear on it. I've
> produced the same graph in Mathematica 5
> (http://econserv2.bess.tcd.ie/cullens/R/DOverDelta.eps), and it was
> smooth  over this range (it takes ages to run, though). Is this a
> numerical  precision problem? Any suggestions on how to improve the
> precision?
> 
> I'm running R 1.9 on WinXP, PIII. I haven't changed any R parameters
> that  I know of.

This is probably related to integrating a non-smooth function across
the singularity. It works better like this:

> mdevstdev
function(a){
  l <- dnorm(a)/(1-pnorm(a))
  integrand  <- function(z)(abs(z-l)*dnorm(z))
  inted <- if (l < a) integrate(integrand, a, Inf)[[1]] else
    integrate(integrand, a, l)[[1]] + integrate(integrand, l, Inf)[[1]]
  inted/((1- pnorm(a))*sqrt((1 + a*l - l^2)))
}

(as you see, I was too lazy to think about whether l >= a always...)


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