R help - May 2004 - Discontinuities in a simple graph (machine precision?)

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.

-- 
SC

Simon Cullen
Room 3030
Dept. Of Economics
Trinity College Dublin

Ph. (608)3477
Email cullens at tcd.ie

Note, this is about integrate, nothing else.  You are looking for an
answer with high absolute precision for larger a because of your divisor
(about 1e-3). If I add rel.tol=1e-12 it is fine.


On Wed, 5 May 2004, Simon Cullen wrote:
> 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.
> 
> 
-- 
Brian D. Ripley,                  ripley at stats.ox.ac.uk
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford,             Tel:  +44 1865 272861 (self)
1 South Parks Road,                     +44 1865 272866 (PA)
Oxford OX1 3TG, UK                Fax:  +44 1865 272595

"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

Simon Cullen wrote:
> 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)
It's a matter of numerical accuracy, you might want to use, e.g.:

     inted <- integrate(integrand, a, Inf,
         rel.tol = .Machine$double.eps^0.5)

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

Thanks to Prof. Ripley, Peter Dalgaard and Uwe Ligges for pointing out  
that all I needed to do to improve the precision was
replace
> inted <- integrate(integrand, a, Inf)
with
> inted <- integrate(integrand, a, Inf,rel.tol = .Machine$double.eps^0.5)
Thanks for the advice!

-- 
SC

Simon Cullen
Room 3030
Dept. Of Economics
Trinity College Dublin

Ph. (608)3477
Email cullens at tcd.ie