[R] Creating q and p functions from a self-defined distribution

[Eli Gurarie]

Hello all,

I am fishing for some suggestions on efficient ways to make qdist and 
pdist type functions from an arbitrary distribution whose probability 
density function I've defined myself.

For example, let's say I have a distribution whose pdf is:

dRN <- function(x,d,v,s)
# d, v, and s are parameters
    return(d/x^2/sv/sqrt(2*pi)*exp(-(d-v*x)^2/2/(sv^2*x^2)))

this is a legitimate distribution over the reals (though it has a 
singularity at x=0)

at the moment, to get a "p" value, I sum the dRN over some arbitrary 
interval dt:

pRN<-function(x,d,v,s,dt=.001)
{
   xs<-seq(0,x,dt)
   dRN<-dRN(xs,d,v,s)
   return(sum(dRN)*dt)
}

which is OK, except I sometimes want a vector of p values for a vector 
of x's, say: pRN(xs=seq(0,10,.1) ... ) which involves an unwieldly loop.

Getting a quantile value is a real nightmare!  Right now I have a thing 
that involves using approx() and matching rounded numbers and requires a 
loop and is slow.

It seems surprising that it would be so hard to invert a function that 
is perfectly defined!

Are there packages/functions/algorithms that allow one to manipulate 
arbitrarily defined distributions?

Any suggestions appreciated,
Thanks,
   Eli




*************************************************
  Eli Gurarie
  Quantitative Ecology and Resource Management
  University of Washington, Seattle