Code covered by the BSD License

### Highlights from Random Sample from Discrete PDF

5.0
5.0 | 2 ratings Rate this file 5 Downloads (last 30 days) File Size: 2.07 KB File ID: #37698 Version: 1.0

# Random Sample from Discrete PDF

### Joshua Stough (view profile)

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

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)
28 Aug 2014 Martin Eisemann

### Martin Eisemann (view profile)

Maybe I am mistaken here, so please correct me if I am wrong, but wouldn't this help?

cdf = cumsum(px);
cdf = cdf - min(cdf);
cdf = cdf/max(cdf);

rnd = rand(sampleSize, 1);
X = interp1(cdf, x, rnd, 'linear', 0);

15 Dec 2013 Joshua Stough

### Joshua Stough (view profile)

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.

Comment only
13 Dec 2013 Dan Chavas

### Dan Chavas (view profile)

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');

Comment only
12 Dec 2013 Dan Chavas

### Dan Chavas (view profile)

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.