Probability distributions
At the heart of probabilistic methods and Bayesian decision making are the probability
distributions and and therefore we need to study them to gain true power to this
magic.
Bernoulli distribution
Bernoulli distribution is defined as
and its parameter (probability of ) is estimated from
but where that comes from?
P(x) p(x)
Bern(x; μ) = μx(1 − μ)1−x , (1)
μ x = 1
μ =
N
∑
n=1
xn , (2)
1
N
In [1]:
import matplotlib.pyplot as plt
import numpy as np
#
# 1. Bernoulli parameter estimation (ML)
np.random.seed(666) # to always get the same points
mu = 0.3
N_seq = np.array([1,2,4,6,8,10,12,14,16,18,20,25,30,35,40,45,50])
mu_ML = np.empty(N_seq.shape)
a=4
b=6
a2=40
b2=60
mu_MAP = np.empty(N_seq.shape)
mu_MAP2 = np.empty(N_seq.shape)
for N_num, N in enumerate(N_seq):
x_n = np.random.binomial(1,mu,size=N)
mu_ML[N_num] = sum(x_n)/N
mu_MAP[N_num] = (sum(x_n)+a)/(N+a+b)
mu_MAP2[N_num] = (sum(x_n)+a2)/(N+a2+b2)
#
# Plot 1 – only ML
fig = plt.figure(figsize =(8, 4))
plt.plot(N_seq, mu_ML,’c-‘,label=’mu ML’)
plt.plot([0,N],[mu,mu],’k–‘)
plt.xticks(N_seq)
plt.xlabel(‘N’)
plt.ylabel(‘mu’)
plt.legend()
#plt.gcf().subplots_adjust(top=0.8)
plt.show()
lec05_probability_distributions http://localhost:8888/nbconvert/html/lec05_probability_…
1 of 21 9/13/21, 16:17
is estimated from
Gaussian distribution
μMAP
μMAP = . (3)
m + a
N + a + b
In [2]:
#
# Plot 2 – ML plus 2x MAP
fig = plt.figure(figsize =(8, 4))
plt.plot(N_seq, mu_ML,’c-‘,label=’mu ML’)
plt.plot(N_seq, mu_MAP,’m-‘,label=’mu MAP (a=4,b=6)’)
plt.plot(N_seq, mu_MAP2,’y-‘,label=’mu MAP (a=40,b=60)’)
plt.plot([0,N],[mu,mu],’k–‘)
plt.xticks(N_seq)
plt.xlabel(‘N’)
plt.ylabel(‘mu’)
plt.legend()
plt.show()
N (x; μ, σ2) = e
−
, (4)
1
√2πσ2
(x−μ)2
2σ2
lec05_probability_distributions http://localhost:8888/nbconvert/html/lec05_probability_…
2 of 21 9/13/21, 16:17
Poisson distribution
Notably:
2D Gaussian
Let’s assume we have two dimensional measurements , for example is
height and is weight.
In [3]:
# 1. Bernoulli parameter estimation (ML)
np.random.seed(666) # to always get the same points
x = np.linspace(-6.0,6.0,101)
mu_1 = 0.0
sigma2_1 = 2
gauss_1 = 1/np.sqrt(2*np.pi*sigma2_1)*np.exp(-1/(2*sigma2_1)*(x-mu_1)**2)
mu_2 = 2.0
sigma2_2 = 2
gauss_2 = 1/np.sqrt(2*np.pi*sigma2_2)*np.exp(-1/(2*sigma2_2)*(x-mu_2)**2)
mu_3 = 2.0
sigma2_3 = 4
gauss_3 = 1/np.sqrt(2*np.pi*sigma2_3)*np.exp(-1/(2*sigma2_3)*(x-mu_3)**2)
#
# Plot 1
plt.plot(x, gauss_1,’c-‘,label=’mu=0.0, sigma2=2.0’)
plt.plot(x, gauss_2,’m-‘,label=’mu=2.0, sigma2=2.0’)
plt.plot(x, gauss_3,’y-‘,label=’mu=2.0, sigma2=4.0’)
plt.xlabel(‘x’)
plt.ylabel(‘y’)
plt.legend()
plt.show()
p(x = k; λ) = (5)
λke−λ
k!
λ = E(x) = V ar(x) (6)
→x = (x1, x2)T x1
x2
lec05_probability_distributions http://localhost:8888/nbconvert/html/lec05_probability_…
3 of 21 9/13/21, 16:17
In [4]:
#
# Let’s make 2D correlated data of hobbits and elfs
# Hobits
np.random.seed(42) # to always get the same points
N_h = 80
x1_h = np.random.normal(1.1,0.3,N_h)
a_h = 50.0
b_h = 20.0
x2_h_noise = np.random.normal(0,8,N_h)
x2_h = a_h*x1_h+b_h+x2_h_noise
# Elves
N_e = 20
x1_e = np.random.normal(1.9,0.4,N_e)
a_e = 30.0
b_e = 30.0
x2_e_noise = np.random.normal(0,8,N_e)
x2_e = a_e*x1_e+b_e+x2_e_noise
#
# PLOT 1 both classes in 1D
plt.plot(x1_h,np.zeros([N_h,1]),’co’, label=”hobit”)
plt.plot(x1_e,np.zeros([N_e,1]),’mo’, label=”elf”)
plt.title(‘Training samples from two classes c1 and c2’)
plt.legend()
plt.xlabel(‘height [m]’)
plt.show()
plt.hist(x1_h, bins = 10, alpha=0.5, color=’cyan’)
plt.hist(x1_e, bins = 10, alpha=0.5, color=’magenta’)
plt.title(‘Training samples from two classes c1 and c2’)
plt.xlabel(‘height [m]’)
plt.ylabel(‘number’)
plt.show()
#
# PLOT 2 both classes in 1D
plt.plot(x2_h,np.zeros([N_h,1]),’co’, label=”hobitti”)
plt.plot(x2_e,np.zeros([N_e,1]),’mo’, label=”haltija”)
plt.title(‘Training samples from two classes c1 and c2’)
plt.legend()
plt.xlabel(‘weight [kg]’)
plt.show()
plt.hist(x2_h, bins = 10, alpha=0.5, color=’cyan’)
plt.hist(x2_e, bins = 10, alpha=0.5, color=’magenta’)
plt.title(‘Training samples from two classes c1 and c2’)
plt.xlabel(‘weight [kg]’)
plt.ylabel(‘number’)
plt.show()
#
# PLOT 2 both classes in 2D
plt.plot(x1_h,x2_h,’co’, label=”hobitti”)
plt.plot(x1_e,x2_e,’mo’, label=”haltija”)
plt.title(‘Training samples from two classes c1 and c2’)
plt.legend()
plt.xlabel(‘height [m]’)
plt.ylabel(‘weight [kg]’)
#plt.axis([0.5,2.5,-1.1,+1.1])
lec05_probability_distributions http://localhost:8888/nbconvert/html/lec05_probability_…
4 of 21 9/13/21, 16:17
x2_h_test_noise = np.random.normal(0,8,N_h_test)
x2_h_test = a_h*x1_h_test+b_h+x2_h_test_noise
c_h_test = np.zeros(N_h_test)
c_h_test[:] = 1
N_e_test = 10
x1_e_test = np.random.normal(1.9,0.4,N_e_test)
x2_e_test_noise = np.random.normal(0,8,N_e_test)
x2_e_test = a_e*x1_e_test+b_e+x2_e_test_noise
c_e_test = np.zeros(N_e_test)
c_e_test[:] = 2
plt.plot(x1_h,x2_h,’co’, label=”hobit”)
plt.plot(x1_e,x2_e,’mo’, label=”elf”)
plt.plot(x1_h_test,x2_h_test,’ko’, label=”test (hob)”)
plt.plot(x1_e_test,x2_e_test,’kd’, label=”test (elf)”)
plt.title(‘Training and test samples from two classes c1 ja c2’)
plt.legend()
plt.xlabel(‘height [m]’)
plt.ylabel(‘weight [kg]’)
#plt.axis([0.5,2.5,-1.1,+1.1])
plt.show()
x1_test = np.concatenate((x1_h_test,x1_e_test))
x2_test = np.concatenate((x2_h_test,x2_e_test))
c_test = np.concatenate((c_h_test,c_e_test))
c_test_hat = np.zeros(c_test.shape)
priori_h = N_h/(N_h+N_e)
mu1_h = np.mean(x1_h)
mu2_h = np.mean(x2_h)
sigma2_1_h = np.var(x1_h)
sigma2_2_h = np.var(x2_h)
p_h_1 = 1/np.sqrt(2*np.pi*sigma2_1_h)*np.exp(-1/(2*sigma2_1_h)*(x1_test-mu1_h
P_h_1 = priori_h*p_h_1
p_h_2 = 1/np.sqrt(2*np.pi*sigma2_2_h)*np.exp(-1/(2*sigma2_2_h)*(x2_test-mu2_h
P_h_2 = priori_h*p_h_2
P_h = priori_h*p_h_1*p_h_2
priori_e = N_e/(N_h+N_e)
mu1_e = np.mean(x1_e)
mu2_e = np.mean(x2_e)
sigma2_1_e = np.var(x1_e)
sigma2_2_e = np.var(x2_e)
p_e_1 = 1/np.sqrt(2*np.pi*sigma2_1_e)*np.exp(-1/(2*sigma2_1_e)*(x1_test-mu1_e
P_e_1 = priori_e*p_e_1
p_e_2 = 1/np.sqrt(2*np.pi*sigma2_2_e)*np.exp(-1/(2*sigma2_2_e)*(x2_test-mu2_e
P_e_2 = priori_e*p_e_2
P_e = priori_e*p_e_1*p_e_2
c_test_hat[np.argwhere(np.greater(P_h_1,P_e_1) == True)] = 1
c_test_hat[np.argwhere(np.greater(P_h_1,P_e_1) == False)] = 2
corr = np.count_nonzero(np.equal(c_test,c_test_hat))
success_rate corr/(N_h_test+N_e_test)
lec05_probability_distributions http://localhost:8888/nbconvert/html/lec05_probability_…
5 of 21 9/13/21, 16:17
lec05_probability_distributions http://localhost:8888/nbconvert/html/lec05_probability_…
6 of 21 9/13/21, 16:17
We can now use 2D Gaussian distribution
lec05_probability_distributions http://localhost:8888/nbconvert/html/lec05_probability_…
7 of 21 9/13/21, 16:17
N (→x; →μ, Σ) = e− (→x−→μ)
T Σ−1(→x−→μ) . (7)
1
2π|Σ|
1
2
1
2
lec05_probability_distributions http://localhost:8888/nbconvert/html/lec05_probability_…
8 of 21 9/13/21, 16:17
In [5]:
from numpy.linalg import inv
from mpl_toolkits import mplot3d
# Hobits
np.random.seed(42) # to always get the same points
N_h = 80
x1_h = np.random.normal(1.1,0.3,N_h)
a_h = 50.0
b_h = 20.0
x2_h_noise = np.random.normal(0,8,N_h)
x2_h = a_h*x1_h+b_h+x2_h_noise
# Elves
N_e = 20
x1_e = np.random.normal(1.9,0.4,N_e)
a_e = 30.0
b_e = 30.0
x2_e_noise = np.random.normal(0,8,N_e)
x2_e = a_e*x1_e+b_e+x2_e_noise
print(‘Mean and covariance for hobits’)
X_h = np.concatenate(([x1_h],[x2_h]),axis=0)
mu_h = np.mean(X_h,axis=1)
print(mu_h)
Sigma_h = np.cov(X_h)
print(Sigma_h)
priori_h = N_h/(N_h+N_e)
print(‘Mean and covariance for elfs’)
X_e = np.concatenate(([x1_e],[x2_e]),axis=0)
mu_e = np.mean(X_e,axis=1)
print(mu_e)
Sigma_e = np.cov(X_e)
print(Sigma_e)
priori_e = N_e/(N_h+N_e)
def gaussian2d(X,Y, mu, Sigma, priori):
Z = np.zeros(X.shape)
for i in range(X.shape[0]):
for j in range(X.shape[1]):
Z[i,j] = priori*1/np.sqrt(2*np.pi)*1/np.sqrt(np.linalg.det(Sigma
return Z
x_1 = np.linspace(0,3.0,50)
x_2 = np.linspace(0,150.0,50)
X, Y = np.meshgrid(x_1,x_2)
Z_h = gaussian2d(X,Y, mu_h, Sigma_h, priori_h)
Z_e = gaussian2d(X,Y, mu_e, Sigma_e, priori_e)
fig = plt.figure()
ax = plt.axes(projection=’3d’)
ax.contour3D(X, Y, Z_h, 50, cmap=’Greens’)
ax.contour3D(X, Y, Z_e, 50, cmap=’Reds’)
ax.set_xlabel(‘height [m]’)
ax.set_ylabel(‘weight [kg]’)
ax.set_zlabel(‘Probability density’);
plt.legend()
plt.show()
lec05_probability_distributions http://localhost:8888/nbconvert/html/lec05_probability_…
9 of 21 9/13/21, 16:17
Mean and covariance for hobits
[ 1.06284795 73.05149079]
[[8.25096665e-02 4.10250733e+00]
[4.10250733e+00 2.57865233e+02]]
Mean and covariance for elfs
[ 2.036617 89.43100536]
[[ 0.16147336 4.32103664]
[ 4.32103664 148.50026465]]
No handles with labels found to put in legend.
[1.08673514e-01 5.31471966e-02 3.87556653e-02 5.42005210e-02
6.04980863e-02 7.47124289e-02 8.63787750e-02 1.17863114e-01
1.06417861e-01 2.69055703e-05 1.56804323e-20 4.48521482e-02
c_h_test = np.zeros(N_h_test)
c_h_test[:] = 1
N_e_test = 10
x1_e_test = np.random.normal(1.9,0.4,N_e_test)
x2_e_test_noise = np.random.normal(0,8,N_e_test)
x2_e_test = a_e*x1_e_test+b_e+x2_e_test_noise
c_e_test = np.zeros(N_e_test)
c_e_test[:] = 2
x1_test = np.concatenate((x1_h_test,x1_e_test))
x2_test = np.concatenate((x2_h_test,x2_e_test))
c_test = np.concatenate((c_h_test,c_e_test))
c_test_hat = np.zeros(c_test.shape)
P_h = np.zeros(x1_test.shape)
P_e = np.zeros(x1_test.shape)
for i in range(x1_test.shape[0]):
P_h[i] = priori_h*1/np.sqrt(2*np.pi)*1/np.sqrt(np.linalg.det(Sigma_h))*
P_e[i] = priori_e*1/np.sqrt(2*np.pi)*1/np.sqrt(np.linalg.det(Sigma_e))*
print(P_h)
c_test_hat[np.argwhere(np.greater(P_h,P_e) == True)] = 1
c_test_hat[np.argwhere(np.greater(P_h,P_e) == False)] = 2
corr = np.count_nonzero(np.equal(c_test,c_test_hat))
success_rate = corr/(N_h_test+N_e_test)
print(f’Success rate using both: {success_rate:.2f}’)
lec05_probability_distributions http://localhost:8888/nbconvert/html/lec05_probability_…
10 of 21 9/13/21, 16:17
4.09119894e-09 2.77936689e-03 7.71264106e-03 9.89061138e-10
1.91243096e-11 3.06220284e-03 2.35139734e-05 1.49007940e-07]
Success rate using both: 0.90
D-dimensional Gaussian
Naive Bayes classifier
Gaussian mixture model
where
N (→x; →μ, Σ) = e− (→x−→μ)
T Σ−1(→x−→μ) (8)
1
(2π) |Σ|
D
2
1
2
1
2
Σ3×3 =
⎡⎢⎢⎣
σ21 σ1σ2 σ1σ3
σ2σ1 σ
2
2 σ2σ2
σ3σ1 σ3σ2 σ
2
3
⎤⎥⎥⎦
Σnaive3×3 =
⎡⎢⎢⎣
σ21 0 0
0 σ22 0
0 0 σ23
⎤⎥⎥⎦
. (9)
P(→x; {αc, →μc, Σc}c=1,…,C) =
C
∑
c=1
αcN (→x; →μc, Σc) . (10)
C
∑
c=1
αc = 1 . (11)
lec05_probability_distributions http://localhost:8888/nbconvert/html/lec05_probability_…
11 of 21 9/13/21, 16:17
In [6]:
from sklearn.mixture import GaussianMixture
#
# Let’s make 2D correlated data of hobbits and elfs
# Hobits
np.random.seed(42) # to always get the same points
N_h = 80
x1_h = np.random.normal(1.1,0.3,N_h)
a_h = 50.0
b_h = 20.0
x2_h_noise = np.random.normal(0,8,N_h)
x2_h = a_h*x1_h+b_h+x2_h_noise
# Elves
N_e = 20
x1_e = np.random.normal(1.9,0.4,N_e)
a_e = 30.0
b_e = 30.0
x2_e_noise = np.random.normal(0,8,N_e)
x2_e = a_e*x1_e+b_e+x2_e_noise
# Dwarfs
N_d = 60
x1_d = np.random.normal(1.0,0.1,N_d)
a_d = 50.0
b_d = 40.0
x2_d_noise = np.random.normal(0,6,N_d)
x2_d = a_d*x1_d+b_d+x2_d_noise
print(‘Mean and covariance for hobits’)
X_h = np.concatenate(([x1_h],[x2_h]),axis=0)
mu_h = np.mean(X_h,axis=1)
print(mu_h)
Sigma_h = np.cov(X_h)
print(Sigma_h)
priori_h = N_h/(N_h+N_e+N_d)
print(‘Mean and covariance for elfs’)
X_e = np.concatenate(([x1_e],[x2_e]),axis=0)
mu_e = np.mean(X_e,axis=1)
print(mu_e)
Sigma_e = np.cov(X_e)
print(Sigma_e)
priori_e = N_e/(N_h+N_e+N_d)
print(‘Mean and covariance for dwarfs’)
X_d = np.concatenate(([x1_d],[x2_d]),axis=0)
mu_d = np.mean(X_d,axis=1)
print(mu_d)
Sigma_d = np.cov(X_d)
print(Sigma_d)
priori_d = N_d/(N_h+N_e+N_d)
print(‘Mean and covariance for dwarfs and elves combined’)
X_ed = np.concatenate((X_d, X_e),axis=1)
mu_ed = np.mean(X_ed,axis=1)
print(mu_ed)
Sigma_ed = np.cov(X_ed)
print(Sigma_ed)
priori_ed = (N_d+N_e)/(N_h+N_e+N_d)
lec05_probability_distributions http://localhost:8888/nbconvert/html/lec05_probability_…
12 of 21 9/13/21, 16:17
x_1 = np.linspace(0,3.0,70)
x_2 = np.linspace(0,150.0,70)
X, Y = np.meshgrid(x_1,x_2)
Z_h = gaussian2d(X,Y, mu_h, Sigma_h, priori_h)
Z_e = gaussian2d(X,Y, mu_e, Sigma_e, priori_e)
Z_d = gaussian2d(X,Y, mu_d, Sigma_d, priori_d)
Z_ed = gaussian2d(X,Y, mu_ed, Sigma_ed, priori_ed)
fig = plt.figure()
ax = plt.axes(projection=’3d’)
ax.contour3D(X, Y, Z_h, 50, cmap=’Greens’)
ax.contour3D(X, Y, Z_e, 50, cmap=’Reds’)
ax.contour3D(X, Y, Z_d, 50, cmap=’Blues’)
ax.set_xlabel(‘height [m]’)
ax.set_ylabel(‘weight [kg]’)
ax.set_zlabel(‘Probability density’);
plt.legend()
plt.show()
fig = plt.figure()
ax = plt.axes(projection=’3d’)
ax.contour3D(X, Y, Z_h, 50, cmap=’Greens’)
ax.contour3D(X, Y, Z_ed, 50, cmap=’Greys’)
ax.set_xlabel(‘height [m]’)
ax.set_ylabel(‘weight [kg]’)
ax.set_zlabel(‘Probability density’);
plt.legend()
plt.show()
#
# Run Gaussian mixture model
gm = GaussianMixture(n_components=2, random_state=0).fit(X_ed.T)
mu_ed1 = gm.means_[0]
mu_ed2 = gm.means_[1]
Sigma_ed1 = gm.covariances_[0]
Sigma_ed2 = gm.covariances_[0]
priori_ed1 = gm.weights_[0]/(gm.weights_[0]+gm.weights_[1])*priori_ed
priori_ed2 = gm.weights_[1]/(gm.weights_[0]+gm.weights_[1])*priori_ed
Z_ed_1 = gaussian2d(X,Y, mu_ed1, Sigma_ed1, priori_ed1)
Z_ed_2 = gaussian2d(X,Y, mu_ed2, Sigma_ed2, priori_ed2)
fig = plt.figure()
ax = plt.axes(projection=’3d’)
ax.contour3D(X, Y, Z_h, 50, cmap=’Greens’)
ax.contour3D(X, Y, Z_ed_1, 50, cmap=’Greys’)
ax.contour3D(X, Y, Z_ed_2, 50, cmap=’Greys’)
ax.set_xlabel(‘height [m]’)
ax.set_ylabel(‘weight [kg]’)
ax.set_zlabel(‘Probability density’);
plt.legend()
plt.show()
#
lec05_probability_distributions http://localhost:8888/nbconvert/html/lec05_probability_…
13 of 21 9/13/21, 16:17
Mean and covariance for hobits
[ 1.06284795 73.05149079]
[[8.25096665e-02 4.10250733e+00]
[4.10250733e+00 2.57865233e+02]]
Mean and covariance for elfs
[ 2.036617 89.43100536]
[[ 0.16147336 4.32103664]
[ 4.32103664 148.50026465]]
Mean and covariance for dwarfs
[ 1.01689026 91.03087247]
[[1.20847128e-02 5.89419803e-01]
[5.89419803e-01 6.10677983e+01]]
Mean and covariance for dwarfs and elves combined
[ 1.27182195 90.63090569]
[[ 0.24529913 1.16967159]
[ 1.16967159 81.80884492]]
No handles with labels found to put in legend.
c_e_test[:] = 2
x1_test = np.concatenate((x1_h_test,x1_e_test))
x2_test = np.concatenate((x2_h_test,x2_e_test))
c_test = np.concatenate((c_h_test,c_e_test))
c_test_hat = np.zeros(c_test.shape)
P_h = np.zeros(x1_test.shape)
P_e = np.zeros(x1_test.shape)
P_ed = np.zeros(x1_test.shape)
for i in range(x1_test.shape[0]):
P_h[i] = priori_h*1/np.sqrt(2*np.pi)*1/np.sqrt(np.linalg.det(Sigma_h))*
P_e[i] = priori_e*1/np.sqrt(2*np.pi)*1/np.sqrt(np.linalg.det(Sigma_e))*
P_ed[i] = priori_ed*1/np.sqrt(2*np.pi)*1/np.sqrt(np.linalg.det(Sigma_ed
print(P_h)
print(P_e)
print(P_ed)
print(np.greater(P_h,P_ed))
c_test_hat[np.argwhere(np.greater(P_h,P_ed) == True)] = 1
c_test_hat[np.argwhere(np.greater(P_h,P_ed) == False)] = 2
corr = np.count_nonzero(np.equal(c_test,c_test_hat))
success_rate = corr/(N_h_test+N_e_test)
print(f’Success rate using mixed: {success_rate:.2f}’)
lec05_probability_distributions http://localhost:8888/nbconvert/html/lec05_probability_…
14 of 21 9/13/21, 16:17
No handles with labels found to put in legend.
No handles with labels found to put in legend.
[9.21197738e-02 4.51291121e-02 7.24522596e-02 5.08163971e-03
3.50946409e-02 4.88350601e-02 2.63797425e-02 5.67671207e-02
6.24075010e-02 5.82669904e-02 7.69816839e-04 2.51620559e-05
3.88315692e-06 2.55153307e-05 3.35142767e-05 1.05801437e-09
5.20983222e-03 2.19385296e-03 4.56587802e-04 1.59562422e-07]
[2.49720838e-04 1.51546861e-03 1.79432288e-04 4.13142292e-05
2.57817248e-05 1.13016183e-04 7.95754816e-04 2.88086788e-04
2.34812202e-05 7.14214666e-04 9.94704637e-03 1.43562418e-03
7.35897431e-03 1.50357980e-02 1.55558129e-02 1.13275420e-02
4.38554046e-03 4.18114072e-03 1.09394373e-02 1.13957439e-02]
[1.31698789e-02 1.29562571e-02 4.84565865e-04 1.95044274e-02
1.77551333e-03 2.55549665e-05 3.50023330e-02 3.42411080e-02
lec05_probability_distributions http://localhost:8888/nbconvert/html/lec05_probability_…
15 of 21 9/13/21, 16:17
Non-parametric probability density models
For samples we defineN xn
P(x) =
N
∑
n=1
k ( ) =
N
∑
n=1
N (xn; x, h) (12)
1
N
1
h
x − xn
h
lec05_probability_distributions http://localhost:8888/nbconvert/html/lec05_probability_…
16 of 21 9/13/21, 16:17
In [7]:
#
# 1. Generate and plot random points for training
np.random.seed(66) # to always get the same points
N_h = 1000
N_e = 200
mu_h_gt = 1.1
mu_e_gt = 1.9
sigma_h_gt = 0.3
sigma_e_gt = 0.4
x_h = np.random.normal(mu_h_gt,sigma_h_gt,N_h)
x_e = np.random.normal(mu_e_gt,sigma_e_gt,N_e)
plt.plot(x_h,np.zeros([N_h,1]),’co’, label=”hobit”)
plt.plot(x_e,np.zeros([N_e,1]),’mo’, label=”elf”)
plt.title(f'{N_h} hobit height’)
plt.legend()
plt.xlabel(‘height [m]’)
plt.axis([0.0,3.0,-1.1,+1.1])
plt.show()
#
#
mu_h_est = np.mean(x_h)
mu_e_est = np.mean(x_e)
sigma_h_est = np.std(x_h)
sigma_e_est = np.std(x_e)
print(f’Avg height hobits {mu_h_est:0.2f} (GT: {mu_h_gt:0.2f})’)
print(f’Avg height elves {mu_e_est:0.2f} (GT: {mu_e_gt:0.2f})’)
print(f’Height st. deviation hobits {sigma_h_est:0.3f} (GT: {sigma_h_gt:0.3f
print(f’Height st. deviation elves {sigma_e_est:0.3f} (GT: {sigma_e_gt:0.3f
def gaussian(x, mu, sigma, priori):
z = np.zeros(x.shape)
for i in range(x.shape[0]):
z[i] = priori*1/np.sqrt(2*np.pi)*1/sigma*np.exp((-1/2*(x[i]-mu)**2)
return z
#
#
[x, step_size] = np.linspace(0,3.0,70,retstep=True)
lhood_h_est = gaussian(x, mu_h_est, sigma_h_est, 1)
lhood_e_est = gaussian(x, mu_e_est, sigma_e_est, 1)
lhood_h_gt = gaussian(x, mu_h_gt, sigma_h_gt, 1)
lhood_e_gt = gaussian(x, mu_e_gt, sigma_e_gt, 1)
plt.plot(x_h,np.zeros([N_h,1]),’co’, label=”hobit”)
plt.plot(x_e,np.zeros([N_e,1]),’mo’, label=”elf”)
plt.plot(x,lhood_h_gt,’c-‘, label=”hobit (GT)”)
plt.plot(x,lhood_e_gt,’m-‘, label=”elf (GT)”)
plt.plot(x,lhood_h_est,’c–‘, label=”hobit (est)”)
plt.plot(x,lhood_e_est,’m–‘, label=”elf (est)”)
plt.legend()
plt.xlabel(‘pituus [m]’)
#plt.axis([0.0,3.0,-1.1,+5])
plt.show()
#
#
kern_width = 0.2
[x, step_size] = np.linspace(0,3.0,70,retstep=True)
plt.plot(x_h,np.zeros([N_h,1]),’co’, label=”hobit”)
plt.plot(x_e,np.zeros([N_e,1]),’mo’, label=”elf”)
[x_kern, step_size_kern] = np.linspace(0,3.0,11,retstep=True)
lec05_probability_distributions http://localhost:8888/nbconvert/html/lec05_probability_…
17 of 21 9/13/21, 16:17
#plt.axis([0.0,3.0,-1.1,+5])
plt.show()
# Output value is Gaussian kernel multiplied by all positive samples
lhood_h_est_kern = np.zeros(len(x))
for xind, xval in enumerate(x):
lhood_h_est_kern[xind] = sum(gaussian(x_h,xval,kern_width,1))/len(x_h)
#lhood_h_est_kern[xind] = sum(stats.norm.pdf(x_h, xval, kernel_width))
plt.plot(x_h,np.zeros([N_h,1]),’co’, label=”hobit”)
plt.plot(x_e,np.zeros([N_e,1]),’mo’, label=”elf”)
plt.plot(x,lhood_h_gt,’c-‘, label=”hobit (GT)”)
plt.plot(x,lhood_h_est_kern,’c–‘, label=”hobit (est)”)
#plt.plot(x,lhood_e_est,’m–‘, label=”haltija (est)”)
plt.legend()
plt.xlabel(‘height [m]’)
#plt.axis([0.0,3.0,-1.1,+5])
plt.show()
#
#
kern_width = 0.02
[x, step_size] = np.linspace(0,3.0,70,retstep=True)
plt.plot(x_h,np.zeros([N_h,1]),’co’, label=”hobit”)
plt.plot(x_e,np.zeros([N_e,1]),’mo’, label=”elf”)
[x_kern, step_size_kern] = np.linspace(0,3.0,11,retstep=True)
x_kern_plot = np.linspace(-1.0,+4.0,201)
for foo_ind, foo_val in enumerate(x_kern):
foo_kern = gaussian(x_kern_plot, foo_val, kern_width, 1)
plt.plot(x_kern_plot, foo_kern,’y–‘,label=’kerneli’)
break
plt.legend()
plt.xlabel(‘height [m]’)
#plt.axis([0.0,3.0,-1.1,+5])
plt.show()
# Output value is Gaussian kernel multiplied by all positive samples
lhood_h_est_kern = np.zeros(len(x))
for xind, xval in enumerate(x):
lhood_h_est_kern[xind] = sum(gaussian(x_h,xval,kern_width,1))/len(x_h)
#lhood_h_est_kern[xind] = sum(stats.norm.pdf(x_h, xval, kernel_width))
plt.plot(x_h,np.zeros([N_h,1]),’co’, label=”hobit”)
plt.plot(x_e,np.zeros([N_e,1]),’mo’, label=”elf”)
plt.plot(x,lhood_h_gt,’c-‘, label=”hobit (GT)”)
plt.plot(x,lhood_h_est_kern,’c–‘, label=”hobit (est)”)
#plt.plot(x,lhood_e_est,’m–‘, label=”haltija (est)”)
plt.legend()
plt.xlabel(‘height [m]’)
#plt.axis([0.0,3.0,-1.1,+5])
plt.show()
#
#
kern_width = 0.8
[x, step_size] = np.linspace(0,3.0,70,retstep=True)
plt plot(x_h np zeros([N_h 1]),’co’ label “hobit”)
lec05_probability_distributions http://localhost:8888/nbconvert/html/lec05_probability_…
18 of 21 9/13/21, 16:17
Avg height hobits 1.11 (GT: 1.10)
Avg height elves 1.85 (GT: 1.90)
Height st. deviation hobits 0.296 (GT: 0.300)
Height st. deviation elves 0.428 (GT: 0.400)
for xind, xval in enumerate(x):
lhood_h_est_kern[xind] = sum(gaussian(x_h,xval,kern_width,1))/len(x_h)
#lhood_h_est_kern[xind] = sum(stats.norm.pdf(x_h, xval, kernel_width))
plt.plot(x_h,np.zeros([N_h,1]),’co’, label=”hobit”)
plt.plot(x_e,np.zeros([N_e,1]),’mo’, label=”elf”)
plt.plot(x,lhood_h_gt,’c-‘, label=”hobit (GT)”)
plt.plot(x,lhood_h_est_kern,’c–‘, label=”hobit (est)”)
#plt.plot(x,lhood_e_est,’m–‘, label=”haltija (est)”)
plt.legend()
plt.xlabel(‘height [m]’)
#plt.axis([0.0,3.0,-1.1,+5])
plt.show()
lec05_probability_distributions http://localhost:8888/nbconvert/html/lec05_probability_…
19 of 21 9/13/21, 16:17
lec05_probability_distributions http://localhost:8888/nbconvert/html/lec05_probability_…
20 of 21 9/13/21, 16:17
References
C.M. Bishop (2006): Pattern Recognition and Machine Learning, Chapter 1-2.
lec05_probability_distributions http://localhost:8888/nbconvert/html/lec05_probability_…
21 of 21 9/13/21, 16:17