In [2]:
import numpy as np
from numpy import transpose as tr
import scipy.stats as stats
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import matplotlib.colors as col
import matplotlib.patches as patches
import sys
import copy
sys.path.append('../src')
np.set_printoptions(precision=4, suppress=True)
%matplotlib inline
%load_ext autoreload
%autoreload 2
import model
import data
In [47]:
def plot_fit(costs):
fig, ax = plt.subplots(figsize=(10, 4))
ax.plot(costs, linewidth=2, marker='s',markersize=5, markerfacecolor='red')
ax.grid()
ax.set_xlabel('Iteration')
ax.set_ylabel('MSE')
def angle(v):
a = np.arccos(v[0] / np.linalg.norm(v)) / np.pi * 180
if v[1] >= 0:
return a
else:
return 360 - a
# https://onlinecourses.science.psu.edu/stat505/node/36
def ellipse_shape(cov, alpha=0.05):
l, e = np.linalg.eig(cov)
df = len(l)
widths = [np.sqrt(np.maximum(0.1, l[i]*stats.chi2.ppf(1-alpha, df=df))) for i in range(df)]
return [widths, angle(e[:,0])]
def lm_to_ellipse(lm):
widths, angle = ellipse_shape(lm.l.mean.dot(tr(lm.l.mean)))
return patches.Ellipse(lm.m.mean, widths[0], widths[1], angle)
def plot_grid(n, ncols=4, size=(5, 5)):
nrows = int(np.ceil(n/float(ncols)))
fig, ax = plt.subplots(nrows=nrows, ncols=ncols, figsize=(size[0]*ncols, size[1]*nrows))
ax = ax.ravel()
return [fig, ax]
def plot_predictions(models, xlim=None, ylim=None, ellipse=True):
n = len(models)
fig, ax = plot_grid(n)
for i in range(n):
model = models[i]
if ellipse:
for s in range(model.S):
e = ax[i].add_artist(lm_to_ellipse(model.q_lm[s]))
e.set_edgecolor('orange')
e.set_facecolor('orange')
e.set_alpha(max(0.1, model.q_pi.expectation()[s]))
e.set_linewidth(1)
e.set_zorder(0)
ax[i].scatter(model.y[0,:], model.y[1,:], color='blue')
yp = model.predict_y()
ax[i].scatter(yp[0,:], yp[1,:], color='red')
ax[i].set_title('Iteration {:d}'.format(i))
ax[i].set_xlim(xlim)
ax[i].set_ylim(ylim)
def plot_lambda(ax, lm, norm=None, cmap=cm.jet):
X = lm.l.mean
if norm is None:
norm = col.Normalize(np.min(X), np.max(X))
im = ax.imshow(X, interpolation='nearest', cmap=cmap, norm=norm, aspect='auto')
ax.set_yticks(range(X.shape[0]))
ax.set_xticks(range(X.shape[1]))
return im
def plot_lambdas(models, width=3, height=3):
M = len(models)
S = models[0].S
l_min = 0.0
l_max = 0.0
for i in range(M):
for j in range(S):
l_min = np.minimum(l_min, np.min(models[i].q_lm[j].l.mean))
l_max = np.maximum(l_max, np.max(models[i].q_lm[j].l.mean))
norm = col.Normalize(l_min, l_max)
fig, ax = plt.subplots(nrows=M, ncols=S, figsize=(width*S, height*M))
ax = ax.reshape(M, S)
for i in range(M):
for j in range(S):
lm = models[i].q_lm[j]
im = plot_lambda(ax[i, j], lm, norm)
fig.colorbar(im, orientation='horizontal')
return [fig, ax]
In [51]:
np.random.seed(0)
mr = model.Model(model.Hyper(2, 2, 1), 20)
mr.q_lm[0].l.mean = np.array([[10.0, -2.0],
[10.0, 2.0]])
mr.q_lm[0].m.mean = np.array([5.0, 5.0])
mr.y = mr.predict_y()
m0 = model.Model(model.Hyper(2, 2, 1), mr.y)
m0.init_rnd()
m = copy.deepcopy(m0)
def update(model):
model.update_lm()
model.update_x()
model.update_nu()
costs, models = m.infer(update=update)
plot_fit(costs)
plot_predictions(models[:8], xlim=(-2, 8), ylim=(-2, 8), ellipse=True)
In [130]:
np.random.seed(0)
mr = model.Model(model.Hyper(2, 2, 4), 80)
mr.q_lm[0].m.mean = np.array([-10.0, -10.0])
mr.q_lm[1].m.mean = np.array([10.0, -10.0])
mr.q_lm[2].m.mean = np.array([-10.0, 10.0])
mr.q_lm[3].m.mean = np.array([10.0, 10.0])
mr.q_s.s[0, 0:20] = 1.0
mr.q_s.s[1, 20:40] = 1.0
mr.q_s.s[2, 40:60] = 1.0
mr.q_s.s[3, 60:80] = 1.0
mr.y = mr.predict_y()
m0 = model.Model(model.Hyper(mr.P, mr.Q, 1), mr.y)
m0.init_rnd()
m = copy.deepcopy(m0)
def update(m, it):
m.update_lm()
m.update_x()
m.update_nu()
#m.update_s()
#m.update_pi()
costs, models = m.infer(update=update, eps=1e-5)
plot_fit(costs)
plot_predictions(models[:8])
In [138]:
np.random.seed(0)
mr = model.Model(model.Hyper(2, 2, 4), 80)
mr.q_lm[0].m.mean = np.array([-10.0, -10.0])
mr.q_lm[1].m.mean = np.array([10.0, -10.0])
mr.q_lm[2].m.mean = np.array([-10.0, 10.0])
mr.q_lm[3].m.mean = np.array([10.0, 10.0])
mr.q_s.s[0, 0:20] = 1.0
mr.q_s.s[1, 20:40] = 1.0
mr.q_s.s[2, 40:60] = 1.0
mr.q_s.s[3, 60:80] = 1.0
mr.y = mr.predict_y()
m0 = model.Model(model.Hyper(mr.P, mr.Q, mr.S), mr.y)
m0.init_rnd()
m0.init_complex()
m = copy.deepcopy(m0)
def update(m, it):
m.update_lm()
m.update_x()
m.update_nu()
m.update_s()
m.update_pi()
costs, models = m.infer(update=update, eps=1e-5)
plot_fit(costs)
plot_predictions(models[:8])
In [132]:
np.random.seed(0)
mr = model.Model(model.Hyper(2, 2, 4), 80)
for s in range(mr.S):
mr.q_lm[s].l.mean = np.random.normal(0.0, 1e-3, mr.P*mr.Q).reshape(mr.P, mr.Q)
#mr.q_lm[s].m.mean.fill(-5.0)
mr.q_lm[0].l.mean[:, 0] = [0.0, 10.0]
mr.q_lm[0].m.mean = np.array([0.0, 0.0])
mr.q_lm[1].l.mean[:, 0] = [4.0, 0.0]
mr.q_lm[1].m.mean = np.array([1.5, 2.0])
mr.q_lm[2].l.mean[:, 0] = [4.0, 0.0]
mr.q_lm[2].m.mean = np.array([1.2, 0.0])
mr.q_lm[3].l.mean[:, 0] = [4.0, 0.0]
mr.q_lm[3].m.mean = np.array([1.5, -2.0])
mr.q_s.s.fill(0.0)
mr.q_s.s[0, 0:20] = 1.0
mr.q_s.s[1, 20:40] = 1.0
mr.q_s.s[2, 40:60] = 1.0
mr.q_s.s[3, 60:80] = 1.0
mr.y = mr.predict_y()
m0 = model.Model(model.Hyper(mr.P, mr.Q, mr.S), mr.y)
#m0.init_comple ytre3w2q=653x(hard=False)
m0.init_rnd()
m = copy.deepcopy(m0)
def update(m, it):
m.update_lm()
m.update_x()
m.update_s()
m.update_nu()
m.update_pi()
costs, models = m.infer(update=update)
plot_predictions(models)