Direct methods make direct use of the definition of the distribution.
For example, consider binomial random numbers. You can think of a binomial random number as the number of heads in N tosses of a coin with probability p of a heads on any 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.
The following function is a simple implementation of a binomial RNG using this 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)
Direct Methods (another binomial example)
You can think of binomial random numbers as the number of heads in n tosses of a coin with probability p of a heads on any toss. If you generate n uniform random numbers and count the number that are greater than p, the result is binomial with parameters n and p.
determinations = 10000;
n = 100;
p = 1 - 0.3;
result = zeros(determinations, 1);
for i = 1:determinations
count = 0;
for j = 1:n
if(rand > p)
count = count + 1;
end
end
result(i,1) = count;
end
hist(result, 101);
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 a random number X from a continuous distribution with specified cdf F can be obtained using X = F-1(U).
Exponential CDF with parameter lambda = F(x, lambda) = p = 1 - exp(-lambda * x)
Inverse exponential CDF with parameter lambda = F-1(p, lambda) = x = -ln(1 - p)/lambda
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:
determinations = 10000;
lambda = 1;
result = zeros(determinations, 1);
for i = 1:determinations
result(i,1) = -log(1 - rand) / lambda;
end
numbins = 101;
hist(result,numbins)
hold on
[bincounts,binpositions] = hist(result,numbins);
binwidth = binpositions(2) - binpositions(1);
histarea = binwidth*sum(bincounts);
x = binpositions(1):0.001:binpositions(end);
y = (1 / lambda) * exp(-x / lambda);
plot(x,histarea*y,'r','LineWidth',2)
Continuous Distribution Inversion Methods
Inversion methods can be adapted to discrete distributions. Suppose you want a random number X from a discrete distribution with a probability mass vector P(X =
xi) = pi, where x0 < x1 < x2 < ... . You could 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:
determinations = 10000;
p = [0.1 0.2 0.3 0.2 0.1 0.1]; % Probability mass vector;
sum_p = zeros(length(p),1);
sum_p(1) = p(1);
for i = 2:length(p)
sum_p(i) = p(i) + sum_p(i - 1);
end
result = zeros(determinations, 1);
for i = 1:determinations
x = rand;
for j = 1:length(p)
if x < sum_p(j);
result(i,1) = j;
break;
end
end
end
hist(result, length(p));
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 can provide a good solution in these cases.
Acceptance-rejection methods also begin with uniform random numbers, but they require the availability of an additional random number generator. If the 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. Choose a density g.
2. Find a constant c such that f (x) / g (x) ≤ c for all x.
3. Generate a uniform random number u.
4. Generate a random number v from g.
5. If c*u ≤ f (v) / g (v) , accept and return v.
6. Otherwise, reject v and go to 3.
For efficiency, you need a cheap method for generating random numbers from g, and the scalar c should be small. The expected number of iterations to produce a random number is c.
For example, the function f (x) = xe�x^2/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
to generate random numbers from f :
determinations = 10000;
result = zeros(determinations, 1);
c = 2.2;
% f = x.*exp(-(x.^2)/2);
% g = exp(-x);
X = zeros(determinations, 1); % Preallocate memory
for i = 1:determinations
accept = false;
while accept == false
u = rand();
v = -log(1 - rand) / lambda; % exponential rng using the inversion method
if c*u <= v.*exp(-(v.^2)/2) / exp(-v)
result(i) = v;
accept = true;
end
end
end
hist(result, 101);
* portions from here.