Calculating a spatial distribution from a probability density

[trmcclain]

I'm hoping this will be the last time I call for help, but in any case, here it goes.

I thought I had a handle on this before, but in all of my attempts, my code diverges within a few iterations.
My problem is creating a spatial distribution of particles given a probability density. I've tried to start out with an easier function, the normal distribution f(x)=e^{-x^2}

The method I've tried to implement, which has been failing, is the Newton-Raphson method:
x_{new}=x_{guess}-\frac{U-G(x)}{f(x)} where U is a random number between 0 and 1 and G(x)=\int_{-∞}^x f(x)\,dx

My code is in C++ and is as follows:
dr = rand() % 100;
U=dr/100.;
do{
Gfunc=(1./R)*(sqrt(PI)/2.)*(erf(rg)+1);
func=exp(-(rg*rg));
whole=(U-Gfunc)/func;

rn=rg-whole;
if(rn > 1E-5)
rg=rn;
}while(rn>1E-5);

(Sorry if there's a way to make code pretty on the forums)
**Note, I do use srand(time(NULL)) above this loop.**

 

[pasmith]

I don't think your method works even in theory.

Ultimately what you want to do is produce a sample from the given distribution.

There are ways to take output from rand() and produce a sample from an arbitrary distribution; the relevant section of Numerical Recipes (3rd Edition) is a good place to start. (You can view a limited number of pages of the online edition per month for free; the relevant section starts at page 361.)

The usual way to get a U[0,1] random variable from rand() is:

double u01;

u01 = ((double)rand())/((double)[url=http://www.cplusplus.com/reference/cstdlib/RAND_MAX/]RAND_MAX[/url]);

 

[Stephen Tashi]

the normal distribution f(x)=e^{-x^2}

That isn't a probability density unless you put the correct factors on it (which apparently you did try to do in your code).

Your statement of what you are trying to do isn't clear. I think you are trying to write a function to compute the inverse of a cumulative probability distribution (not "generate a distribution").

So your problem is:

Solve for x in the equation v - \int_{-\infty}^x f(s) ds = 0

where v is a given probability and f(s) is a given probabiity density.