Generating Random Numbers from Gaussian Distribution¶
Faisal Qureshi
faisal.qureshi@ontariotechu.ca
Copyright information¶
© Faisal Qureshi
License¶
This work is licensed under a Creative Commons Attribution-NonCommercial 4.0 International License.
In [1]:
import sys
import numpy as np
import matplotlib.pyplot as plt
Uniform distribution¶
In [2]:
x = np.random.rand(100000)
plt.figure(figsize=(4,4))
plt.title('Uniform distribution')
plt.hist(x, bins=np.linspace(0,1,100));
Gaussian (Normal) distribution¶
In [3]:
mu = -2
sigma = 3
x = np.random.randn(100000) * sigma + mu
plt.figure(figsize=(4,4))
plt.title(f'Normal distribution: $\\mu$={mu}, $\\sigma$={sigma}')
plt.hist(x, bins=np.linspace(-10,10,100));
Playing with $\mu$ and $\sigma$¶
In [4]:
m = [0, 0, 3]
s = [1, 3, 2]
c = ['red', 'blue', 'green']
x = []
for i in range(3):
x.append(np.random.randn(100000) * s[i] + m[i])
plt.figure(figsize=(4,4))
for i in range(3):
label_str = f'mu={m[i]}, sigma={s[i]}$'
plt.hist(x[i], bins=np.linspace(-10,10,100), color=c[i%3], alpha=0.2, label=label_str);
plt.legend()
Out[4]:
<matplotlib.legend.Legend at 0x11669f8c0>
Polar method¶
A method for generating samples for Normal distribution with zero mean and unit variance
In [5]:
def polar_method():
"""
Returns:
- a tuple containing two samples.
- True if some samples were rejected, False otherwise
"""
nloop = 1000
while(nloop > 0):
u = np.random.rand(2)
v = 2*u - 1
s = np.sum(np.square(v))
if s >= 1:
nloop = nloop - 1
continue
x = v * np.sqrt(-2*np.log(s)/s)
break
return x, nloop<1000
Polar method generates two numbers¶
In [6]:
(n1, n2), rejection = polar_method()
print(f' sample 1: {n1}\n sample 2: {n2}\n rejection: {rejection}')
sample 1: 0.04396358153992705 sample 2: -0.027219382319545633 rejection: True
In [7]:
tries = 10000
number_of_rejections = 0
for i in range(tries):
_, rejection = polar_method()
if rejection:
number_of_rejections += 1
print(f' tries: {tries}\n rejections: {number_of_rejections}')
tries: 10000 rejections: 2156
Plotting samples generated using Polar method¶
In [8]:
def randn_polar_method(n):
x = np.empty(n)
for i in range(1, n//2):
(n1, n2), _ = polar_method()
x[2*i] = n1
x[2*i + 1] = n2
return x
In [9]:
mu = 0
sigma = 3
x = randn_polar_method(100000) * sigma + mu
plt.figure(figsize=(4,4))
plt.title(f'Uniform distribution (Polar method): $\\mu$={mu}, $\\sigma$={sigma}')
plt.hist(x, bins=np.linspace(-10,10,100));
Ratio method¶
A method for generating samples from Normal distribution with zero mean and unit variance
In [10]:
def ratio_method():
"""
Generates a single sample
"""
nloop = 1000
while(nloop > 0):
u = np.random.rand(2)
x = np.sqrt(8/(2*np.exp(1)))*(u[0]-0.5)/u[1]
if x*x <= 5 - 4*np.exp(1/4)*u[1]:
return x, nloop < 1000
elif x*x >= 4*np.exp(-1.35)/u[1] + 1.4:
nloop = nloop - 1
continue
elif x*x <= -4/np.log(u[1]):
return x, nloop < 1000
else:
nloop = nloop - 1
continue
return 0
Ratio method generates a single number¶
In [11]:
ratio_method()
Out[11]:
(np.float64(0.7390788687671123), False)
Ploting samples generated using Ratio method¶
In [12]:
def randn_ratio_method(n):
x = np.empty(n)
for i in range(1, n):
n, _ = ratio_method()
x[i] = n
return x
In [13]:
mu = 0
sigma = 3
x = randn_ratio_method(1000000) * sigma + mu
plt.figure(figsize=(4,4))
plt.title(f'Uniform distribution (Polar method): $\\mu$={mu}, $\\sigma$={sigma}')
plt.hist(x, bins=np.linspace(-10,10,100));
Exponential distribution¶
Samples
In [14]:
def plot_exponential(n):
mu = 1
x = np.empty(n)
for i in range(0,n):
u = - mu * np.log( np.random.rand() )
x[i] = u
return x
In [15]:
n = 100000
x = plot_exponential(n)
plt.figure(figsize=(4,4))
plt.title(f'Exponential')
plt.hist(x, bins=np.linspace(0,10,100));
In [ ]: