apunkt/GPy
1
0
Fork 0
mirror of https://github.com/SheffieldML/GPy.git synced 2026-05-13 22:12:38 +02:00
Code Issues Projects Releases Packages Wiki Activity Actions

Merge branch 'devel' of github.com:SheffieldML/GPy into devel

Browse source
This commit is contained in:
Nicolo Fusi 2013-06-17 15:56:21 +01:00
commit bbca026a21
7 changed files with 70 additions and 20 deletions
Download patch file Download diff file Expand all files Collapse all files

33
GPy/core/sparse_gp.py
View file

@ -63,6 +63,7 @@ class SparseGP(GPBase):
def _computations(self):
# factor Kmm
self.Lm = jitchol(self.Kmm)
@ -89,17 +90,18 @@ class SparseGP(GPBase):
self.B = np.eye(self.num_inducing) + self.A
self.LB = jitchol(self.B)
# TODO: make a switch for either first compute psi1V, or VV.T
self.psi1V = np.dot(self.psi1.T, self.likelihood.V)
#VVT_factor is a matrix such that tdot(VVT_factor) = VVT...this is for efficiency!
self.psi1Vf = np.dot(self.psi1.T, self.likelihood.VVT_factor)
# back substutue C into psi1V
tmp, info1 = dtrtrs(self.Lm, np.asfortranarray(self.psi1V), lower=1, trans=0)
self._LBi_Lmi_psi1V, _ = dtrtrs(self.LB, np.asfortranarray(tmp), lower=1, trans=0)
# back substutue C into psi1Vf
tmp, info1 = dtrtrs(self.Lm, np.asfortranarray(self.psi1Vf), lower=1, trans=0)
self._LBi_Lmi_psi1Vf, _ = dtrtrs(self.LB, np.asfortranarray(tmp), lower=1, trans=0)
tmp, info2 = dpotrs(self.LB, tmp, lower=1)
self.Cpsi1V, info3 = dtrtrs(self.Lm, tmp, lower=1, trans=1)
self.Cpsi1Vf, info3 = dtrtrs(self.Lm, tmp, lower=1, trans=1)
# Compute dL_dKmm
tmp = tdot(self._LBi_Lmi_psi1V)
tmp = tdot(self._LBi_Lmi_psi1Vf)
self.data_fit = np.trace(tmp)
self.DBi_plus_BiPBi = backsub_both_sides(self.LB, self.output_dim * np.eye(self.num_inducing) + tmp)
tmp = -0.5 * self.DBi_plus_BiPBi
tmp += -0.5 * self.B * self.output_dim
@ -108,7 +110,7 @@ class SparseGP(GPBase):
# Compute dL_dpsi # FIXME: this is untested for the heterscedastic + uncertain inputs case
self.dL_dpsi0 = -0.5 * self.output_dim * (self.likelihood.precision * np.ones([self.num_data, 1])).flatten()
self.dL_dpsi1 = np.dot(self.Cpsi1V, self.likelihood.V.T).T
self.dL_dpsi1 = np.dot(self.likelihood.VVT_factor, self.Cpsi1Vf.T)
dL_dpsi2_beta = 0.5 * backsub_both_sides(self.Lm, self.output_dim * np.eye(self.num_inducing) - self.DBi_plus_BiPBi)
if self.likelihood.is_heteroscedastic:
@ -138,18 +140,18 @@ class SparseGP(GPBase):
# likelihood is not heterscedatic
self.partial_for_likelihood = -0.5 * self.num_data * self.output_dim * self.likelihood.precision + 0.5 * self.likelihood.trYYT * self.likelihood.precision ** 2
self.partial_for_likelihood += 0.5 * self.output_dim * (self.psi0.sum() * self.likelihood.precision ** 2 - np.trace(self.A) * self.likelihood.precision)
self.partial_for_likelihood += self.likelihood.precision * (0.5 * np.sum(self.A * self.DBi_plus_BiPBi) - np.sum(np.square(self._LBi_Lmi_psi1V)))
self.partial_for_likelihood += self.likelihood.precision * (0.5 * np.sum(self.A * self.DBi_plus_BiPBi) - self.data_fit)
def log_likelihood(self):
""" Compute the (lower bound on the) log marginal likelihood """
if self.likelihood.is_heteroscedastic:
A = -0.5 * self.num_data * self.output_dim * np.log(2.*np.pi) + 0.5 * np.sum(np.log(self.likelihood.precision)) - 0.5 * np.sum(self.likelihood.V * self.likelihood.Y)
A = -0.5 * self.num_data * self.output_dim * np.log(2.*np.pi) + 0.5 * np.sum(np.log(self.likelihood.precision)) - 0.5 * np.sum(self.likelihood.V*self.likelihood.Y)
B = -0.5 * self.output_dim * (np.sum(self.likelihood.precision.flatten() * self.psi0) - np.trace(self.A))
else:
A = -0.5 * self.num_data * self.output_dim * (np.log(2.*np.pi) - np.log(self.likelihood.precision)) - 0.5 * self.likelihood.precision * self.likelihood.trYYT
B = -0.5 * self.output_dim * (np.sum(self.likelihood.precision * self.psi0) - np.trace(self.A))
C = -self.output_dim * (np.sum(np.log(np.diag(self.LB)))) # + 0.5 * self.num_inducing * np.log(sf2))
D = 0.5 * np.sum(np.square(self._LBi_Lmi_psi1V))
D = 0.5 * self.data_fit
return A + B + C + D + self.likelihood.Z
def _set_params(self, p):
@ -158,6 +160,7 @@ class SparseGP(GPBase):
self.likelihood._set_params(p[self.Z.size + self.kern.num_params:])
self._compute_kernel_matrices()
self._computations()
self.Cpsi1V = None
def _get_params(self):
return np.hstack([self.Z.flatten(), self.kern._get_params_transformed(), self.likelihood._get_params()])
@ -224,6 +227,14 @@ class SparseGP(GPBase):
symmetrify(Bi)
Kmmi_LmiBLmi = backsub_both_sides(self.Lm, np.eye(self.num_inducing) - Bi)
if self.Cpsi1V is None:
psi1V = np.dot(self.psi1.T,self.likelihood.V)
tmp, _ = dtrtrs(self.Lm, np.asfortranarray(psi1V), lower=1, trans=0)
tmp, _ = dpotrs(self.LB, tmp, lower=1)
self.Cpsi1V, _ = dtrtrs(self.Lm, tmp, lower=1, trans=1)
if X_variance_new is None:
Kx = self.kern.K(self.Z, Xnew, which_parts=which_parts)
mu = np.dot(Kx.T, self.Cpsi1V)

39
GPy/examples/dimensionality_reduction.py
View file

@ -261,6 +261,7 @@ def bgplvm_simulation(optimize='scg',
k = kern.linear(Q, ARD=True) + kern.bias(Q, np.exp(-2)) + kern.white(Q, np.exp(-2)) # + kern.bias(Q)
m = BayesianGPLVM(Y, Q, init="PCA", num_inducing=num_inducing, kernel=k, _debug=True)
# m.constrain('variance|noise', logexp_clipped())
m.ensure_default_constraints()
m['noise'] = Y.var() / 100.
@ -327,28 +328,56 @@ def brendan_faces():
data_show = GPy.util.visualize.image_show(y[None, :], dimensions=(20, 28), transpose=True, invert=False, scale=False)
lvm_visualizer = GPy.util.visualize.lvm(m.X[0, :].copy(), m, data_show, ax)
raw_input('Press enter to finish')
lvm_visualizer.close()
return m
def stick_play(range=None, frame_rate=15):
data = GPy.util.datasets.stick()
# optimize
if range==None:
Y = data['Y'].copy()
else:
Y = data['Y'][range[0]:range[1], :].copy()
y = Y[0, :]
data_show = GPy.util.visualize.stick_show(y[None, :], connect=data['connect'])
GPy.util.visualize.data_play(Y, data_show, frame_rate)
return Y
def stick():
data = GPy.util.datasets.stick()
m = GPy.models.GPLVM(data['Y'], 2)
# optimize
m = GPy.models.GPLVM(data['Y'], 2)
m.ensure_default_constraints()
m.optimize(messages=1, max_f_eval=10000)
m._set_params(m._get_params())
plt.clf
ax = m.plot_latent()
y = m.likelihood.Y[0, :]
data_show = GPy.util.visualize.stick_show(y[None, :], connect=data['connect'])
lvm_visualizer = GPy.util.visualize.lvm(m.X[0, :].copy(), m, data_show, ax)
raw_input('Press enter to finish')
lvm_visualizer.close()
return m
def stick_bgplvm(model=None):
data = GPy.util.datasets.stick()
Q = 6
kernel = GPy.kern.rbf(Q, ARD=True) + GPy.kern.bias(Q, np.exp(-2)) + GPy.kern.white(Q, np.exp(-2))
m = BayesianGPLVM(data['Y'], Q, init="PCA", num_inducing=20,kernel=kernel)
# optimize
m.ensure_default_constraints()
m.optimize(messages=1, max_f_eval=3000,xtol=1e-300,ftol=1e-300)
m._set_params(m._get_params())
plt.clf, (latent_axes, sense_axes) = plt.subplots(1, 2)
plt.sca(latent_axes)
m.plot_latent()
y = m.likelihood.Y[0, :].copy()
data_show = GPy.util.visualize.stick_show(y[None, :], connect=data['connect'])
lvm_visualizer = GPy.util.visualize.lvm_dimselect(m.X[0, :].copy(), m, data_show, latent_axes=latent_axes, sense_axes=sense_axes)
raw_input('Press enter to finish')
return m
def cmu_mocap(subject='35', motion=['01'], in_place=True):
data = GPy.util.datasets.cmu_mocap(subject, motion)

2
GPy/kern/Brownian.py
View file

@ -21,7 +21,7 @@ class Brownian(Kernpart):
def __init__(self,input_dim,variance=1.):
self.input_dim = input_dim
assert self.input_dim==1, "Brownian motion in 1D only"
self.num_params = 1.
self.num_params = 1
self.name = 'Brownian'
self._set_params(np.array([variance]).flatten())

6
GPy/likelihoods/ep.py
View file

@ -34,6 +34,8 @@ class EP(likelihood):
self.Z = 0
self.YYT = None
self.V = self.precision * self.Y
self.VVT_factor = self.V
self.trYYT = 0.
def restart(self):
self.tau_tilde = np.zeros(self.N)
@ -44,6 +46,8 @@ class EP(likelihood):
self.Z = 0
self.YYT = None
self.V = self.precision * self.Y
self.VVT_factor = self.V
self.trYYT = 0.
def predictive_values(self,mu,var,full_cov):
if full_cov:
@ -71,6 +75,8 @@ class EP(likelihood):
self.covariance_matrix = np.diag(1./self.tau_tilde)
self.precision = self.tau_tilde[:,None]
self.V = self.precision * self.Y
self.VVT_factor = self.V
self.trYYT = np.trace(self.YYT)
def fit_full(self,K):
"""

5
GPy/likelihoods/gaussian.py
View file

@ -40,9 +40,11 @@ class Gaussian(likelihood):
if D > self.N:
self.YYT = np.dot(self.Y, self.Y.T)
self.trYYT = np.trace(self.YYT)
self.YYT_factor = jitchol(self.YYT)
else:
self.YYT = None
self.trYYT = np.sum(np.square(self.Y))
self.YYT_factor = self.Y
def _get_params(self):
return np.asarray(self._variance)
@ -53,12 +55,13 @@ class Gaussian(likelihood):
def _set_params(self, x):
x = np.float64(x)
if np.all(self._variance != x):
if x == 0.:
if x == 0.:#special case of zero noise
self.precision = np.inf
self.V = None
else:
self.precision = 1. / x
self.V = (self.precision) * self.Y
self.VVT_factor = self.precision * self.YYT_factor
self.covariance_matrix = np.eye(self.N) * x
self._variance = x

3
GPy/util/visualize.py
View file

@ -203,6 +203,7 @@ class lvm_dimselect(lvm):
self.sense_axes = sense_axes
self.labels = labels
lvm.__init__(self,vals,Model,data_visualize,latent_axes,sense_axes,latent_index)
self.show_sensitivities()
print "use left and right mouse butons to select dimensions"
@ -506,5 +507,5 @@ def data_play(Y, visualizer, frame_rate=30):
for y in Y:
visualizer.modify(y)
visualizer.modify(y[None, :])
time.sleep(1./float(frame_rate))

2
setup.py
View file

@ -5,7 +5,7 @@ import os
from setuptools import setup
# Version number
version = '0.4.3'
version = '0.4.5'
def read(fname):
return open(os.path.join(os.path.dirname(__file__), fname)).read()