15,918,333 members
See more:
Hi,

I am trying to write some code to run a Monte Carlo simulation. I am stuck on the first bit which requires the generation of input values to the calculations. If each input value has a probability distribution, for each run I need to pick values from this distribution by choosing an area value between 0 and 1 each time, where the area under the distribution is normalized to be 1.

Can anyone help me? I have been hunting for several hours and have not really got anywhere. I have looked at some examples in numerical recipes in C but am not sure if what I am looking at is what I need.

I have the following thoughts, which may or may not have anything correct in them:

If I have a probability distribution f(x), where x goes from xmin to xmax, then f(x)dx over all x = 1.

If I pick a random area value, e.g. 0.3, I need to find the value of x, xval, which is the point at which the area under the curve is 0.3.
So I need to integrate f(x) from xmin to xval and equate that to 0.3 in order to find xval

Eg if f(x) = x^2 then f(x)dx = x^3/3
So (xval^3 – xmin^3)/3 = 0.3 then I just need to find xval

Is this right?
If so, I just need to find the integrated formula for all the probability distributions
e.g. if my distribution is x^2 then I don’t need to know this, I just need to know that my integrated formula is x^3/3

Thank you for taking the time to read this.
Posted
Updated 25-Jan-13 3:58am
v2
Sergey Alexandrovich Kryukov 25-Jan-13 10:04am
No, it's f(x)dx ≡ 1 :-)

So, before doing Monte-Carlo on some distribution, you need to get the distribution itself.

You further considerations are very non-mathematical, using some examples, etc., so they are hard to understand.
Usually, the libraries can give you only the uniform distribution (with some good or not so quality), sometimes, something else, like Gauss. Your task is to convert this distribution into some other distribution. For this, you need to write some function which converts uniform output into non-uniform, according to your requirements. Now, you have a simple mathematical problem (actually its simplicity depends on the distribution function you require :-). Now, it's your turn...

—SA
Jackie Lloyd 27-Jan-13 5:53am
Thank you for your reply. I don't know how to put in the mathematical symbols for integration so sorry that my examples seem non-mathematical. I think what you say is what I was trying to say, but I was not very clear. I have been a bit puzzled by your last few sentences, but I think I understand a bit now. Could you please confirm that 'you need to write some function which converts uniform output into non-uniform' means that the function that I need will use random numbers with a uniform distribution to generate input values for the monte carlo simulation so that more likely values are chosen more often (non-uniform). This function will be different for each distribution type. I think David has given me the Box-Muller algorithm (below) to achieve this for the Gaussian distribution, so I will research that more now. Thank you :)

Solution 1

In effect, what you have pointed out is correct.

The formal way to draw a value x from a probability distribution with a cumulative distribution function F(x) is by randomly choosing a value y between 0 and 1 in a uniform law of probability and then doing:

x = F(-1)(y)

where F(-1) is the inverse function of F, i.e:

x = F(-1)(y) <=> y = F(x)

But it is not always necessary to obtain F in order to draw a value from a probability distribution.

If your population comes from a normal (gaussian) law of probability, it is much easier to use the Box-Müller algorithm:

1.- Choose a value u from a U[0,1] (uniform law of probability in [0, 1]).
2.- Choose another value v from a U[0,1] (uniform law of probability in [0, 1]).
3.- Draw a value x from a N(mu, sigma) (gaussian with mean mu and standard deviation sigma) according to the following formula:

x = mu + sigma * square_root(2 * ln(1 / (1 - u))) * cos(2 * PI * v)

ln is neperian logarithm an PI is 3.141592...

I hope this helps

Jackie Lloyd 28-Jan-13 6:39am
Thank you so much for your help David, both for letting me know I was thinking along the right lines and helping me with the next step. I understand what the Box-Muller method is for and have managed to find C++ libraries which do this for a range of probabilty distributions, which is what I need. I'm really grateful to you :).

I have found some links to C++ libraries which generate random numbers for various PDs (pasted below incase anyone else needs some). If you happen to know of a good one that you have used, I would be very grateful if you forward me a link.
http://www.robertnz.net/nr02doc.htm
http://www.netlib.org/random/
http://www.johndcook.com/cpp_TR1_random.html