1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
|
import numpy as np
import seaborn
import matplotlib.pyplot as plt
from math import sqrt
from scipy.stats import multivariate_normal
seaborn.set_style("white")
x = np.array([-1.87, -1.76, -1.67, -1.22, -0.07, 0.11, 0.67, 1.60, 2.22, 2.51])
y = np.array([0.06, 1.67, 0.54, -1.45, -0.18, -0.67, 0.92, 2.95, 5.13, 5.18])
D = 20
eps = 1
base = np.linspace(-2, 3, num=D)
def map_features(x, base=base):
x = x.reshape(len(x), 1)
phi = np.exp(-(x - base) ** 2 / (2 * eps ** 2))
phi = np.concatenate((phi, np.ones((phi.shape[0], 1))), 1)
return phi
def draw_prior(x):
phi = map_features(x)
D = phi.shape[1]
w = np.random.normal(0, 1, D)
return np.dot(phi, w)
def posterior():
phi = map_features(x)
D = phi.shape[1]
invsigma = np.identity(D) + np.dot(phi.T, phi)
sigma = np.linalg.inv(invsigma)
mu = np.dot(sigma, np.dot(phi.T, y))
return mu, sigma
def plot_posteriors():
plt.figure(figsize=(9, 6))
mu, sigma = posterior()
px = np.linspace(-2, 3, num=200)
phi = map_features(px)
plt.plot(x, y, "o", label="data")
for i in xrange(5):
w = np.random.multivariate_normal(mu, sigma)
plt.plot(px, np.dot(phi, w), label="post. " + str(i))
plt.legend(loc="lower right")
plt.savefig("posterior.pdf")
def plot_tube():
fig = plt.figure(figsize=(9, 6))
mu, sigma = posterior()
px = np.linspace(-2, 3, num=200)
phi = map_features(px)
plt.plot(px, np.dot(phi, mu), label="post. pred. mean")
ax = fig.get_axes()[0]
inf = [np.dot(row, mu) - 1.96 * sqrt(1 + np.dot(row, np.dot(sigma, row)))
for row in phi]
sup = [np.dot(row, mu) + 1.96 * sqrt(1 + np.dot(row, np.dot(sigma, row)))
for row in phi]
ax.fill_between(px, inf, sup, alpha=0.2)
plt.plot(x, y, "o", label="data")
plt.legend(loc="lower right")
plt.savefig("predictive.pdf")
def plot_priors():
plt.figure(figsize=(9, 6))
px = np.linspace(-2, 3, num=200)
for i in xrange(5):
plt.plot(px, draw_prior(px), label="prior " + str(i))
plt.legend(loc="lower right")
plt.savefig("prior.pdf")
def marginal_likelihood():
eps = 1
ml = []
px = np.linspace(1, 50, num=10)
for D in px:
base = np.linspace(-2, 3, num=D)
phi = map_features(x, base)
z = np.zeros(phi.shape[0])
rv = multivariate_normal(mean=z, cov=np.identity(phi.shape[0])
+ np.dot(phi, phi.T))
ml.append(rv.pdf(y))
plt.figure(figsize=(9, 6))
plt.bar(px, ml)
plt.savefig("ml.pdf")
if __name__ == "__main__":
plot_priors()
plot_posteriors()
plot_tube()
marginal_likelihood()
#plt.plot(x, y, "o")
#plt.show()
|
