Allowing EP in BGPLVM and MRD

This commit is contained in:
Ricardo 2013-05-20 16:22:43 +01:00
parent afcb30dfbe
commit c3ab4b7979
3 changed files with 46 additions and 38 deletions

View file

@ -26,11 +26,11 @@ class Bayesian_GPLVM(sparse_GP, GPLVM):
:type init: 'PCA'|'random' :type init: 'PCA'|'random'
""" """
def __init__(self, Y, Q, X=None, X_variance=None, init='PCA', M=10, def __init__(self, likelihood, Q, X=None, X_variance=None, init='PCA', M=10,
Z=None, kernel=None, oldpsave=10, _debug=False, Z=None, kernel=None, oldpsave=10, _debug=False,
**kwargs): **kwargs):
if X == None: if X == None:
X = self.initialise_latent(init, Q, Y) X = self.initialise_latent(init, Q, likelihood.Y)
if X_variance is None: if X_variance is None:
X_variance = np.clip((np.ones_like(X) * 0.5) + .01 * np.random.randn(*X.shape), 0.001, 1) X_variance = np.clip((np.ones_like(X) * 0.5) + .01 * np.random.randn(*X.shape), 0.001, 1)
@ -56,7 +56,7 @@ class Bayesian_GPLVM(sparse_GP, GPLVM):
self._savedpsiKmm = [] self._savedpsiKmm = []
self._savedABCD = [] self._savedABCD = []
sparse_GP.__init__(self, X, Gaussian(Y), kernel, Z=Z, X_variance=X_variance, **kwargs) sparse_GP.__init__(self, X, likelihood, kernel, Z=Z, X_variance=X_variance, **kwargs)
@property @property
def oldps(self): def oldps(self):

View file

@ -15,11 +15,11 @@ import pylab
class MRD(model): class MRD(model):
""" """
Do MRD on given Datasets in Ylist. Do MRD on given Datasets in Ylist.
All Ys in Ylist are in [N x Dn], where Dn can be different per Yn, All Ys in likelihood_list are in [N x Dn], where Dn can be different per Yn,
N must be shared across datasets though. N must be shared across datasets though.
:param Ylist...: observed datasets :param likelihood_list...: likelihoods of observed datasets
:type Ylist: [np.ndarray] :type likelihood_list: [GPy.likelihood]
:param names: names for different gplvm models :param names: names for different gplvm models
:type names: [str] :type names: [str]
:param Q: latent dimensionality (will raise :param Q: latent dimensionality (will raise
@ -41,8 +41,9 @@ class MRD(model):
:param kernel: :param kernel:
kernel to use kernel to use
""" """
#TODO allow different kernels for different outputs
def __init__(self, *Ylist, **kwargs): #def __init__(self, *Ylist, **kwargs):
def __init__(self, *likelihood_list, **kwargs):
if kwargs.has_key("_debug"): if kwargs.has_key("_debug"):
self._debug = kwargs['_debug'] self._debug = kwargs['_debug']
del kwargs['_debug'] del kwargs['_debug']
@ -52,7 +53,7 @@ class MRD(model):
self.names = kwargs['names'] self.names = kwargs['names']
del kwargs['names'] del kwargs['names']
else: else:
self.names = ["{}".format(i + 1) for i in range(len(Ylist))] self.names = ["{}".format(i + 1) for i in range(len(likelihood_list))]
if kwargs.has_key('kernel'): if kwargs.has_key('kernel'):
kernel = kwargs['kernel'] kernel = kwargs['kernel']
k = lambda: kernel.copy() k = lambda: kernel.copy()
@ -80,9 +81,10 @@ class MRD(model):
self.M = 10 self.M = 10
self._init = True self._init = True
X = self._init_X(initx, Ylist) X = self._init_X(initx, likelihood_list)
Z = self._init_Z(initz, X) Z = self._init_Z(initz, X)
self.bgplvms = [Bayesian_GPLVM(Y, kernel=k(), X=X, Z=Z, M=self.M, **kwargs) for Y in Ylist] self.bgplvms = [Bayesian_GPLVM(Y, kernel=k(), X=X, Z=Z, M=self.M, **kwargs) for Y in likelihood_list]
del self._init del self._init
self.gref = self.bgplvms[0] self.gref = self.bgplvms[0]
@ -126,11 +128,11 @@ class MRD(model):
if not self._init: if not self._init:
raise AttributeError("bgplvm list not initialized") raise AttributeError("bgplvm list not initialized")
@property @property
def Ylist(self): def likelihood_list(self):
return [g.likelihood.Y for g in self.bgplvms] return [g.likelihood.Y for g in self.bgplvms]
@Ylist.setter @likelihood_list.setter
def Ylist(self, Ylist): def likelihood_list(self, likelihood_list):
for g, Y in itertools.izip(self.bgplvms, Ylist): for g, Y in itertools.izip(self.bgplvms, likelihood_list):
g.likelihood.Y = Y g.likelihood.Y = Y
@property @property
@ -152,7 +154,7 @@ class MRD(model):
def randomize(self, initx='concat', initz='permute', *args, **kw): def randomize(self, initx='concat', initz='permute', *args, **kw):
super(MRD, self).randomize(*args, **kw) super(MRD, self).randomize(*args, **kw)
self._init_X(initx, self.Ylist) self._init_X(initx, self.likelihood_list)
self._init_Z(initz, self.X) self._init_Z(initz, self.X)
def _get_param_names(self): def _get_param_names(self):
@ -225,6 +227,10 @@ class MRD(model):
# g._computations() # g._computations()
def update_likelihood_approximation(self):#TODO: object oriented vs script base
for bgplvm in self.bgplvms:
bgplvm.update_likelihood_approximation()
def log_likelihood(self): def log_likelihood(self):
ll = -self.gref.KL_divergence() ll = -self.gref.KL_divergence()
for g in self.bgplvms: for g in self.bgplvms:
@ -246,17 +252,18 @@ class MRD(model):
partial=g.partial_for_likelihood)]) \ partial=g.partial_for_likelihood)]) \
for g in self.bgplvms]))) for g in self.bgplvms])))
def _init_X(self, init='PCA', Ylist=None): def _init_X(self, init='PCA', likelihood_list=None):
if Ylist is None: if likelihood_list is None:
Ylist = self.Ylist likelihood_list = self.likelihood_list
if init in "PCA_single": if init in "PCA_single":
X = numpy.zeros((Ylist[0].shape[0], self.Q)) X = numpy.zeros((likelihood_list[0].Y.shape[0], self.Q))
for qs, Y in itertools.izip(numpy.array_split(numpy.arange(self.Q), len(Ylist)), Ylist): for qs, Y in itertools.izip(numpy.array_split(numpy.arange(self.Q), len(likelihood_list)), likelihood_list):
X[:, qs] = PCA(Y, len(qs))[0] X[:, qs] = PCA(Y.Y, len(qs))[0]
elif init in "PCA_concat": elif init in "PCA_concat":
X = PCA(numpy.hstack(Ylist), self.Q)[0] X = PCA(numpy.hstack([l.Y for l in likelihood_list]), self.Q)[0]
#X = PCA(numpy.hstack(likelihood_list), self.Q)[0]
else: # init == 'random': else: # init == 'random':
X = numpy.random.randn(Ylist[0].shape[0], self.Q) X = numpy.random.randn(likelihood_list[0].Y.shape[0], self.Q)
self.X = X self.X = X
return X return X
@ -294,8 +301,8 @@ class MRD(model):
fig = self._handle_plotting(fig_num, axes, lambda i, g, ax: ax.imshow(g.X)) fig = self._handle_plotting(fig_num, axes, lambda i, g, ax: ax.imshow(g.X))
return fig return fig
def plot_predict(self, fig_num="MRD Predictions", axes=None): def plot_predict(self, fig_num="MRD Predictions", axes=None, **kwargs):
fig = self._handle_plotting(fig_num, axes, lambda i, g, ax: ax.imshow(g.predict(g.X)[0])) fig = self._handle_plotting(fig_num, axes, lambda i, g, ax: ax.imshow(g.predict(g.X)[0],**kwargs))
return fig return fig
def plot_scales(self, fig_num="MRD Scales", axes=None, *args, **kwargs): def plot_scales(self, fig_num="MRD Scales", axes=None, *args, **kwargs):

View file

@ -8,6 +8,7 @@ from ..util.plot import gpplot
from .. import kern from .. import kern
from GP import GP from GP import GP
from scipy import linalg from scipy import linalg
from ..likelihoods import Gaussian
class sparse_GP(GP): class sparse_GP(GP):
""" """
@ -172,19 +173,19 @@ class sparse_GP(GP):
For a Gaussian likelihood, no iteration is required: For a Gaussian likelihood, no iteration is required:
this function does nothing this function does nothing
""" """
if self.has_uncertain_inputs: if not isinstance(self.likelihood,Gaussian): #Updates not needed for Gaussian likelihood
self.likelihood.restart() #TODO check consistency with pseudo_EP
Lmi = chol_inv(self.Lm) if self.has_uncertain_inputs:
Kmmi = tdot(Lmi.T) Lmi = chol_inv(self.Lm)
diag_tr_psi2Kmmi = np.array([np.trace(psi2_Kmmi) for psi2_Kmmi in np.dot(self.psi2,Kmmi)]) Kmmi = tdot(Lmi.T)
diag_tr_psi2Kmmi = np.array([np.trace(psi2_Kmmi) for psi2_Kmmi in np.dot(self.psi2,Kmmi)])
self.likelihood.fit_FITC(self.Kmm,self.psi1,diag_tr_psi2Kmmi) #This uses the fit_FITC code, but does not perfomr a FITC-EP.#TODO solve potential confusion
#raise NotImplementedError, "EP approximation not implemented for uncertain inputs"
else:
self.likelihood.fit_DTC(self.Kmm, self.psi1)
# self.likelihood.fit_FITC(self.Kmm,self.psi1,self.psi0)
self._set_params(self._get_params()) # update the GP
self.likelihood.fit_FITC(self.Kmm,self.psi1,diag_tr_psi2Kmmi) #This uses the fit_FITC code, but does not perfomr a FITC-EP.#TODO solve potential confusion
#raise NotImplementedError, "EP approximation not implemented for uncertain inputs"
else:
self.likelihood.fit_DTC(self.Kmm, self.psi1)
# self.likelihood.fit_FITC(self.Kmm,self.psi1,self.psi0)
self._set_params(self._get_params()) # update the GP
def _log_likelihood_gradients(self): def _log_likelihood_gradients(self):
return np.hstack((self.dL_dZ().flatten(), self.dL_dtheta(), self.likelihood._gradients(partial=self.partial_for_likelihood))) return np.hstack((self.dL_dZ().flatten(), self.dL_dtheta(), self.likelihood._gradients(partial=self.partial_for_likelihood)))