From 92248b62e632c3474f41d19612d12a1a4b3d36e6 Mon Sep 17 00:00:00 2001 From: Zhenwen Dai Date: Tue, 3 Nov 2015 15:09:59 +0000 Subject: [PATCH] add sample_W to SSGPLVM --- GPy/models/ss_gplvm.py | 36 +++++++++++++++++++++++++++++++++++- 1 file changed, 35 insertions(+), 1 deletion(-) diff --git a/GPy/models/ss_gplvm.py b/GPy/models/ss_gplvm.py index 4ccf58ac..0123d629 100644 --- a/GPy/models/ss_gplvm.py +++ b/GPy/models/ss_gplvm.py @@ -190,4 +190,38 @@ class SSGPLVM(SparseGP_MPI): if self.kern.ARD: return self.kern.input_sensitivity() else: - return self.variational_prior.pi \ No newline at end of file + return self.variational_prior.pi + + def sample_W(self, nSamples, raw_samples=False): + """ + Sample the loading matrix if the kernel is linear. + """ + assert isinstance(self.kern, kern.Linear) + from ..util.linalg import pdinv + N, D = self.Y.shape + Q = self.X.shape[1] + noise_var = self.likelihood.variance.values + + # Draw samples for X + Xs = np.random.randn(*((nSamples,)+self.X.shape))*np.sqrt(self.X.variance.values)+self.X.mean.values + b = np.random.rand(*((nSamples,)+self.X.shape)) + Xs[b>self.X.gamma.values] = 0 + + invcov = (Xs[:,:,:,None]*Xs[:,:,None,:]).sum(1)/noise_var+np.eye(Q) + cov = np.array([pdinv(invcov[s_idx])[0] for s_idx in xrange(invcov.shape[0])]) + Ws = np.empty((nSamples, Q, D)) + tmp = (np.transpose(Xs, (0,2,1)).reshape(nSamples*Q,N).dot(self.Y)).reshape(nSamples,Q,D) + mean = (cov[:,:,:,None]*tmp[:,None,:,:]).sum(2)/noise_var + zeros = np.zeros((Q,)) + for s_idx in xrange(Xs.shape[0]): + Ws[s_idx] = (np.random.multivariate_normal(mean=zeros,cov=cov[s_idx],size=(D,))).T+mean[s_idx] + + if raw_samples: + return Ws + else: + return Ws.mean(0), Ws.std(0) + + + + + \ No newline at end of file