what called this method of generation random sample

 Hi all

is the code below represent inverse cdf method?

y is pdf of any distribution

cdf_y = cumsum(y);

sum_y = sum(y);

for j = 1:N

randx = sum_y*rand();

i = 1;

while cdf_y(i) < randx

i = i + 1;

end

f(j) = x(i); end

please give me explain

[José-Luis]

Please edit your post.

http://www.mathworks.com/matlabcentral/answers/13205-tutorial-how-to-format-your-question-with-markup

Answer 1

[Roger Stafford]

In the code you show, Mutah, strictly speaking the y vector is not a pdf - a probability distribution function. It is a pmf - a probability mass function. Actually the author of the code has allowed y to be merely proportional to these probabilities, so y is proportional to a pmf. That is because it is dealing with a discrete set of values, not a continuous distribution. For each i, y(i) must be proportional to the probability of a random variable, f, being equal to the corresponding x(i).

With that understanding, the code is correct and the N-length array, f, will have a statistical distribution in accordance with the y values. The reason the author finds the cumsum(y) and sum(y) is that cumsum(y) is the cumulative sum of the y values and therefore proportional to the cumulative probability that is required of the result, and sum(y) would be the proportionality constant. If on the i-th trip through the for-loop, we have 0 < randx <= cdf_y(1), the while loop will not execute at all and you will get i = 1 which will insert x(1) into f. If cdf_y(1) < randx <= cdf_y(2), the while loop will make just one trip and i will emerge equal to 2, which inserts x(2) into f. This continues to be true all the way to the case cdf_y(end-1) < randx < cdf_y(end) which would then put x(end) into f. Therefore this has the effect of finding the inverse function of cdf_y, that is for each cumulative probability value, randx, it provides the appropriate x(i) value.

However, it must be said that this code is not very efficient. It does not take advantage of the fact that the cdf_y values are an ascending sequence of values. It should not be necessary to possibly have to run through the entire cdf_y sequence with the while loop seeking the point where each randx fits into the correct interval of cdf_y values. The following code uses the matlab function 'histc' which does take advantage of the ordering in cdf_f and accordingly should be much faster code for large x and y arrays.

 cdf_y = cumsum(y);
 cdf_y = cdf_y/cdf_y(end);
 [~,ix] = histc(rand(N,1),[0;cdf_y]);
 f = x(ix);

I have assumed here that y and x are column vectors. To understand why this works, you should make a careful study of the 'histc' function at:

 http://www.mathworks.com/help/matlab/ref/histc.html

[mutah]

thank you. very nice explain

[mutah]

Do you mean we replace while loop by

cdf_y = cumsum(y);

 cdf_y = cdf_y/cdf_y(end);

 [~,ix] = histc(rand(N,1),[0;cdf_y]);

 f = x(ix);

or what ?

[Roger Stafford]

Yes, that is what I mean. Matlab's 'histc' function does the searching for you instead of using a 'while' loop and it generates N samples instead of using a 'for' loop. Try it out!

[mutah]

 Error using histc
Edge vector must be monotonically non-decreasing.

[Roger Stafford]

"y is pdf of any distribution" The only way I can see that your edge vector failed to be monotonically non-decreasing is that one of your 'y' values was negative. However, the pdf of distributions should always be non-negative.

You should be able to check that for yourself. Find the point where cdf_y decreases and check the corresponding y value.

Answer 2

[Image Analyst]

Not sure exactly what you mean and how specific your question is to that exact code, but in general the process of generating a bunch of random numbers and running an "experiment" N times (in a loop like you did) is called a "Monte Carlo Simulation". For example, check out the attached Monty Hall Problem that I coded up as a Monte Carlo simulation.

[mutah]

 there are many methods for generation random samples 
like :1- inverse cdf method
2- acceptance rejection method

the code above which method represent

[Image Analyst]

That would be inverse transform sampling, http://en.wikipedia.org/wiki/Inverse_transform_sampling. See my attached demo which works with a Rayleigh distribution.

[mutah]

Do you mean the code above represent inverse cdf method?

if yes please give explain

[Image Analyst]

Well, sort of. I think that may be the intent of the author but they really messed up the code, and that's not even talking about the formatting. I did attach an example in my prior comment.

[mutah]

please see the code above

and give me your opinion about code

is this code represent inverse cdf ?

[Image Analyst]

It's a Monte Carlo simulation but the problem is that x is not defined. A skilled MATLAB programmer would also use find() rather than while.

[mutah]

x=0:0.001:10

y is pdf(x)

why the author use cumsum(y) and sum(y)

[Image Analyst]

cumsum gives the cdf of y. For that x, the pdf if flat since it's a uniform distribution so the cdf will be a straight ramp. They use sum(y) to get the max value that the cdf could be. They just as well could have used cdf(end), because the cdf and pdf are not normalized and the sum of all y values is the number of x elements you have, which is 10,001.

[mutah]

what about this

N=10000

for j = 1:N

randx = sum_y*rand();

i = 1;

while cdf_y(i) < randx

i = i + 1;

end

f(j) = x(i);

end

[Image Analyst]

Well, like I said, the while (and line before and two lines after) could be replaced:

i = find(cdf_y > randx, 1, 'first')