[R] followup: lookup function for density(...) objects

[Prasad, Rajiv]

Thanks to Ross Ihaka and Bob Wheeler for responding to my earlier question.  I
looked into the Johnson system functions in SuppDists package.  For now, I want
to stick with the density(...) estimator, and so still need the variate lookup
function.

As per Ross' suggestion, I just did a numerical integration on the density
object and used approxfun/splinefun to "lookup" the variate value for specified
cumulative probability.  Here's the function.  Comments are welcome.

--------
#
# return the variate(s) for specified cumulative
# probability(ies) for a density(...) object
#
# Rajiv Prasad (rajiv.prasad at pnl.gov) 05/14/2001
#
RPDensityLookup <- function(density.obj, p=0.5, interp="l")
{
  if(missing(density.obj))
  {
    cat("\nUsage: RPDensityLookup(density.obj, p=0.5, interp=\"l\")")
    cat("\n  density.obj: an object returned from the density(...) function")
    cat("\n            p: cumulative probability value(s) for which variate")
    cat("\n               values are needed")
    cat("\n       interp: \"l\" for linear interpolation between data points")
    cat("\n               \"s\" for spline fit using data points\n\n")
  }

  n <- length(density.obj$x)
  dx <- density.obj$x[2:n] - density.obj$x[1:(n-1)]
  x <- rep(NA, length(p))

  # midpoints
  cumx <- (density.obj$x[2:n] + density.obj$x[1:(n-1)]) / 2
  # numerical integration
  cumy <- cumsum((density.obj$y[2:n] + density.obj$y[1:(n-1)]) / 2 * dx)

  if(interp == "l")
  {
    Fx <- approxfun(cumy, cumx)  # linear interpolation function
    x <- Fx(p)                   # the variates corresponding to p's
  }
  else if(interp == "s")
  {
    Fx <- splinefun(cumy, cumx)  # spline fit
    x <- Fx(p)                   # the variates corresponding to p's
  }
  else
  {
    cat(paste("\nError: unrecognized interp method (\"", interp, "\").\n\n",
sep=""))
  }

  x
}
--------

Example:

> x _ rnorm(1000)
> d.x _ density(x)
> p _ seq(0.001,0.999,by=0.001)
> x.l _ RPDensityLookup(d.x, p, "l")
> x.s _ RPDensityLookup(d.x, p, "s")
> hist(x,probability=T,ylim=c(0,1)); box()
> lines(d.x,col="blue")
> lines(x.l,p,col="green")
> lines(x.s,p,col="red")
> RPDensityLookup(d.x, 0.5, "l")
[1] -0.001036816
> RPDensityLookup(d.x, 0.5, "s")
[1] -0.001036673
> RPDensityLookup(d.x, 0.95, "l")
[1] 1.751115
> RPDensityLookup(d.x, 0.95, "s")
[1] 1.751107


Rajiv.
--------
Rajiv Prasad
Postdoctoral Research Associate, Hydrology Group
Pacific Northwest National Laboratory, P.O. Box 999, MSIN K9-33
Richland, WA 99352
Voice: (509) 375-2096  Fax: (509) 372-6089  Email: rajiv.prasad at pnl.gov


-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._