In [9]:
import numpy as np
import matplotlib.pyplot as plt
import math
import scipy.special as spsp
import scipy.stats as spst
%matplotlib inline
import sympy
From the previous derivation, we know that the best $\lambda$ value is 1, while the best $c$ value is $\sqrt{\frac{2e}{\pi}}$.
Now, let’s start our process to generate samples from a normal distribution.
In [10]:
def pdf_p(x):
return np.exp(-x)
def pdf_t(x):
return 2*1/np.sqrt(2*np.pi)*np.exp(-x**2/2) #normalized distribution
# if we don’t multiply by 2, the plot is problematic
# the area under the organe line is much smaller than under the blue line, and is not close to 1
# multiply 2 is half, so we could get normalized
# unnormalized pdf:1/np.sqrt(2*np.pi)*np.exp(-x**2/2)
c=np.sqrt(2*np.exp(1)/np.pi)
def half_sample():
u=np.random.rand() # uniform distribution between 0 and 1
#-1/lm*log(u) generate a sample from lmdda distribution
proposal=-1/1*np.log(u)
while np.random.rand()>pdf_t(proposal)/pdf_p(proposal)/c:
u=np.random.rand()
#-1/lm*log(u)
proposal=-1/1*np.log(u) #in -1/1*np.log(u): 1 is lmbda
return proposal
half_sample()
Out[10]:
0.585146774858902
In [11]:
x=np.linspace(0,10,1000)
pdf_px=pdf_p(x)
pdf_tx=pdf_t(x)
plt.plot(x,pdf_px,label=”proposal”)
plt.plot(x,pdf_tx,label=”target”)
plt.legend()
plt.show()

In [12]:
#how to generate samples from a normalized distribution
half_samples=[half_sample() for i in range(100)]
half_samples=np.array(half_samples)
half_samples
#convert all positive values that should be negative to negative values
#the reason of this problem is that it is a half normalized distribution
FlipB=np.random.rand(100)<0.5
half_samples[FlipB]*=-1
#or
#half_samples[FlipB]*=-half_samples[FlipB]
sn_samples=half_samples
plt.hist(sn_samples,density=True,bins=30)
Out[12]:
(array([0.0655142 , 0. , 0. , 0.1310284 , 0.1310284 ,
0.1965426 , 0. , 0.1310284 , 0.0655142 , 0.1310284 ,
0.26205681, 0.26205681, 0.26205681, 0.26205681, 0.32757101,
0.58962781, 0.32757101, 0.32757101, 0.0655142 , 0.45859941,
0.1310284 , 0.39308521, 0.52411361, 0.58962781, 0.1965426 ,
0.0655142 , 0.32757101, 0.0655142 , 0.1310284 , 0.1310284 ]),
array([-2.68772923, -2.53509057, -2.38245191, -2.22981325, -2.07717459,
-1.92453592, -1.77189726, -1.6192586 , -1.46661994, -1.31398128,
-1.16134262, -1.00870396, -0.8560653 , -0.70342664, -0.55078797,
-0.39814931, -0.24551065, -0.09287199, 0.05976667, 0.21240533,
0.36504399, 0.51768265, 0.67032131, 0.82295998, 0.97559864,
1.1282373 , 1.28087596, 1.43351462, 1.58615328, 1.73879194,
1.8914306 ]),

In [13]:
sigma=10
mu=5
samples=sigma*sn_samples+mu
In [14]:
x=np.arange(-20,100)
plt.plot(x,spst.norm.pdf(x,mu,sigma))
plt.hist(samples,bins=30,density=True)
plt.show()
# for continuous distribution, it is very normal to have values exceed to 1, because
# it is likelihood instead of prob, prob is very close to 0

In [15]: