Generate number from a probability distribution

4 vues (au cours des 30 derniers jours)
Bernoulli Lizard
Bernoulli Lizard le 26 Juin 2013
I am trying to generate numbers from a probability distribution using the inverse cdf method. However, I can not calculate the CDF analytically. I can instead calculate the discrete PDF and use cumsum to get the discrete CDF. After generating random numbers from the uniform distribution from 0 to 1, how can I use the discrete CDF to generate numbers from my PDF? How would the random numbers on 0 to 1 corrospond to the values from my CDF?
thanks

Réponses (1)

Roger Stafford
Roger Stafford le 27 Juin 2013
(I presume that since you took a simple cumsum of your pdf values, these were equally spaced with respect to the random variable's values. Otherwise it should be a weighted cumsum, that is, each density value should be weighted by the length of interval it represents.)
Regardless of how you have obtained the discrete cdf values, you can use 'histc' to take the place of an inverse of the cdf function. We assume that you have two column vectors, c and F of the same length, n, where the c values are the successive random variable values and F the corresponding cumulative probability values - that is, that the probability that the random variable is less than or equal to c(i) is equal to F(i). Also we assume that F(1) = 0 and F(n) = 1 and that F is ascending with no two successive values equal. Then do the following:
u = rand(N,1);
[~,b] = histc(u,F);
x = (u-F(b))./(F(b+1)-F(b)).*(c(b+1)-c(b)) + c(b);
Then x will be a vector of length N which approximates the required random variable distribution in accordance with the F cdf vector. The third line gives a linear interpolation between c(b) and c(b+1) since the corresponding u is between F(b) and F(b+1). The value of index b should never be 0 or n if you have F(1) = 0 and F(n) = 1 because 'rand' should never be less than 0 or greater than or equal to 1.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by