(New page: == '''1. Introduction''' == In this slecture, we will explain the principle of how to generate Gaussian random samples. Even though there are more general methods to generate random samp...)
 
Line 1: Line 1:
 
 
== '''1. Introduction''' ==
 
== '''1. Introduction''' ==
  
Line 7: Line 6:
 
== '''2. Inverse transform sampling''' ==
 
== '''2. Inverse transform sampling''' ==
  
Let U be a random variable which is uniformly distributed on the interval [0, 1]. And let F be a continuous CDF(cumulative distribution function). Then, inverse CDF is defined by
+
Let U be a random variable which is uniformly distributed on the interval [0, 1]. And let F be a continuous CDF(cumulative distribution function) of a random variable, X which we want to generate. Then, inverse CDF is defined by
 +
 
 +
<center><math>
 +
F^{-1}(u)=inf\{ x|F(x)\leq u, u\in [0, 1] \}
 +
</math></center>
 +
 
 +
Using this inverse CDF, we can generate random variable X as following
 +
 
 +
<center><math>
 +
X <- F^{-1}(U)\quad
 +
</math></center>
 +
 
 +
Now, we will apply this to generate a random variable X which has the exponential distribution
 +
 
 +
<center><math>
 +
f(x) = \lambda exp(-\lambda x) \quad where \quad x \geq 0
 +
</math></center>
 +
 
 +
Let's compute CDF of X
 +
 
 +
<center><math>
 +
F(x) = \int_{-\infty}^x \lambda exp(-\lambda x') dx' = \int_0^x \lambda exp(-\lambda x') dx' = [-exp(-\lambda x')]_0^x = 1-exp(-\lambda x) \leq u
 +
</math></center>
 +
 
 +
Then, we can compute x as
 +
<center><math>
 +
x = -\lambda^{-1}ln(1-u) \quad
 +
</math></center>
 +
 
 +
In other words, we can get an exponentially distributed random variable X using an uniformly distributed random variable U
 +
 
 +
<center><math>
 +
X = -\lambda^{-1}ln(1-U) \quad
 +
</math></center>
 +
 
 +
This method can be applied to generate any random variable if the CDF of it is easily calculated. If it is not, we cannot use this method. Instead, we can use other methods such as rejection sampling, importance sampling, Hastings-Metropolis sampling, etc. However, we do not cover those sampling methods here.
 +
 
 +
 
 +
== 3. Box Muller Transform ==
 +
 
 +
Before explaing Box Muller transform, we will apply Inverse transform sampling to generate a normal distributed random variable.
 +
 
 +
<center><math>
 +
F(x) = \int_{-\infty}^x \frac{1}{\sqrt{2\pi}} exp(-\frac{1}{2} x'^2) dx'
 +
</math></center>
 +
 
 +
This integration is not straightforward to make an algebraic expression in terms of elementary functions such as exponentials and logs. Therefore, we use Box Muller transform method to generate a normal distributed random variable.
 +
 
 +
Let (X, Y) be a pair of independent normal distributed random variable. Then, we can express F(x,y) as follows
 +
 
 +
<center><math>
 +
F(x,y) = \int_{-\infty}^y \int_{-\infty}^x \{ \frac{1}{\sqrt{2\pi}} exp(-\frac{1}{2} x'^2) \} \{ \frac{1}{\sqrt{2\pi}} exp(-\frac{1}{2} y'^2)\} dx'dy' = \int_{-\infty}^y \int_{-\infty}^x \frac{1}{2\pi} exp(-\frac{1}{2} (x'^2 + y'^2) ) dx'dy'
 +
</math></center>
 +
 
 +
Now, define random variable R and θ as follows
 +
<center><math>
 +
R = \sqrt{X^2 + Y^2}, \quad  \theta = tan(\frac{Y}{X}) \quad where \quad \theta \in [0, 2\pi]
 +
</math></center>
 +
 
 +
X and Y can be expressed using R and θ as follows
 +
 
 +
<center><math>
 +
X = Rcos(\theta) ,\quad Y = Rsin(\theta)
 +
</math></center>
 +
 
 +
Then, we can express F(x,y) using r and θ (dxdy = rdrdθ)
 +
 
 +
<center><math>
 +
F(x,y) = F(r,\theta) = \int_{-\infty}^\theta \int_{-\infty}^r \frac{1}{2\pi} exp(-\frac{1}{2} r'^2 ) r' dr'd\theta'  = \int_{0}^\theta \int_{0}^r \frac{1}{2\pi} exp(-\frac{1}{2} r'^2 ) r' dr'd\theta' =  \frac{1}{2\pi} \int_{-\infty}^\theta [-exp(-\frac{1}{2} r'^2 )]_0^r d\theta'
 +
</math></center>
 +
<center><math>
 +
= \frac{\theta}{2\pi}  \{ 1-exp(-\frac{1}{2} r^2 ) \} = F(r)F(\theta)
 +
</math></center>
 +
 
 +
Let <math>U_1</math> and <math>U_2</math> be uniformly distributed random variables such that <math> U_1,U_2 \in [0, 1] </math> . Then we can say as follows
 +
 
 +
<center><math>
 +
F(\theta) \leq u_1, \quad F(r) \leq u_2
 +
</math></center>
 +
 
 +
Since <math> U_2 \in [0, 1] </math>, we can say <math> 1-U_2 \in [0, 1]  </math>. Then
 +
 
 +
<center><math>
 +
F(\theta) \leq u_1, \quad F(r) \leq 1-u_2
 +
</math></center>
 +
 
 +
Therefore, we can get <math>u_1</math> and <math>u_2</math> as follows
 +
 
 +
<center><math>
 +
u_1 = \frac{\theta}{2\pi}, \quad \theta = 2\pi u_1
 +
</math></center>
 +
<center><math>
 +
1-u_2 = 1-exp(-\frac{1}{2} r^2 ), \quad u_2 = exp(-\frac{1}{2} r^2), \quad r = \sqrt{-2ln(u_2)}
 +
</math></center>
 +
 
 +
Now, we can generate X and Y using <math>U1</math> and <math>U2</math>.
 +
<center><math>
 +
X = \sqrt{-2ln(U_2)}cos( 2\pi U_1 ),
 +
\quad u_2 = exp(-\frac{1}{2} r^2), \quad r = \sqrt{-2ln(u_2)}
 +
</math></center>

Revision as of 00:45, 1 May 2014

1. Introduction

In this slecture, we will explain the principle of how to generate Gaussian random samples. Even though there are more general methods to generate random samples which have any distribution, we will focus on the simple method such as Box Muller transform to generate Gaussian random samples in this slecture. To do this, we will introduce inverse transform sampling which is a bas method for pseudo random number sampling first. Then, we will explain Gaussian random sample generation method based on Box Muller transform. Finally, we will introduce Marsaglia polar method which improves Box Muller transform.


2. Inverse transform sampling

Let U be a random variable which is uniformly distributed on the interval [0, 1]. And let F be a continuous CDF(cumulative distribution function) of a random variable, X which we want to generate. Then, inverse CDF is defined by

$ F^{-1}(u)=inf\{ x|F(x)\leq u, u\in [0, 1] \} $

Using this inverse CDF, we can generate random variable X as following

$ X <- F^{-1}(U)\quad $

Now, we will apply this to generate a random variable X which has the exponential distribution

$ f(x) = \lambda exp(-\lambda x) \quad where \quad x \geq 0 $

Let's compute CDF of X

$ F(x) = \int_{-\infty}^x \lambda exp(-\lambda x') dx' = \int_0^x \lambda exp(-\lambda x') dx' = [-exp(-\lambda x')]_0^x = 1-exp(-\lambda x) \leq u $

Then, we can compute x as

$ x = -\lambda^{-1}ln(1-u) \quad $

In other words, we can get an exponentially distributed random variable X using an uniformly distributed random variable U

$ X = -\lambda^{-1}ln(1-U) \quad $

This method can be applied to generate any random variable if the CDF of it is easily calculated. If it is not, we cannot use this method. Instead, we can use other methods such as rejection sampling, importance sampling, Hastings-Metropolis sampling, etc. However, we do not cover those sampling methods here.


3. Box Muller Transform

Before explaing Box Muller transform, we will apply Inverse transform sampling to generate a normal distributed random variable.

$ F(x) = \int_{-\infty}^x \frac{1}{\sqrt{2\pi}} exp(-\frac{1}{2} x'^2) dx' $

This integration is not straightforward to make an algebraic expression in terms of elementary functions such as exponentials and logs. Therefore, we use Box Muller transform method to generate a normal distributed random variable.

Let (X, Y) be a pair of independent normal distributed random variable. Then, we can express F(x,y) as follows

$ F(x,y) = \int_{-\infty}^y \int_{-\infty}^x \{ \frac{1}{\sqrt{2\pi}} exp(-\frac{1}{2} x'^2) \} \{ \frac{1}{\sqrt{2\pi}} exp(-\frac{1}{2} y'^2)\} dx'dy' = \int_{-\infty}^y \int_{-\infty}^x \frac{1}{2\pi} exp(-\frac{1}{2} (x'^2 + y'^2) ) dx'dy' $

Now, define random variable R and θ as follows

$ R = \sqrt{X^2 + Y^2}, \quad \theta = tan(\frac{Y}{X}) \quad where \quad \theta \in [0, 2\pi] $

X and Y can be expressed using R and θ as follows

$ X = Rcos(\theta) ,\quad Y = Rsin(\theta) $

Then, we can express F(x,y) using r and θ (dxdy = rdrdθ)

$ F(x,y) = F(r,\theta) = \int_{-\infty}^\theta \int_{-\infty}^r \frac{1}{2\pi} exp(-\frac{1}{2} r'^2 ) r' dr'd\theta' = \int_{0}^\theta \int_{0}^r \frac{1}{2\pi} exp(-\frac{1}{2} r'^2 ) r' dr'd\theta' = \frac{1}{2\pi} \int_{-\infty}^\theta [-exp(-\frac{1}{2} r'^2 )]_0^r d\theta' $
$ = \frac{\theta}{2\pi} \{ 1-exp(-\frac{1}{2} r^2 ) \} = F(r)F(\theta) $

Let $ U_1 $ and $ U_2 $ be uniformly distributed random variables such that $ U_1,U_2 \in [0, 1] $ . Then we can say as follows

$ F(\theta) \leq u_1, \quad F(r) \leq u_2 $

Since $ U_2 \in [0, 1] $, we can say $ 1-U_2 \in [0, 1] $. Then

$ F(\theta) \leq u_1, \quad F(r) \leq 1-u_2 $

Therefore, we can get $ u_1 $ and $ u_2 $ as follows

$ u_1 = \frac{\theta}{2\pi}, \quad \theta = 2\pi u_1 $
$ 1-u_2 = 1-exp(-\frac{1}{2} r^2 ), \quad u_2 = exp(-\frac{1}{2} r^2), \quad r = \sqrt{-2ln(u_2)} $

Now, we can generate X and Y using $ U1 $ and $ U2 $.

$ X = \sqrt{-2ln(U_2)}cos( 2\pi U_1 ), \quad u_2 = exp(-\frac{1}{2} r^2), \quad r = \sqrt{-2ln(u_2)} $

Alumni Liaison

To all math majors: "Mathematics is a wonderfully rich subject."

Dr. Paul Garrett