View License

Download apps, toolboxes, and other File Exchange content using Add-On Explorer in MATLAB.

» Watch video

Highlights from
Sampling from a discrete distribution

4.8 | 8 ratings Rate this file 41 Downloads (last 30 days) File Size: 2.25 KB File ID: #21912 Version: 1.1

Sampling from a discrete distribution


Dahua Lin (view profile)


27 Oct 2008 (Updated )

The function is to draw samples from an arbitrary discrete distribution.

| Watch this File

File Information

There are a lot of cases that you might need to sample from a discrete distribution in Monte Carlo simulations. Here are some typical examples that this function may help

(1) you want to sample from a discrete distribution over finitely many categories(labels). However, the pmf is not easy to directly sample from. Then you can just compute p(k) for each category, and then use this function to do the sampling.

(2) For a general (possibly non-parametric) continuous distribution defined in a 1D/2D space, you may discretized the sample space into small regions, and compute the probability mass for these regions, and then use this function to do the sampling.

(3) In many models, the probability distribution is expressed as a weighted sum of several "modes". To sample from this kind of distributions, you may need to first choose which mode to generate the sample according to their weights, then this function may fit in.

In sum, this function can be used directly in simple sampling, and may also be used as a building block of complicated simulation procedure.

The use of this function is easy:

x = discretesample(p, n).

You just input the probability mass, and tell the function how many sample you want to sample, then it returns the samples in form of a 1 x n vector.


This file inspired Lynx Matlab Toolbox and Fast Sampling From A Discrete Distribution.

MATLAB release MATLAB 7.6 (R2008a)
Tags for This File   Please login to tag files.
Please login to add a comment or rating.
Comments and Ratings (14)
21 Feb 2017 Rani Barbara


I have a joint probability distribution of 4 variables. and I want to sample from it using this function. when I pass my 4D array to the function it seems it collapses the 4D array into a one long vector. I get the indexes of the sampled but I am unable to figure out how to take the indexes back to the 4D form to know which combination is sampled.
alternatively I thought to sample the marginals and each 4 samples of the marginals will be one single sample from the joint. is the latter equivalent? or maybe you have some other advise.
Thank you

Comment only
21 Oct 2016 Narayanan Rengaswamy

10 Feb 2016 FM

FM (view profile)

To further Mo Chen's approach: histogram(discretize(rand(1,n),[0;cumsum(p(:))/sum(p)]))

Comment only
02 Jan 2015 Kevin

Kevin (view profile)

Excellent product. Thank you.

08 Nov 2013 Justin

Justin (view profile)

Very helpful! Thx!

10 Jul 2013 Christos

Very useful, saves a lot of time when it comes in large scale simulation problems.

21 Jun 2012 Woojae Kim

Mo Chen, that is brilliant! That should be the fastest way to get the job done in MATLAB. Thanks!

13 Jan 2012 Thomas Clark

Very useful, a good utility function which saves some though when implementing MC models.

22 Feb 2011 Reza Fathzadeh

05 Nov 2009 Mo Chen

Mo Chen (view profile)

your function seems to complicate the problem a little bit.
The following line is enough to do the jod
[~,x] = histc(rand(1,n),[0;cumsum(p(:))/sum(p)]);

Comment only
21 Jun 2009 Carlos Baiz

29 Oct 2008 Dahua Lin

Dahua Lin (view profile)


I have checked with randp. Though they seem offering similar functionalities, however, the efficiency is drastically different, especially in very large scale monte carlo simulation, say you need to draw thousands or millions of samples from a distribution over thousands or even millions of states, which is not unusual in real engineering applications.

With randp, it would be incur obvious latency when you want to draw thousands of samples from thousands of states, or even run out of memory (thus resulting an empty matrix), as the algorithm implemented by randp is of complexity O(k n), where k is the number of states in the sample space. However, even million-state-level sampling can be accomplished by this function within milliseconds, as its complexity is only O(n logk).

Hope my explanation can clarify the differences between these two files.

Comment only
28 Oct 2008 Jos (10584)

Jos (10584) (view profile)

if you'd searched before submitting you would have found RANDP on the FEX. What does this function add?

Comment only
28 Oct 2008 Dahua Lin

Dahua Lin (view profile)

There's a bug in the initial submission, which has been fixed, and the new version has been uploaded, which may come out in a day.

Comment only
28 Oct 2008 1.1

Fixing a bug that may lead to incorrect sampling.

Contact us