Random Number Generation
Random number generation
These methods only generate pseudo-random numbers for trivial reasons. Two desirable property for random numbers is:
- each \(U_i\) is uniformly distributed between \(0\) and \(1\);
- the \(U_i\)’s are mutually independent.
Note that property 1 is a convenient but arbitrary normalization. In layman terms, the value of \(U_i\) should not be predictable from \(U_1, \cdots, U_{i-1}\).
Linear congruential generator
A linear congruential generator is a recurrence of the form
\[x_{i+1} = ax_i \mod m\]
\[u_{i+1} = x_{i+1}/m \]
The value of \(x_0\) is called seed and \(a\), \(m\) are integers. L.C.G. is the most widely used generators.
One can see that the pseudo-random numbers generated by L.C.G. generated has a period of at most \(m\), i.e., repeats. For some choice of \(a\) and \(m\), they attain full-period or maximum period and in some other cases, they attain a period smaller than \(m-1\). In general full-period is a desired criteria.
The general form of L.C.G takes the following form \[x_{i+1} = (ax_i+c) \mod m\] \[u_{i+1} = x_{i+1}/m \]
There are various conditions for which full period can be generated and can be found in the book Art of computer programming, semi numeric algorithms by Donald Knuth. Marsaglia was able to prove that the plot of the \(k\)-tuple \((u_i,\cdots, u_{i+k})\) will be hyper-planes in \(k\)-dimension and was also able to arrive at a closed form for maximum number of hyper-planes; this is known as lattice structure of linear congruential generators. The algorithm for implementing linear congruential generators is easy and therefore skipped.
Inverse transform method
Suppose we want to generate a sample from a cumulative distribution function \(F\), i.e., we want to generate a random variable \(X\) with the property that \(P(X\le x) = F(x)\) for all \(X\). The inverse transform method set \[ X = F^{-1}(U) \quad U \sim \text{Unif}[0,1],\] where \(F^{-1}\) is the inverse of \(F\) and \(\text{Unif}[0,1]\) denotes the uniform distribution on \([0,1]\).
The algorithm for this method is pretty straightforward and hence omitted.
Inverse transform method for discrete distributions.
In case of discrete random variables, evaluation of \(F^{-1}\) reduces to a table lookup. Consider a discrete random variable whose possible values \(c_1<\cdots c_n\). Let \(p_i\) be the probability attached to \(c_i, i= 1,\cdots,n\), and set \(q_0 = 0\), \[q_i = \sum\limits_{j=1}^{i} p_j, \quad i,\cdots,n.\] These are the cumulative probabilities associated with \(q_i = F(c_i), i = 1,\cdots, n\). To sample from this distribution,
- generate a uniform \(U\);
- find \(K\in\{1,\cdots,n\}\) such that \(q_{K-1} < U \le q_k\);
- set \(X = c_K\).
The second step may be implemented using binary search as \(q_k\) is an increasing sequence.
Acceptance Rejection method
Acceptance rejection method is a widely applicable method. The rejection mechanism is designed so that the accepted samples are indeed according to the target distribution.
Suppose that we wish to generate samples from a density \(f\) defined on some set \(X\), let \(g\) be the density on \(X\) from which we know how to generate samples and with the property that \[f(x) \le c g(x), \text{ for all } x\in X\] for some constant \(c\). In acceptance rejection method we generate a sample \(X\) from \(G\) and accept the sample with probability \(f(X)/cg(X)\); this can be implemented by sampling \(U\) uniformly over \((0,1)\) and accepting \(X\) if \(U\le F(X)/cg(X)\). If \(X\) is rejected, a new sample is sampled from \(g\) and the acceptance test is applied again.
Algorithm
{% gist b56a1c8f5d961e98d77e %}
Some typical examples are Beta distribution, normal from double exponential,
Normal random variables and vectors
Suppose we have a sample from \(N(0,1)\), i.e., normal distribution from parameters \(\mu = 0\) and \(\sigma = 1\), then from this sample we can generate \(N(\mu, \sigma)\) by observing that \(\mu + \sigma Z \sim N(\mu, \sigma^2)\).
Box-Muller method
A simple method to generate sample from normal distribution is Box-Muller. The algorithm for this method is given below
{% gist 572520bf1f070295fa85 %}
Marsaglia-Bray method
Marsaglia and Bray developed a modification of the Box-Muller method that reduces computing time by avoiding evaluation of sine and cosine functions. (I’m told that calculation of \(log\), \(\sqrt{}\), requires lesser amount of time than \(\sin\) and \(\cos\))
Algorithm
{% gist a1718c75604807350c57 %}
References
- Monte Carlo methods in financial engineering, Paul Glasserman.