In [1]:
import numpy as np
import scipy.stats
import matplotlib.pylab as plt
import os, sys
sys.path.insert(0, "../")
import geepee.aep_models as aep
import geepee.ep_models as ep
%matplotlib inline
np.random.seed(42)
import pdb
In [2]:
# We first define several utility functions
def lincos_true(x):
fx = -0.5*x + 5*np.cos(0.5*x)
return fx
def lincos(T, process_noise, obs_noise, xprev=None):
if xprev is None:
xprev = np.random.randn()
y = np.zeros([T, ])
x = np.zeros([T, ])
xtrue = np.zeros([T, ])
for t in range(T):
fx = -0.5*xprev + 5*np.cos(0.5*xprev)
xtrue[t] = fx
x[t] = fx + np.sqrt(process_noise)*np.random.randn()
xprev = x[t]
y[t] = 0.5*x[t] + np.sqrt(obs_noise)*np.random.randn()
return xtrue, x, y
def plot_latent(model, y, plot_title=''):
# make prediction on some test inputs
N_test = 500
x_test = np.linspace(-10, 8, N_test)
x_test = np.reshape(x_test, [N_test, 1])
if isinstance(model, aep.SGPSSM_Linear):
zu = model.dyn_layer.zu
else:
zu = model.sgp_layer.zu
mu, vu = model.predict_f(zu)
mf, vf = model.predict_f(x_test)
# plot function
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(x_test[:,0], lincos_true(x_test[:,0]), '-', color='k')
ax.plot(zu, mu, 'ob')
ax.plot(x_test[:,0], mf[:,0], '-', color='b', label='f, alpha=%.3f' % alpha)
ax.fill_between(
x_test[:,0],
mf[:,0] + 2*np.sqrt(vf[:,0]),
mf[:,0] - 2*np.sqrt(vf[:,0]),
alpha=0.2, edgecolor='b', facecolor='b')
ax.plot(
y[0:model.N-1],
y[1:model.N],
'r+', alpha=0.5)
ax.set_xlabel(r'$x_{t-1}$')
ax.set_ylabel(r'$x_{t}$')
ax.set_xlim([-10, 8])
plt.title(plot_title)
plt.savefig('/tmp/lincos_'+plot_title+'.pdf')
In [ ]:
# generate a dataset from the lincos function above
T = 200
process_noise = 0.3
obs_noise = 0.2
(xtrue, x, y) = lincos(T, process_noise, obs_noise)
y_train = np.reshape(y, [y.shape[0], 1])
# init hypers
Dlatent = 1
Dobs = 1
M = 30
C = 0.5*np.ones((1, 1))
R = np.ones(1)*np.log(obs_noise)/2
lls = np.reshape(np.log(2), [Dlatent, ])
lsf = np.reshape(np.log(2), [1, ])
zu = np.linspace(-8, 6, M)
zu = np.reshape(zu, [M, 1])
lsn = np.log(0.01)/2
params = {'ls': lls, 'sf': lsf, 'sn': lsn, 'R': R, 'C': C, 'zu': zu}
In [ ]:
alphas = [0.001, 0.05, 0.2, 0.5, 1.0]
for alpha in alphas:
print 'alpha = %.3f' % alpha
# create AEP model
model_aep = aep.SGPSSM_Linear(y_train, Dlatent, M,
lik='Gaussian', prior_mean=0, prior_var=1000)
hypers = model_aep.init_hypers_old(y_train)
for key in params.keys():
hypers[key] = params[key]
model_aep.update_hypers(hypers)
# optimise
model_aep.set_fixed_params(['C', 'R'])
model_aep.optimise(method='L-BFGS-B', alpha=alpha, maxiter=10000, reinit_hypers=False)
opt_hypers = model_aep.get_hypers()
plot_latent(model_aep, x, 'AEP %.3f'%alpha)
# create EP model
model_ep = ep.SGPSSM(y_train, Dlatent, M,
lik='Gaussian', prior_mean=0, prior_var=1000)
# init EP model using the AEP solution
model_ep.update_hypers(opt_hypers)
# run EP
if alpha == 1.0:
decay = 0.999
parallel = True
no_epochs = 200
elif alpha == 0.001 or alpha == 0.05 or alpha ==0.2:
decay = 0.9
parallel = True
no_epochs = 1000
else:
decay = 0.99
parallel = True
no_epochs = 500
model_ep.inference(no_epochs=no_epochs, alpha=alpha, parallel=parallel, decay=decay)
plot_latent(model_ep, x, 'PEP %.3f'%alpha)
# create EP model
model_ep = ep.SGPSSM(y_train, Dlatent, M,
lik='Gaussian', prior_mean=0, prior_var=1000)
# init EP model using the AEP solution
model_ep.update_hypers(opt_hypers)
aep_sgp_layer = model_aep.dyn_layer
Nm1 = aep_sgp_layer.N
model_ep.sgp_layer.t1 = 1.0/Nm1 * np.tile(
aep_sgp_layer.theta_2[np.newaxis, :, :], [Nm1, 1, 1])
model_ep.sgp_layer.t2 = 1.0/Nm1 * np.tile(
aep_sgp_layer.theta_1[np.newaxis, :, :, :], [Nm1, 1, 1, 1])
model_ep.x_prev_1 = np.copy(model_aep.x_factor_1)
model_ep.x_prev_2 = np.copy(model_aep.x_factor_2)
model_ep.x_next_1 = np.copy(model_aep.x_factor_1)
model_ep.x_next_2 = np.copy(model_aep.x_factor_2)
model_ep.x_up_1 = np.copy(model_aep.x_factor_1)
model_ep.x_up_2 = np.copy(model_aep.x_factor_2)
model_ep.x_prev_1[0, :] = 0
model_ep.x_prev_2[0, :] = 0
model_ep.x_next_1[-1, :] = 0
model_ep.x_next_2[-1, :] = 0
# run EP
if alpha == 1.0:
decay = 0.999
parallel = True
no_epochs = 200
elif alpha == 0.001 or alpha == 0.05 or alpha == 0.2:
decay = 0.9
parallel = True
no_epochs = 1000
else:
decay = 0.99
parallel = True
no_epochs = 500
model_ep.inference(no_epochs=no_epochs, alpha=alpha, parallel=parallel, decay=decay)
plot_latent(model_ep, x, 'PEP (AEP init) %.3f'%alpha)
In [ ]: