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

About this list Date view Thread view Subject view Author view Attachment view

From: Peter Dalgaard (p.dalgaard@biostat.ku.dk)
Date: Wed 05 May 2004 - 19:42:35 EST


Message-id: <x2wu3r9pmc.fsf@biostat.ku.dk>


"Simon Cullen" <cullens@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@biostat.ku.dk)             FAX: (+45) 35327907

______________________________________________ R-help@stat.math.ethz.ch mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


About this list Date view Thread view Subject view Author view Attachment view

This archive was generated by hypermail 2.1.3 : Mon 31 May 2004 - 23:05:07 EST