From 169fd9b8d4852b29f77ac05535f1cb8a311a7008 Mon Sep 17 00:00:00 2001 From: Alan Saul Date: Fri, 14 Mar 2014 12:42:36 +0000 Subject: [PATCH 01/13] Stablised exp link_function and quadrature variances --- GPy/likelihoods/likelihood.py | 8 +++++++- GPy/likelihoods/link_functions.py | 11 +++++++---- GPy/likelihoods/poisson.py | 2 +- 3 files changed, 15 insertions(+), 6 deletions(-) diff --git a/GPy/likelihoods/likelihood.py b/GPy/likelihoods/likelihood.py index 3eafedb1..4b8881de 100644 --- a/GPy/likelihoods/likelihood.py +++ b/GPy/likelihoods/likelihood.py @@ -178,7 +178,13 @@ class Likelihood(Parameterized): #E( E(Y_star|f_star)**2 ) def int_pred_mean_sq(f,m,v,predictive_mean_sq): - return self.conditional_mean(f)**2*np.exp(-(0.5/v)*np.square(f - m)) + p = np.exp(-(0.5/v)*np.square(f - m)) + #If p is zero then conditional_mean**2 will overflow + if p < 1e-10: + return 0. + else: + return self.conditional_mean(f)**2*p + scaled_exp_exp2 = [quad(int_pred_mean_sq, -np.inf, np.inf,args=(mj,s2j,pm2j))[0] for mj,s2j,pm2j in zip(mu,variance,predictive_mean_sq)] exp_exp2 = np.array(scaled_exp_exp2)[:,None] / normalizer diff --git a/GPy/likelihoods/link_functions.py b/GPy/likelihoods/link_functions.py index 2a1bf147..942fe2f4 100644 --- a/GPy/likelihoods/link_functions.py +++ b/GPy/likelihoods/link_functions.py @@ -6,6 +6,9 @@ from scipy import stats import scipy as sp from GPy.util.univariate_Gaussian import std_norm_pdf,std_norm_cdf,inv_std_norm_cdf +_exp_lim_val = np.finfo(np.float64).max +_lim_val = np.log(_exp_lim_val) + class GPTransformation(object): """ Link function class for doing non-Gaussian likelihoods approximation @@ -92,16 +95,16 @@ class Log(GPTransformation): """ def transf(self,f): - return np.exp(f) + return np.exp(np.clip(f, -_lim_val, _lim_val)) def dtransf_df(self,f): - return np.exp(f) + return np.exp(np.clip(f, -_lim_val, _lim_val)) def d2transf_df2(self,f): - return np.exp(f) + return np.exp(np.clip(f, -_lim_val, _lim_val)) def d3transf_df3(self,f): - return np.exp(f) + return np.exp(np.clip(f, -_lim_val, _lim_val)) class Log_ex_1(GPTransformation): """ diff --git a/GPy/likelihoods/poisson.py b/GPy/likelihoods/poisson.py index 419514d1..c67a7e12 100644 --- a/GPy/likelihoods/poisson.py +++ b/GPy/likelihoods/poisson.py @@ -21,7 +21,7 @@ class Poisson(Likelihood): """ def __init__(self, gp_link=None): if gp_link is None: - gp_link = link_functions.Log_ex_1() + gp_link = link_functions.Log() super(Poisson, self).__init__(gp_link, name='Poisson') From 7cd75ccbf58c5f01441669503d5da82a5d1f33d4 Mon Sep 17 00:00:00 2001 From: Alan Saul Date: Fri, 14 Mar 2014 12:52:52 +0000 Subject: [PATCH 02/13] Stablised other quadrature (should speed things up also), added sampling ability to poisson --- GPy/likelihoods/likelihood.py | 14 ++++++++++++-- GPy/likelihoods/poisson.py | 2 +- 2 files changed, 13 insertions(+), 3 deletions(-) diff --git a/GPy/likelihoods/likelihood.py b/GPy/likelihoods/likelihood.py index 4b8881de..98a856fd 100644 --- a/GPy/likelihoods/likelihood.py +++ b/GPy/likelihoods/likelihood.py @@ -142,7 +142,12 @@ class Likelihood(Parameterized): """ #conditional_mean: the edpected value of y given some f, under this likelihood def int_mean(f,m,v): - return self.conditional_mean(f)*np.exp(-(0.5/v)*np.square(f - m)) + p = np.exp(-(0.5/v)*np.square(f - m)) + #If p is zero then conditional_mean will overflow + if p < 1e-10: + return 0. + else: + return self.conditional_mean(f)*p scaled_mean = [quad(int_mean, -np.inf, np.inf,args=(mj,s2j))[0] for mj,s2j in zip(mu,variance)] mean = np.array(scaled_mean)[:,None] / np.sqrt(2*np.pi*(variance)) @@ -165,7 +170,12 @@ class Likelihood(Parameterized): # E( V(Y_star|f_star) ) def int_var(f,m,v): - return self.conditional_variance(f)*np.exp(-(0.5/v)*np.square(f - m)) + p = np.exp(-(0.5/v)*np.square(f - m)) + #If p is zero then conditional_variance will overflow + if p < 1e-10: + return 0. + else: + return self.conditional_variance(f)*p scaled_exp_variance = [quad(int_var, -np.inf, np.inf,args=(mj,s2j))[0] for mj,s2j in zip(mu,variance)] exp_var = np.array(scaled_exp_variance)[:,None] / normalizer diff --git a/GPy/likelihoods/poisson.py b/GPy/likelihoods/poisson.py index c67a7e12..c0e2c81f 100644 --- a/GPy/likelihoods/poisson.py +++ b/GPy/likelihoods/poisson.py @@ -143,7 +143,7 @@ class Poisson(Likelihood): """ return self.gp_link.transf(gp) - def samples(self, gp): + def samples(self, gp, Y_metadata=None): """ Returns a set of samples of observations based on a given value of the latent variable. From a12bdafad7d03e57bd40818cbcb49d671205d7b7 Mon Sep 17 00:00:00 2001 From: James Hensman Date: Fri, 14 Mar 2014 14:28:34 +0000 Subject: [PATCH 03/13] added jitter to fitc --- GPy/inference/latent_function_inference/fitc.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/GPy/inference/latent_function_inference/fitc.py b/GPy/inference/latent_function_inference/fitc.py index 9294e25d..c4147d06 100644 --- a/GPy/inference/latent_function_inference/fitc.py +++ b/GPy/inference/latent_function_inference/fitc.py @@ -3,6 +3,7 @@ from posterior import Posterior from ...util.linalg import jitchol, tdot, dtrtrs, dpotri, pdinv +from ...util import diag import numpy as np log_2_pi = np.log(2*np.pi) @@ -14,8 +15,7 @@ class FITC(object): the posterior. """ - def __init__(self): - self.const_jitter = 1e-6 + const_jitter = 1e-6 def inference(self, kern, X, Z, likelihood, Y): @@ -33,6 +33,7 @@ class FITC(object): U = Knm #factor Kmm + diag.add(Kmm, self.const_jitter) Kmmi, L, Li, _ = pdinv(Kmm) #compute beta_star, the effective noise precision From 4757265b240924af09ac0032984c5c8d85f4d805 Mon Sep 17 00:00:00 2001 From: Alan Saul Date: Fri, 14 Mar 2014 15:23:32 +0000 Subject: [PATCH 04/13] Added a hack fix as suggested by max, zeroing any negative values (should really be numerically negative values on diagonal) --- GPy/kern/_src/independent_outputs.py | 3 ++- GPy/kern/_src/stationary.py | 28 +++++++++++++++------------- GPy/testing/kernel_tests.py | 19 +++++++++++++------ 3 files changed, 30 insertions(+), 20 deletions(-) diff --git a/GPy/kern/_src/independent_outputs.py b/GPy/kern/_src/independent_outputs.py index 1848bf6a..85cb95bc 100644 --- a/GPy/kern/_src/independent_outputs.py +++ b/GPy/kern/_src/independent_outputs.py @@ -8,7 +8,7 @@ import itertools def index_to_slices(index): """ - take a numpy array of integers (index) and return a nested list of slices such that the slices describe the start, stop points for each integer in the index. + take a numpy array of integers (index) and return a nested list of slices such that the slices describe the start, stop points for each integer in the index. e.g. >>> index = np.asarray([0,0,0,1,1,1,2,2,2]) @@ -40,6 +40,7 @@ class IndependentOutputs(CombinationKernel): The index of the functions is given by the last column in the input X the rest of the columns of X are passed to the underlying kernel for computation (in blocks). + Kern is wrapped with a slicer metaclass """ def __init__(self, kern, index_dim=-1, name='independ'): assert isinstance(index_dim, int), "IndependentOutputs kernel is only defined with one input dimension being the indeces" diff --git a/GPy/kern/_src/stationary.py b/GPy/kern/_src/stationary.py index df7ba058..a9e837a9 100644 --- a/GPy/kern/_src/stationary.py +++ b/GPy/kern/_src/stationary.py @@ -15,21 +15,21 @@ class Stationary(Kern): """ Stationary kernels (covariance functions). - Stationary covariance fucntion depend only on r, where r is defined as + Stationary covariance fucntion depend only on r, where r is defined as r = \sqrt{ \sum_{q=1}^Q (x_q - x'_q)^2 } - The covariance function k(x, x' can then be written k(r). + The covariance function k(x, x' can then be written k(r). In this implementation, r is scaled by the lengthscales parameter(s): - r = \sqrt{ \sum_{q=1}^Q \frac{(x_q - x'_q)^2}{\ell_q^2} }. - + r = \sqrt{ \sum_{q=1}^Q \frac{(x_q - x'_q)^2}{\ell_q^2} }. + By default, there's only one lengthscale: seaprate lengthscales for each - dimension can be enables by setting ARD=True. + dimension can be enables by setting ARD=True. To implement a stationary covariance function using this class, one need - only define the covariance function k(r), and it derivative. + only define the covariance function k(r), and it derivative. ... def K_of_r(self, r): @@ -37,10 +37,10 @@ class Stationary(Kern): def dK_dr(self, r): return bar - The lengthscale(s) and variance parameters are added to the structure automatically. - + The lengthscale(s) and variance parameters are added to the structure automatically. + """ - + def __init__(self, input_dim, variance, lengthscale, ARD, active_dims, name): super(Stationary, self).__init__(input_dim, active_dims, name) self.ARD = ARD @@ -57,7 +57,7 @@ class Stationary(Kern): if lengthscale.size != input_dim: lengthscale = np.ones(input_dim)*lengthscale else: - lengthscale = np.ones(self.input_dim) + lengthscale = np.ones(self.input_dim) self.lengthscale = Param('lengthscale', lengthscale, Logexp()) self.variance = Param('variance', variance, Logexp()) assert self.variance.size==1 @@ -95,7 +95,9 @@ class Stationary(Kern): #X2, = self._slice_X(X2) X1sq = np.sum(np.square(X),1) X2sq = np.sum(np.square(X2),1) - return np.sqrt(-2.*np.dot(X, X2.T) + (X1sq[:,None] + X2sq[None,:])) + r2 = -2.*np.dot(X, X2.T) + X1sq[:,None] + X2sq[None,:] + r2[r2<0] = 0. # A bit hacky + return np.sqrt(r2) @Cache_this(limit=5, ignore_args=()) def _scaled_dist(self, X, X2=None): @@ -133,7 +135,7 @@ class Stationary(Kern): if self.ARD: #rinv = self._inv_dis# this is rather high memory? Should we loop instead?t(X, X2) #d = X[:, None, :] - X2[None, :, :] - #x_xl3 = np.square(d) + #x_xl3 = np.square(d) #self.lengthscale.gradient = -((dL_dr*rinv)[:,:,None]*x_xl3).sum(0).sum(0)/self.lengthscale**3 tmp = dL_dr*self._inv_dist(X, X2) if X2 is None: X2 = X @@ -247,7 +249,7 @@ class Matern52(Stationary): .. math:: - k(r) = \sigma^2 (1 + \sqrt{5} r + \\frac53 r^2) \exp(- \sqrt{5} r) + k(r) = \sigma^2 (1 + \sqrt{5} r + \\frac53 r^2) \exp(- \sqrt{5} r) """ def __init__(self, input_dim, variance=1., lengthscale=None, ARD=False, active_dims=None, name='Mat52'): super(Matern52, self).__init__(input_dim, variance, lengthscale, ARD, active_dims, name) diff --git a/GPy/testing/kernel_tests.py b/GPy/testing/kernel_tests.py index b057f8ef..3b17a3d0 100644 --- a/GPy/testing/kernel_tests.py +++ b/GPy/testing/kernel_tests.py @@ -120,6 +120,8 @@ def check_kernel_gradient_functions(kern, X=None, X2=None, output_ind=None, verb if verbose: print("Checking covariance function is positive definite.") + #if isinstance(kern, GPy.kern.IndependentOutputs): + #import ipdb; ipdb.set_trace() # XXX BREAKPOINT result = Kern_check_model(kern, X=X).is_positive_semi_definite() if result and verbose: print("Check passed.") @@ -306,17 +308,22 @@ class KernelTestsNonContinuous(unittest.TestCase): D = self.D self.X = np.random.randn(N,D) self.X2 = np.random.randn(N1,D) - self.X_block = np.zeros((N+N1, D+D+1)) + #self.X_block = np.zeros((N+N1, D+D+1)) + #self.X_block[0:N, 0:D] = self.X + #self.X_block[N:N+N1, D:D+D] = self.X2 + #self.X_block[0:N, -1] = 0 + #self.X_block[N:N+N1, -1] = 1 + self.X_block = np.zeros((N+N1, D+1)) self.X_block[0:N, 0:D] = self.X - self.X_block[N:N+N1, D:D+D] = self.X2 - self.X_block[0:N, -1] = 1 - self.X_block[N:N+1, -1] = 2 + self.X_block[N:N+N1, 0:D] = self.X2 + self.X_block[0:N, -1] = 0 + self.X_block[N:N+N1, -1] = 1 self.X_block = self.X_block[self.X_block.argsort(0)[:, -1], :] - + def test_IndependentOutputs(self): k = GPy.kern.RBF(self.D) kern = GPy.kern.IndependentOutputs(k, -1) - self.assertTrue(check_kernel_gradient_functions(kern, X=self.X_block, X2=self.X_block, verbose=verbose)) + self.assertTrue(check_kernel_gradient_functions(kern, X=self.X_block, verbose=verbose)) if __name__ == "__main__": print "Running unit tests, please be (very) patient..." From bfa18ef6625770733f24d154f12505b5353e69a3 Mon Sep 17 00:00:00 2001 From: Alan Saul Date: Fri, 14 Mar 2014 15:34:00 +0000 Subject: [PATCH 05/13] Fixed Y_metadata bug --- GPy/core/gp.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/GPy/core/gp.py b/GPy/core/gp.py index 05e3e671..70b7d695 100644 --- a/GPy/core/gp.py +++ b/GPy/core/gp.py @@ -42,7 +42,10 @@ class GP(Model): assert Y.shape[0] == self.num_data _, self.output_dim = self.Y.shape - self.Y_metadata = Y_metadata or {} + if Y_metadata is None: + Y_metadata = {} + else: + self.Y_metadata = Y_metadata assert isinstance(kernel, kern.Kern) #assert self.input_dim == kernel.input_dim From 7ac0689156df48a798126bb6be67049a649a4a67 Mon Sep 17 00:00:00 2001 From: Ricardo Date: Mon, 17 Mar 2014 10:27:10 +0000 Subject: [PATCH 06/13] Changes in kernel parameters definition --- GPy/util/multioutput.py | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/GPy/util/multioutput.py b/GPy/util/multioutput.py index 79022a5f..d9e8b704 100644 --- a/GPy/util/multioutput.py +++ b/GPy/util/multioutput.py @@ -35,7 +35,8 @@ def build_likelihood(Y_list,noise_index,likelihoods_list=None): likelihoods_list = [GPy.likelihoods.Gaussian(name="Gaussian_noise_%s" %j) for y,j in zip(Y_list,range(Ny))] else: assert len(likelihoods_list) == Ny - likelihood = GPy.likelihoods.mixed_noise.MixedNoise(likelihoods_list=likelihoods_list, noise_index=noise_index) + #likelihood = GPy.likelihoods.mixed_noise.MixedNoise(likelihoods_list=likelihoods_list, noise_index=noise_index) + likelihood = GPy.likelihoods.mixed_noise.MixedNoise(likelihoods_list=likelihoods_list) return likelihood @@ -43,7 +44,7 @@ def ICM(input_dim, num_outputs, kernel, W_rank=1,W=None,kappa=None,name='X'): """ Builds a kernel for an Intrinsic Coregionalization Model - :input_dim: Input dimensionality + :input_dim: Input dimensionality (does not include dimension of indices) :num_outputs: Number of outputs :param kernel: kernel that will be multiplied by the coregionalize kernel (matrix B). :type kernel: a GPy kernel @@ -54,7 +55,8 @@ def ICM(input_dim, num_outputs, kernel, W_rank=1,W=None,kappa=None,name='X'): kernel.input_dim = input_dim warnings.warn("kernel's input dimension overwritten to fit input_dim parameter.") - K = kernel.prod(GPy.kern.Coregionalize([input_dim], num_outputs,W_rank,W,kappa,name='B'),name=name) + K = kernel.prod(GPy.kern.Coregionalize(1, num_outputs, active_dims=[input_dim], rank=W_rank,W=W,kappa=kappa,name='B'),name=name) + #K = kernel * GPy.kern.Coregionalize(1, num_outputs, active_dims=[input_dim], rank=W_rank,W=W,kappa=kappa,name='B') #K = kernel ** GPy.kern.Coregionalize(input_dim, num_outputs,W_rank,W,kappa, name= 'B') K['.*variance'] = 1. K['.*variance'].fix() @@ -65,7 +67,7 @@ def LCM(input_dim, num_outputs, kernels_list, W_rank=1,name='X'): """ Builds a kernel for an Linear Coregionalization Model - :input_dim: Input dimensionality + :input_dim: Input dimensionality (does not include dimension of indices) :num_outputs: Number of outputs :param kernel: kernel that will be multiplied by the coregionalize kernel (matrix B). :type kernel: a GPy kernel From bb7c26b41699b084915dabd4114a4a48b7a79aa3 Mon Sep 17 00:00:00 2001 From: Ricardo Date: Mon, 17 Mar 2014 10:27:53 +0000 Subject: [PATCH 07/13] Y_metadata definition changed --- GPy/plotting/matplot_dep/models_plots.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/GPy/plotting/matplot_dep/models_plots.py b/GPy/plotting/matplot_dep/models_plots.py index 7507c376..ae79569b 100644 --- a/GPy/plotting/matplot_dep/models_plots.py +++ b/GPy/plotting/matplot_dep/models_plots.py @@ -6,6 +6,7 @@ import numpy as np import Tango from base_plots import gpplot, x_frame1D, x_frame2D from ...util.misc import param_to_array +from ...models.gp_coregionalized_regression import GPCoregionalizedRegression def plot_fit(model, plot_limits=None, which_data_rows='all', @@ -85,8 +86,9 @@ def plot_fit(model, plot_limits=None, which_data_rows='all', lower = m - 2*np.sqrt(v) upper = m + 2*np.sqrt(v) else: - m, v = model.predict(Xgrid, full_cov=False, Y_metadata=Y_metadata) - lower, upper = model.predict_quantiles(Xgrid, Y_metadata=Y_metadata) + meta = {'output_index': Xgrid[:,-1:].astype(np.int)} if isinstance(model,GPCoregionalizedRegression) else None + m, v = model.predict(Xgrid, full_cov=False, Y_metadata=meta) + lower, upper = model.predict_quantiles(Xgrid, Y_metadata=meta) for d in which_data_ycols: From d4476b76c54969be105be43d207fe621e26e0283 Mon Sep 17 00:00:00 2001 From: Ricardo Date: Mon, 17 Mar 2014 10:28:24 +0000 Subject: [PATCH 08/13] Changes in likelihood definition --- GPy/models/gp_coregionalized_regression.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/GPy/models/gp_coregionalized_regression.py b/GPy/models/gp_coregionalized_regression.py index 313e09d4..6d478fd9 100644 --- a/GPy/models/gp_coregionalized_regression.py +++ b/GPy/models/gp_coregionalized_regression.py @@ -31,7 +31,7 @@ class GPCoregionalizedRegression(GP): def __init__(self, X_list, Y_list, kernel=None, likelihoods_list=None, name='GPCR',W_rank=1,kernel_name='X'): #Input and Output - X,Y,self.noise_index = util.multioutput.build_XY(X_list,Y_list) + X,Y,self.output_index = util.multioutput.build_XY(X_list,Y_list) Ny = len(Y_list) #Kernel @@ -39,6 +39,6 @@ class GPCoregionalizedRegression(GP): kernel = util.multioutput.ICM(input_dim=X.shape[1]-1, num_outputs=Ny, kernel=GPy.kern.rbf(X.shape[1]-1), W_rank=1,name=kernel_name) #Likelihood - likelihood = util.multioutput.build_likelihood(Y_list,self.noise_index,likelihoods_list) + likelihood = util.multioutput.build_likelihood(Y_list,self.output_index,likelihoods_list) - super(GPCoregionalizedRegression, self).__init__(X,Y,kernel,likelihood, noise_index=self.noise_index) + super(GPCoregionalizedRegression, self).__init__(X,Y,kernel,likelihood, Y_metadata={'output_index':self.output_index}) From f5a5d4b25ea432d9c1347fb083f705b19c9285ee Mon Sep 17 00:00:00 2001 From: Ricardo Date: Mon, 17 Mar 2014 10:29:04 +0000 Subject: [PATCH 09/13] Re-definition of the week --- GPy/likelihoods/mixed_noise.py | 38 ++++++++++++++++++++++++++-------- 1 file changed, 29 insertions(+), 9 deletions(-) diff --git a/GPy/likelihoods/mixed_noise.py b/GPy/likelihoods/mixed_noise.py index 946cbaf6..bfcb5916 100644 --- a/GPy/likelihoods/mixed_noise.py +++ b/GPy/likelihoods/mixed_noise.py @@ -23,22 +23,22 @@ class MixedNoise(Likelihood): def exact_inference_gradients(self, dL_dKdiag, Y_metadata): assert all([isinstance(l, Gaussian) for l in self.likelihoods_list]) - ind = Y_metadata['output_index'] + ind = Y_metadata['output_index'].flatten() return np.array([dL_dKdiag[ind==i].sum() for i in range(len(self.likelihoods_list))]) def predictive_values(self, mu, var, full_cov=False, Y_metadata=None): if all([isinstance(l, Gaussian) for l in self.likelihoods_list]): - ind = Y_metadata['output_index'] + ind = Y_metadata['output_index'].flatten() _variance = np.array([self.likelihoods_list[j].variance for j in ind ]) if full_cov: var += np.eye(var.shape[0])*_variance - d = 2*np.sqrt(np.diag(var)) - low, up = mu - d, mu + d + #d = 2*np.sqrt(np.diag(var)) + #low, up = mu - d, mu + d else: var += _variance - d = 2*np.sqrt(var) - low, up = mu - d, mu + d - return mu, var, low, up + #d = 2*np.sqrt(var) + #low, up = mu - d, mu + d + return mu, var#, low, up else: raise NotImplementedError @@ -52,8 +52,28 @@ class MixedNoise(Likelihood): def covariance_matrix(self, Y, Y_metadata): assert all([isinstance(l, Gaussian) for l in self.likelihoods_list]) + ind = Y_metadata['output_index'].flatten() variance = np.zeros(Y.shape[0]) - for lik, ind in itertools.izip(self.likelihoods_list, self.likelihoods_indices): - variance[ind] = lik.variance + for lik, j in zip(self.likelihoods_list, range(len(self.likelihoods_list))): + variance[ind==j] = lik.variance return np.diag(variance) + + def samples(self, gp, Y_metadata): + """ + Returns a set of samples of observations based on a given value of the latent variable. + + :param gp: latent variable + """ + N1, N2 = gp.shape + Ysim = np.zeros((N1,N2)) + ind = Y_metadata['output_index'].flatten() + for j in np.unique(ind): + flt = ind==j + gp_filtered = gp[flt,:] + n1 = gp_filtered.shape[0] + lik = self.likelihoods_list[j] + _ysim = np.array([np.random.normal(lik.gp_link.transf(gpj), scale=np.sqrt(lik.variance), size=1) for gpj in gp_filtered.flatten()]) + Ysim[flt,:] = _ysim.reshape(n1,N2) + return Ysim + From 8d98652e8ba401cf20ce8dd9e159830a43c45125 Mon Sep 17 00:00:00 2001 From: Ricardo Date: Mon, 17 Mar 2014 10:29:25 +0000 Subject: [PATCH 10/13] Re-definition of the week --- GPy/likelihoods/likelihood.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/GPy/likelihoods/likelihood.py b/GPy/likelihoods/likelihood.py index 3eafedb1..f8fc17d2 100644 --- a/GPy/likelihoods/likelihood.py +++ b/GPy/likelihoods/likelihood.py @@ -387,13 +387,14 @@ class Likelihood(Parameterized): return pred_mean, pred_var - def predictive_quantiles(self, mu, var, quantiles, Y_metadata): + def predictive_quantiles(self, mu, var, quantiles, Y_metadata=None): #compute the quantiles by sampling!!! N_samp = 1000 s = np.random.randn(mu.shape[0], N_samp)*np.sqrt(var) + mu - ss_f = s.flatten() - ss_y = self.samples(ss_f) - ss_y = ss_y.reshape(mu.shape[0], N_samp) + #ss_f = s.flatten() + #ss_y = self.samples(ss_f, Y_metadata) + ss_y = self.samples(s, Y_metadata) + #ss_y = ss_y.reshape(mu.shape[0], N_samp) return [np.percentile(ss_y ,q, axis=1)[:,None] for q in quantiles] From bb6f9378815e18496122bdcfcd400e87f2d89dac Mon Sep 17 00:00:00 2001 From: Ricardo Date: Mon, 17 Mar 2014 10:30:18 +0000 Subject: [PATCH 11/13] Changes for compatiblity with changes in likelihood --- GPy/likelihoods/gaussian.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/GPy/likelihoods/gaussian.py b/GPy/likelihoods/gaussian.py index aaa356b6..101aac4b 100644 --- a/GPy/likelihoods/gaussian.py +++ b/GPy/likelihoods/gaussian.py @@ -56,7 +56,7 @@ class Gaussian(Likelihood): def update_gradients(self, grad): self.variance.gradient = grad - def exact_inference_gradients(self, dL_dKdiag): + def exact_inference_gradients(self, dL_dKdiag,Y_metadata=None): return dL_dKdiag.sum() def _preprocess_values(self, Y): @@ -295,7 +295,7 @@ class Gaussian(Likelihood): """ return self.variance - def samples(self, gp): + def samples(self, gp, Y_metadata=None): """ Returns a set of samples of observations based on a given value of the latent variable. @@ -303,6 +303,8 @@ class Gaussian(Likelihood): """ orig_shape = gp.shape gp = gp.flatten() + #orig_shape = gp.shape + gp = gp.flatten() Ysim = np.array([np.random.normal(self.gp_link.transf(gpj), scale=np.sqrt(self.variance), size=1) for gpj in gp]) return Ysim.reshape(orig_shape) From f2d5ee42eb05d5aaf517656b46f9dd6a59dd378a Mon Sep 17 00:00:00 2001 From: Zhenwen Dai Date: Mon, 17 Mar 2014 16:09:25 +0000 Subject: [PATCH 12/13] fix the bug regarding to the change of the name dL_dthetaL --- GPy/models/bayesian_gplvm.py | 2 +- GPy/models/ss_gplvm.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/GPy/models/bayesian_gplvm.py b/GPy/models/bayesian_gplvm.py index 3617e260..e0818a2f 100644 --- a/GPy/models/bayesian_gplvm.py +++ b/GPy/models/bayesian_gplvm.py @@ -66,7 +66,7 @@ class BayesianGPLVM(SparseGP): super(BayesianGPLVM, self).parameters_changed() self._log_marginal_likelihood -= self.variational_prior.KL_divergence(self.X) - self.X.mean.gradient, self.X.variance.gradient = self.kern.gradients_qX_expectations(variational_posterior=self.X, Z=self.Z, **self.grad_dict) + self.X.mean.gradient, self.X.variance.gradient = self.kern.gradients_qX_expectations(variational_posterior=self.X, Z=self.Z, dL_dpsi0=self.grad_dict['dL_dpsi0'], dL_dpsi1=self.grad_dict['dL_dpsi1'], dL_dpsi2=self.grad_dict['dL_dpsi2']) # update for the KL divergence self.variational_prior.update_gradients_KL(self.X) diff --git a/GPy/models/ss_gplvm.py b/GPy/models/ss_gplvm.py index 5994814b..1c2ecf4c 100644 --- a/GPy/models/ss_gplvm.py +++ b/GPy/models/ss_gplvm.py @@ -61,7 +61,7 @@ class SSGPLVM(SparseGP): super(SSGPLVM, self).parameters_changed() self._log_marginal_likelihood -= self.variational_prior.KL_divergence(self.X) - self.X.mean.gradient, self.X.variance.gradient, self.X.binary_prob.gradient = self.kern.gradients_qX_expectations(variational_posterior=self.X, Z=self.Z, **self.grad_dict) + self.X.mean.gradient, self.X.variance.gradient, self.X.binary_prob.gradient = self.kern.gradients_qX_expectations(variational_posterior=self.X, Z=self.Z, dL_dpsi0=self.grad_dict['dL_dpsi0'], dL_dpsi1=self.grad_dict['dL_dpsi1'], dL_dpsi2=self.grad_dict['dL_dpsi2']) # update for the KL divergence self.variational_prior.update_gradients_KL(self.X) From caf1dc2609259d48af32d46009d09e8906961427 Mon Sep 17 00:00:00 2001 From: James Hensman Date: Mon, 17 Mar 2014 16:38:04 +0000 Subject: [PATCH 13/13] dL_dthetaL in missing data vardtc --- GPy/inference/latent_function_inference/var_dtc.py | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/GPy/inference/latent_function_inference/var_dtc.py b/GPy/inference/latent_function_inference/var_dtc.py index 7a0c14e8..9b6e26c0 100644 --- a/GPy/inference/latent_function_inference/var_dtc.py +++ b/GPy/inference/latent_function_inference/var_dtc.py @@ -65,7 +65,7 @@ class VarDTC(object): _, output_dim = Y.shape #see whether we've got a different noise variance for each datum - beta = 1./np.fmax(likelihood.variance, 1e-6) + beta = 1./np.fmax(likelihood.gaussian_variance(Y_metadata), 1e-6) # VVT_factor is a matrix such that tdot(VVT_factor) = VVT...this is for efficiency! #self.YYTfactor = self.get_YYTfactor(Y) #VVT_factor = self.get_VVTfactor(self.YYTfactor, beta) @@ -221,7 +221,7 @@ class VarDTCMissingData(object): psi2_all = None Ys, traces = self._Y(Y) - beta_all = 1./np.fmax(likelihood.variance, 1e-6) + beta_all = 1./np.fmax(likelihood.gaussian_variance(Y_metadata), 1e-6) het_noise = beta_all.size != 1 import itertools @@ -328,18 +328,20 @@ class VarDTCMissingData(object): diag.add(Bi, 1) woodbury_inv_all[:, :, ind] = backsub_both_sides(Lm, Bi)[:,:,None] + dL_dthetaL = likelihood.exact_inference_gradients(dL_dR) + # gradients: if uncertain_inputs: grad_dict = {'dL_dKmm': dL_dKmm, 'dL_dpsi0':dL_dpsi0_all, 'dL_dpsi1':dL_dpsi1_all, 'dL_dpsi2':dL_dpsi2_all, - 'dL_dR':dL_dR} + 'dL_dthetaL':dL_dthetaL} else: grad_dict = {'dL_dKmm': dL_dKmm, 'dL_dKdiag':dL_dpsi0_all, 'dL_dKnm':dL_dpsi1_all, - 'dL_dR':dL_dR} + 'dL_dthetaL':dL_dthetaL} #get sufficient things for posterior prediction #TODO: do we really want to do this in the loop?