Common Generation Methods

Methods for generating pseudorandom numbers usually start with uniform random numbers, like the MATLAB® rand function produces. The methods described in this section detail how to produce random numbers from other distributions.

Direct Methods

Direct methods directly use the definition of the distribution.

For example, consider binomial random numbers. A binomial random number is the number of heads in N tosses of a coin with probability p of a heads on any single toss. If you generate N uniform random numbers on the interval (0,1) and count the number less than p, then the count is a binomial random number with parameters N and p.

This function is a simple implementation of a binomial RNG using the direct approach:

function X = directbinornd(N,p,m,n)

X = zeros(m,n); % Preallocate memory
for i = 1:m*n
    u = rand(N,1);
    X(i) = sum(u < p);
end

For example:

X = directbinornd(100,0.3,1e4,1);
hist(X,101)
set(get(gca,'Children'),'FaceColor',[.8 .8 1])

The Statistics Toolbox™ function binornd uses a modified direct method, based on the definition of a binomial random variable as the sum of Bernoulli random variables.

You can easily convert the previous method to a random number generator for the Poisson distribution with parameter λ. The Poisson distribution is the limiting case of the binomial distribution as N approaches infinity, p approaches zero, and Np is held fixed at λ. To generate Poisson random numbers, create a version of the previous generator that inputs λ rather than N and p, and internally sets N to some large number and p to λ/N.

The Statistics Toolbox function poissrnd actually uses two direct methods:

  • A waiting time method for small values of λ

  • A method due to Ahrens and Dieter for larger values of λ

Inversion Methods

Inversion methods are based on the observation that continuous cumulative distribution functions (cdfs) range uniformly over the interval (0,1). If u is a uniform random number on (0,1), then using X = F -1(U) generates a random number X from a continuous distribution with specified cdf F.

For example, the following code generates random numbers from a specific exponential distribution using the inverse cdf and the MATLAB uniform random number generator rand:

mu = 1;
X = expinv(rand(1e4,1),mu);

Compare the distribution of the generated random numbers to the pdf of the specified exponential by scaling the pdf to the area of the histogram used to display the distribution:

numbins = 50;
hist(X,numbins)
set(get(gca,'Children'),'FaceColor',[.8 .8 1])
hold on

[bincounts,binpositions] = hist(X,numbins);
binwidth = binpositions(2) - binpositions(1);
histarea = binwidth*sum(bincounts);

x = binpositions(1):0.001:binpositions(end);
y = exppdf(x,mu);
plot(x,histarea*y,'r','LineWidth',2)

Inversion methods also work for discrete distributions. To generate a random number X from a discrete distribution with probability mass vector P(X=xi) = pi where x0<x1< x2<..., generate a uniform random number u on (0,1) and then set X = xi if F(xi–1)<u<F(xi).

For example, the following function implements an inversion method for a discrete distribution with probability mass vector p:

function X = discreteinvrnd(p,m,n)

X = zeros(m,n); % Preallocate memory
for i = 1:m*n
    u = rand;
    I = find(u < cumsum(p));
    X(i) = min(I);
end

Use the function to generate random numbers from any discrete distribution:

p = [0.1 0.2 0.3 0.2 0.1 0.1]; % Probability mass vector
X = discreteinvrnd(p,1e4,1);
[n,x] = hist(X,length(p));
bar(1:length(p),n)
set(get(gca,'Children'),'FaceColor',[.8 .8 1]) 

Acceptance-Rejection Methods

The functional form of some distributions makes it difficult or time-consuming to generate random numbers using direct or inversion methods. Acceptance-rejection methods provide an alternative in these cases.

Acceptance-rejection methods begin with uniform random numbers, but require an additional random number generator. If your goal is to generate a random number from a continuous distribution with pdf f, acceptance-rejection methods first generate a random number from a continuous distribution with pdf g satisfying f(x) ≤ cg(x) for some c and all x.

A continuous acceptance-rejection RNG proceeds as follows:

  1. Chooses a density g.

  2. Finds a constant c such that f(x)/g(x)≤c for all x.

  3. Generates a uniform random number u.

  4. Generates a random number v from g.

  5. If cuf(v)/g (v), accepts and returns v.

  6. Otherwise, rejects v and goes to step 3.

For efficiency, a "cheap" method is necessary for generating random numbers from g, and the scalar c should be small. The expected number of iterations to produce a single random number is c.

The following function implements an acceptance-rejection method for generating random numbers from pdf f, given f, g, the RNG grnd for g, and the constant c:

function X = accrejrnd(f,g,grnd,c,m,n)
               
X = zeros(m,n); % Preallocate memory
for i = 1:m*n
    accept = false;
    while accept == false
        u = rand();
        v = grnd();
        if c*u <= f(v)/g(v)
           X(i) = v;
           accept = true;
        end
    end
end

For example, the function f(x) = xe–x2/2 satisfies the conditions for a pdf on [0,∞) (nonnegative and integrates to 1). The exponential pdf with mean 1, f(x) = e–x, dominates g for c greater than about 2.2. Thus, you can use rand and exprnd to generate random numbers from f:

f = @(x)x.*exp(-(x.^2)/2);
g = @(x)exp(-x);
grnd = @()exprnd(1);
X = accrejrnd(f,g,grnd,2.2,1e4,1);

The pdf f is actually a Rayleigh distribution with shape parameter 1. This example compares the distribution of random numbers generated by the acceptance-rejection method with those generated by raylrnd:

Y = raylrnd(1,1e4,1); 
hist([X Y])
h = get(gca,'Children');
set(h(1),'FaceColor',[.8 .8 1])
legend('A-R RNG','Rayleigh RNG')

The Statistics Toolbox function raylrnd uses a transformation method, expressing a Rayleigh random variable in terms of a chi-square random variable, which you compute using randn.

Acceptance-rejection methods also work for discrete distributions. In this case, the goal is to generate random numbers from a distribution with probability mass Pp(X = i) = pi, assuming that you have a method for generating random numbers from a distribution with probability mass Pq(X = i) = qi. The RNG proceeds as follows:

  1. Chooses a density Pq.

  2. Finds a constant c such that pi/qic for all i .

  3. Generates a uniform random number u.

  4. Generates a random number v from Pq.

  5. If cupv/qv, accepts and returns v.

  6. Otherwise, rejects v and goes to step 3.

Was this topic helpful?