Code covered by the BSD License  

Highlights from
Random Sample from Discrete PDF

5.0

5.0 | 1 rating Rate this file 39 Downloads (last 30 days) File Size: 2.07 KB File ID: #37698
image thumbnail

Random Sample from Discrete PDF

by

 

Inverse transform sampling to generate random sample from pdf given by domain x and p(x).

| Watch this File

File Information
Description

Generate random samples given a pdf. The x and p(x) are specified in the
arguments. See inverse transform sampling, gaussdis, gammadis.

Examples:
X = pdfrnd(0:.1:20, gammadis(3,2, 0:.1:20), 100000);
figure, hist(X, 100);
%X will be 100000 samples from gamma distribution with a/k = 3, b/Theta =
%2.

X = pdfrnd(-10:.1:10, gaussdis(2, 3, -10:.1:10), 100000);
figure, hist(X, 100);
%X will be 100000 samples from a gaussian distribution with mean 2, var 3.
%Should be roughly equivalent to X = sqrt(3)*randn(100000, 1) + 2;

Acknowledgements

Random Numbers From A Discrete Distribution inspired this file.

Required Products MATLAB
MATLAB release MATLAB 7.10 (R2010a)
Tags for This File   Please login to tag files.
Please login to add a comment or rating.
Comments and Ratings (3)
15 Dec 2013 Joshua Stough

Response to Chavas:
Thank you so much for your comment. The problem as I see it is that with the assumption of 0, it is possible to generate samples from outside of the provided domain x. Even for the examples I provide, cdf(1) > 0 and so invalid samples can be generated. However, your proposed solution could also generate such samples. Perhaps to account for binning in x, x(1) - dx/2 should be possible as a minimum, and x(end) + dx/2 as a max (assuming equal sampling of the domain). The other option is to assume the user provides precisely the x(1) and x(end) they want to sample within. I will post new code soon. Thank you again.

13 Dec 2013 Dan Chavas

Whoops, there is a problem: if you input a pdf (px) such that px(1)>0, then cdf(1)>0. However it is possible for rnd < cdf(1) ; in such a case, interp1 returns NaN, which you've forced to be 0. This is wrong.

Consider adding this code to ensure the cdf begins at 0:

%%Account for non-zero probability at start
cdf = cumsum(px);
if(cdf(1)>0)
cdf = [0 cdf];
dx = x(2) - x(1);
x0 = x(1) - dx;
x = [x0 x];
end

rnd = rand(sampleSize, 1);

X = interp1(cdf, x, rnd, 'linear');

12 Dec 2013 Dan Chavas

Just wrote this code myself to apply specifically to a power-law distribution. Then I found this file for arbitrary pdfs, tried it out and they match! Simple but very useful, this should be included as a standard MATLAB function.

Contact us