R help - Mar 2007 - Creating q and p functions from a self-defined distribution

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

On 15-Mar-07 12:09:42, Eli Gurarie wrote:
> 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)
> [...]
> 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?
Do not be surprised! The Normal distribution function itself, with
pdf (1/(sv*sqrt(1*pi)))*exp(-((x - mu)^2)/(2*sv^2)), is perfectly
well defined. Yet the literature of computation over decades has
presented procedure after procedure for computing the cumulative
fucntion, and its inverse, to desired precision. None of these is
simple. Indeed, (to use your own word), the Normal distribution
is a "nightmare" from this point of view.

So being well-defined is no guarantee that computing its p and q
values will be a simple or easy problem. And what kind of method
is good for a particular distribution will depend strongly on
what the distribution is (for instance, whether it has tails
which tend rapidly to 0 like the Normal, whether there are good
asymptotic formulae, etc.).

In the case of the example you give above, the transformation
from x to u = 1/x will translate it into a Normal distribution
problem, after which you can use (circumspectly ... ) the R
functions pnorm and qnorm (which are based on the best from
the above literature) to deal with it.

But you give it only as an example ... and you are asking for
methods for a general user-defined distribution. For the reasons
given above, a good general-purpose method is unlikely to exist.

With best wishes,
Ted.

--------------------------------------------------------------------
E-Mail: (Ted Harding) <Ted.Harding at manchester.ac.uk>
Fax-to-email: +44 (0)870 094 0861
Date: 15-Mar-07                                       Time: 13:26:27
------------------------------ XFMail ------------------------------