OnExtraBit

Generate random number based on a distribution in Octave/Matlab

Trung Nguyen posted on Thu, 31 May 2012 13:26:04 CEST+0200

Updated on Sat, 02 Jun 2012 22:07:16 CEST+0200

Updated: Johannes' gauss_rand version that supports mean and standard deviation.

In some machine algorithms I need to randomly generate weights for a model (hypothesis), normally, using uniform distribution. But given a distribution, generating random numbers that follow the given distribution is also very interesting.

I've just found a method call Inverse transform sampling on Wikipedia. It turns out that it's a quite effective and easy way to implement a random generator.

The idea is: Given a distribution with cumulative distribution function $F(x)$ (CDF), we can generate a random number $x$ from the given distribution by solving $F(x) = u$ with $u$ is a random number generated from uniform distribution.

Let's try to write a Gaussian (or normal distribution) based random generator in Octave/Matlab.

We have Gaussian CDF.

$$F(x; \mu, \sigma^2) = \frac{1}{2}[1 + \mbox{erf}(\frac{x - \mu}{\sqrt{2\sigma^2}})]$$

$\mbox{erf}$ is error function.

Standard Gaussian distribution has $\mu = 0$ and $\sigma^2 = 1$, so we have standard CDF.

$$F(x) = \frac{1}{2}[1 + \mbox{erf}(\frac{x}{\sqrt{2}})]$$

We can find $x$ (our random number) by solving below equation.

$$u = \frac{1}{2}[1 + \mbox{erf}(\frac{x}{\sqrt{2}})]$$

$u$ is a random number generated from uniform distribution. In Octave/Matlab, we can get $u$ by calling function $rand()$. $x$ will be calculated as.

$$x = \sqrt{2}\mbox{erfinv}(2u - 1)$$

$\mbox{erfinv}$ is the inversion of $\mbox{erf}$ function.

The Octave/Matlab code will be.

function r = gauss_rand(n, m)
    % Generate random number based on inverse transform sampling method
    %   standard gauss distribution mean = 0, variance = 1
    % n: number of rows
    % m: number of columns

    r = sqrt(2) * erfinv(2 * rand(n, m) - 1);
end

Test to see how the histogram of our Gaussian based random generator works.

hist(gauss_rand(1000, 1), 50)

Here's the plot.

Standard Gaussian distribution

Thank Johannes for his gauss_rand version that supports mean and standard deviation :)

function r = gauss_rand(n, m, mu, sigma)
    % Generate random number based on inverse transform sampling method
    % n: number of rows
    % m: number of columns
    % mu: mean of distribution
    % sigma: standard deviation of distribution

    if nargin < 4
        mu = 0;
        sigma = 1;
    end
    r = sqrt(2) * sigma * erfinv(2 * rand(n, m) - 1) + mu;
end

Tags: random, gaussian, octave, matlab

comments powered by Disqus

Copyright © 2012 by Trung Nguyen.   Written in Python & powered by Tornado