Random Number Generation AR and Box Muller

Generate Normal Random Numbers using Box- Muller Transformation

Basic Theory

Let \(X\) and \(Y\) are independent standard normal random variables. Then we know that \(X^2\) and \(Y^2\) follows chi-square distribution each with one degress of freedom.

Cleary \(X^2+Y^2 \sim \chi^2_2\) degrees of freedom. Hence we can see that \(X^2+Y^2=R^2\), where \(R\) is a random number from \(\chi^2\) distribution with \(2\) degrees of freedom. It is known that \(\chi^2\) distribution with \(2\) degrees of freedom is expoential random variable with mean \(2\).

The above equation has resemblance with the polar equation of the circle with radius \(R\). Then \(X\) and \(Y\) can be expressed as \(X=R \;Cos \theta\) and \(Y=R \;Sin \theta\) where \(\theta \sim U(0,2\pi)\).

Generating \(\chi^2\) random Variable

As mentioned earlier, \(\chi^2\) distribution with two degress of freedom is the exponential distribution with mean \(2\). We use the inverse transformation method to generate expoential random numbers. Inverse transformation method relies on the principle that the the distribution function of any continuous random variable follows normal distribution. Hence \[ \begin{align*} 1-e^{-(\frac{x}{2})} =u &\implies& e^{-(\frac{x}{2})}=1-u=u\\ & \implies& -(\frac{x}{2}) = \log u \implies x=-2 \log u \end{align*} \]

import random
import numpy as np
random.seed(2345)
u1=[random.random() for i in range(1000)]
exp_rn=(-1)*2*np.log(u1)

Generating Uniform random number between \((0,\pi)\)

As earlier we will use the inversion method to generate random number following uniform distribution over the range \((a,b)\) \[ \begin{align*} \frac{x-a}{b-a} =u &\implies& x-a=u(b-a)\implies x = a+u(b-a) \end{align*} \] Here \(a=0, b=2\pi\), then \(x=2 \pi \times u\)

import random
import numpy as np
random.seed(2345)
u2=[random.random() for i in range(1000)]
uniform_rn=[2*np.pi *u2[i] for i in range(1000)] 

Generate Normal Random Numbers

Now we will generate independent standard normal random numbers

x=np.sqrt(exp_rn)*np.sin(uniform_rn)
y=np.sqrt(exp_rn)*np.cos(uniform_rn)

Note: Generate Standard Cauchy Random Numbers

We know that the ratio of Independent standard normal variates follow standard cauchy distribution

cauchy_rn=x/y

Generate standard cauchy random number using Inversion Method

The distribution function of the Cauchy distribution is given by

\[ F(x)= \dfrac{1}{2}+\dfrac{1}{\pi}\; \tan^{-1}(x), x \in \mathbb{R}\] \[\begin{align*} & u= \dfrac{1}{2}+\dfrac{1}{\pi}\; \tan^{-1}(x) & \implies & u- \dfrac{1}{2}=\dfrac{1}{\pi}\; \tan^{-1}(x) \implies x=\tan \left[\pi\left(u-\dfrac{1}{2}\right)\right] \end{align*}\]

Generating Cauchy Random Numbers

import random
import numpy as np
random.seed(12234)
u=[random.random() for i in range(1000)]
cauchy_rn=[np.tan(np.pi*(u[i]-0.5)) for i in range(1000)]

Illustration that the sample mean of 100 Cauchy samples are not the converging

import random
import numpy as np
import statistics
cauchy_mean=[]
random.seed(12345)
for i in range(100):
  u=[random.random() for i in range(1000)]
  cauchy_rn=[np.tan(np.pi*(u[i]-0.5)) for i in range(1000)]
  cauchy_mean.append(statistics.mean(cauchy_rn))
cauchy_mean

[0.9561667614399473,
 -0.2741123567390719,
 2.2303576593128147,
 5.348081589223743,
 -12.77540819985642,
 0.33069648237566285,
 -1.8232311064450064,
 7.163378043028272,
 0.9820971521189243,
 -1.789623164515212,
 -0.8812434994598319,
 -1.5836028087275154,
 4.984652887947754,
 0.1319881226284066,
 1.860615292094699,
 1.8242809011596965,
 -2.1168667361858615,
 3.9045346102848972,
 -28.87879928308358,
 2.0441939191075114,
 -2.4759975881063268,
 0.09253005828668526,
 12.038202690521324,
 -0.7304520180555217,
 -0.523421277847376,
 -0.46188794051599513,
 -32.35922019096903,
 1.6473711099280504,
 0.6436813132755944,
 0.24242381253533596,
 -13.359927062101539,
 -1.4844963662809332,
 -0.16864021787431047,
 1.4593258004559004,
 3.9416946235943207,
 -2.499685135599515,
 1.0473697652570162,
 -1.4199315778477133,
 0.3965436329365154,
 -1.3105055236590863,
 2.5167347370321167,
 1.3649779531616952,
 0.6950286962797302,
 0.9687633965462925,
 -4.63193361602274,
 1.28250791377469,
 -1.0020116608747398,
 0.27427294970236926,
 1.2632985403261685,
 0.27384255280254877,
 -0.9560630761996312,
 0.583507970929782,
 4.1785155444940525,
 0.6540086476905878,
 -0.7874642252893138,
 0.10989358558962904,
 10.966242309157717,
 -0.6733637619606582,
 -0.23353576603619502,
 1.085751353240093,
 -0.37909947504027025,
 -0.35331557537427544,
 4.663358874145509,
 -0.8322458687216593,
 -0.9976372106632259,
 0.07857529212697828,
 -0.28293707269030816,
 -0.8725395930048265,
 -1.8778214962670927,
 -1.8475032047855966,
 1.8182727227327773,
 1.285963715749989,
 -31.78251427264197,
 -0.45333164138238063,
 0.03228581810727643,
 0.041870050331878655,
 0.5462010365298042,
 16.22641622928821,
 -0.1599014556485397,
 9.148122239047023,
 -0.6868423305035695,
 2.6504073044348386,
 -0.26359839222408854,
 -0.588203585486739,
 0.22292801767400092,
 2.730346135078863,
 0.4639356885199557,
 0.3936827237134467,
 0.8303399922556017,
 -1.0672483360198113,
 5.470985067283734,
 0.5792096648039795,
 -0.5602544400014225,
 0.7369279245884592,
 2.4954940635061944,
 -5.699852680823704,
 0.49161128299854706,
 0.28176289466605586,
 -1.2038962685832464,
 2.201906349105659]

Acceptance Rejection Sampling

Let \(X\) be the random variable for which we have to generate random numbers and has density \(f(.)\). Let \(g(\cdot)\) be a density function such that \(f(x) \leq Cg(x) \forall x\). Then this \(g(\cdot)\) is called proposal density. The \(g(\cdot)\) must be chosen in such a way that

1. the random numbers from \(g(\cdot)\) can be generated easily

2. \(C\) should be minimum as possible

The steps in the generation of random numbers using acceptance rejection procedure are as follows

1. Generate random numbers of the variable $ Y $ having the density function \(g\) or in other words random numbers of proposal density.
2. Generate random numbers from $ U (0, 1), $ independent of $ Y $ ;
3. If $ U < f(Y)/Cg(Y), $ then set $ X = Y $ ; else go back to step 1.

Generate random numbers of the random varaible having \(f(x)=60x^3(1-x)^2, 0\leq x\leq 1\)

As the range of variation of the random variable \(X\) lies between \(0\) and \(1\), standard uniform distribution is the most appropriate candidate for proposal density \(g(\cdots)\).

Now \[ h(x)=\dfrac{f(x)}{g(c)}=60x^3(1-x)^2\]

Now the maximum value of \(h(x)\) is the value for which the first derivative of \(h(x), h'(x)\) becomes zero and where the second derivative \(h''(x)<0\) \[\begin{align*} & h'(x)=0 &\implies& 60\left[x^3 \times 2 \times (1-x)\times (-1)+3 \times x^2\times(1-x)^2 \right] \\ &&\implies & 60x^2(1-x)\left[-2x+3(1-x)\right]=0\implies 60x^2(1-x)(3-5x)=0\\ && \implies& x=0,x=1,x=\dfrac{3}{5}. \end{align*}\] Clearly \(x=\dfrac{3}{5}\) is the only possible candidate for point of maximum or minimum. Now we shall consider \(h''(x)\). \[\begin{align*} & h''(x)&=& \dfrac{d}{dx}60x^2(1-x)(3-5x)=60\left[2x\times(1-x)(3-5x)- x^2(3-5x)-5x^2(1-x)\right] \\ &h''(x)\vert_{x=\dfrac{3}{5}}&=&60 \times (-5) \times \left(\dfrac{3}{5}\right)^2\left(1-\dfrac{3}{5}\right)<0 \end{align*}\]

Therefore \(h(x)\) is maximum at \(x=\dfrac{3}{5}\).

This implies that \(\dfrac{f(x)}{g(x)}\leq h\left(\dfrac{3}{5}\right)\). \[\begin{align*} &c&=&h\left(\dfrac{3}{5}\right)= 60 x^3(1-x)^2 \vert _{\dfrac{3}{5}}\\&&=&60\times \left(\dfrac{3}{5}\right)^3\times \left(1-\dfrac{3}{5}\right)^2\end{align*}\] Hence we can see that \(\dfrac{f(x)}{g(x)} \leq c \implies f(x) \leq c \times g(x)\)

Now let us plot the density functions.

import matplotlib.pyplot as plt
import numpy as np
plt.clf()
x=np.linspace(0,1,1000)
x.shape[0]
beta_cdf=60*x**3*(1-x)**2
c=60*((3/5)**3)*(1-(3/5))**2
print(f"Value of the constant c is {c}")
std_uniform=[1]*(1000)
cgx=[c]*(1000)
plt.plot(x,beta_cdf,label="Given Density",color='red')
plt.plot(x,std_uniform,label="Standard Uniform Density",color='green')
plt.plot(x,cgx,label="C times assumed distribution",color='blue')
plt.legend()
plt.show()

Value of the constant c is 2.0736

We can see that the probability density function for which random numbers has be to generated(red coloured curve) is less than the c times the uniform distribution (blue coloured line). Now we will move on to generating the random number.

Generate random numbers of proposal desnity (Satndard Uniform Distribution)

import random
y=[random.random() for i in range(1000)]

Find the value of \(c\) and define the function for computing the upperbound of the uniform random numbers as per step 3

c=60*((3/5)**3)*(1-(3/5))**2
def upper_bound(x):
  return((60*(x**3)*(1-x)**2)/c)

Generating new random number and checking condition(Step 3)

import random
u2=[random.random() for i in range(1000)]
beta_ran=[]
for i in range(1000):
  if (u2[i]<upper_bound(y[i])):
    beta_ran.append(y[i])

We can see that in this method some values are discarded hence we do not have 1000 random numbers as we expect

len(beta_ran)

487

Suppose we need to have specified number of random numbers. We know that probability of accepting a particular random number is \(\dfrac{1}{c}\). Hence we will initially generate \([C]+1\) number of random numbers for proposal density and uniform random numbers.

Consider the following programme which generates 2000 beta random numbers

import random
import numpy
size=2000
c=60*((3/5)**3)*(1-(3/5))**2
def upper_bound(x):
  return((60*(x**3)*(1-x)**2)/c)

y=[random.random() for i in range(int(np.ceil(c)*size))]
u2=[random.random() for i in range(int(np.ceil(c)*size))]
beta_ran=[]
j=0
for i in range(int(np.ceil(c)*size)):
  if (u2[i]<upper_bound(y[i])):
    beta_ran.append(y[i])
    j=j+1
    if(j==size):
      break

Plotting the histogram and desnity function

import matplotlib.pyplot as plt
import numpy as np
plt.clf()
plt.hist(beta_ran,bins=30,color='blue',density=True)
x=np.linspace(0,1,1000)
beta_den=60*x**3*(1-x)**2
plt.plot(x,beta_den,color='green',label='Given Density')
plt.show()

Generation of Normal random numbers with mean 10 and variance 25

In this section we will consider generating normal random numbers through acceptance rejection sampling method. The steps involved in generating normal random numbers are given below.

1. Generating random numbers of absolute standard normal random variable acceptance rejection sampling.
2. Generate uniform random number.
3. Generate the random numbers from proposal distribution, exponential distribution.
4. Generate random numbers from standard absolute normal distribution.
5. Generate random numbers of \(-1\) and \(1\) with euqal probabilities \(\frac{1}{2}\)
6. Multiply the random numbers generated in 1.3 and 2 to have random numbers of standard normal distribution.
7. Use \(X=\mu+\sigma*Z\) to obatin the normal random numbers with mean \(\mu\) and standard deviation \(\sigma\).

Basic Theory

1. The density function of absolute standard normal random variable is $ f(x)=e^{-}, x $

Proof: Let $ Z $ be a standard normal random variable and let $ X=|Z| $

\[\begin{align*} &F_X(x) &=& P(X\leq x)= P(\left|Z\right|\leq x) = P(-x \leq Z \leq x)\\ &&=& F_Z(x) - F_Z(-x)\\ &f_X(x) &=& \dfrac{d}{dx}\; F_X(x) = \dfrac{d}{dx}\; \left[F_Z(x) - F_Z(-x)\right]\\ &&=& f_Z(x)- f_z(-x)\times (-1)\\ &&=& \dfrac{1}{\pause \sqrt{2\pi}} e^{- \dfrac{x^2}{2}}+\dfrac{1}{ \sqrt{2\pi}} e^{- \dfrac{(-x)^2}{2}}\\ &&=& \dfrac{2}{\sqrt{2\pi}} e^{- \dfrac{x^2}{2}} \end{align**}\]

Acceptance rejection Sampling Method

For the purpose of generating the absolute standard normal random numbers, we use exponential distribution with mean one as the proposal distribution ie., \(g(x)=e^{-x} x\geq 0\)

Now \[ \begin{align*} &h(x)&=&\frac{f(x)}{g(x)}=\dfrac{\dfrac{2}{\sqrt{2\pi}}e^{- \dfrac{x^2}{2}}} {e^{-x}}\\ & &=& e^{x-\frac{x^2}{2}} \sqrt{\frac{2}{\pi}} \end{align*}\] Now we will try to find out the maximum value of \(h(x)\). From the differential calculus we know that at the point maximumization of a function first derivative will be equal to zero and the second derivative at thet point will be negative. \[ \begin{align*} &h'(x)=0 &\implies& \dfrac{d}{dx}e^{x-\frac{x^2}{2}} \sqrt{\frac{2}{\pi}}=0\\ &&\implies&\sqrt{\frac{2}{\pi}} e^{x-\frac{x^2}{2}} \times\left(1-\dfrac{2x}{2}\right)=0 \implies \sqrt{\frac{2}{\pi}} e^{x-\frac{x^2}{2}} \times\left(1-x\right)=0\\ & &\implies& \left(1-x\right)=0 \implies x=1 \end{align*}\] Now \[ \begin{align*} h''(x)&=& \dfrac{d}{dx}\sqrt{\frac{2}{\pi}} e^{x-\frac{x^2}{2}} (1-x)=\sqrt{\frac{2}{\pi}} e^{x-\frac{x^2}{2}}(-1)+\sqrt{\frac{2}{\pi}} e^{x-\frac{x^2}{2}}(1-x)^2\\ & h''(x)\vert_{x=1}&=& -\sqrt{\frac{2}{\pi}} e^{1-\frac{1}{2}} <0 \end{align*}\]

Hence \(h(x)\) is maximum at \(x=1\) and the maximum value of \(h(x)\) is \[h(1)=\sqrt{\frac{2}{\pi}} e^{1-\frac{1}{2}}=\sqrt{\frac{2}{\pi}} e^{\frac{1}{2}}=\sqrt{\frac{2e}{\pi}}\]

Now $= e^{} $

Generating the random numbers from proposal distribution

We will generate the random numbers from standard exponential distribution using the inversion method. To ensure that we have suffient number of random numbers at the end, as we loose certain random numbers during rejection phase, we will generate \([c]+1\) times random numbers of the proposal distribution.

import random
import numpy as np
num_rn=1000
c=np.sqrt(2*np.exp(1)/np.pi)
u1=[random.random() for i in range((int(c)+1)*num_rn)]
y=(-1)*np.log(u1)

Acceptance Rejection Phase

u2=[random.random() for i in range((int(c)+1)*num_rn)]
j=0
i=0
abs_std_norm_rn=[]
while (j<1000):
  if u2[i]<=np.exp(((y[i]-1)**2)/2):
    abs_std_norm_rn.append(y[i])
    j=j+1
  i=i+1

Generating dichotomous random number

import random
u3=[random.random() for i in range(size)]
ber_rn=[]
for i in range(size):
  if u3[i]<0.5:
    ber_rn.append(-1)
  else:
    ber_rn.append(1)

Generating normal random numbers

# Standard normal random numbers
std_norm_rn=[abs_std_norm_rn[i]*ber_rn[i] for i in range(num_rn)]
mu=10
sigma=5
norm_rn=[mu+sigma*std_norm_rn[i] for i in range(num_rn)]

Plotting the histogram and density

import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as stat
x=np.linspace(min(norm_rn),max(norm_rn),100)
plt.hist(norm_rn, bins=40,density=True)
plt.plot(x,stat.norm.pdf(x,loc=mu,scale=sigma),label='Normal Density',color='red')
plt.legend()
plt.show()

Reference

1. Ross. SM, Simulation(4th Edition), Elsevier Academic Press, Amsterdam,2006