diff --git a/GPy/examples/coreg_example.py b/GPy/examples/coreg_example.py index 8f4cfcc6..66ba143d 100644 --- a/GPy/examples/coreg_example.py +++ b/GPy/examples/coreg_example.py @@ -23,13 +23,10 @@ K = Bias.prod(Coreg,name='X') #K.coregion.W = 0 #print K.coregion.W - #print Bias.K(_X,_X) #print K.K(X,X) - #pb.matshow(K.K(X,X)) - Mlist = [GPy.kern.Matern32(1,lengthscale=20.,name="Mat")] kern = GPy.util.multioutput.LCM(input_dim=1,num_outputs=2,kernels_list=Mlist,name='H') kern.B.W = 0 @@ -37,16 +34,22 @@ kern.B.kappa = 1. #kern.B.W.fix() #kern.B.kappa.fix() #m = GPy.models.GPCoregionalizedRegression(X_list=[X1,X2], Y_list=[Y1,Y2], kernel=kern) -m = GPy.models.SparseGPCoregionalizedRegression(X_list=[X1], Y_list=[Y1], kernel=kern) + + +Z1 = np.array([1.5,2.5])[:,None] + +m = GPy.models.SparseGPCoregionalizedRegression(X_list=[X1], Y_list=[Y1], Z_list = [Z1], kernel=kern) #m.optimize() m.checkgrad(verbose=1) +""" fig = pb.figure() ax0 = fig.add_subplot(211) ax1 = fig.add_subplot(212) slices = GPy.util.multioutput.get_slices([Y1,Y2]) m.plot(fixed_inputs=[(1,0)],which_data_rows=slices[0],ax=ax0) #m.plot(fixed_inputs=[(1,1)],which_data_rows=slices[1],ax=ax1) +""" diff --git a/GPy/examples/regression.py b/GPy/examples/regression.py index 190af93b..2a4b91b3 100644 --- a/GPy/examples/regression.py +++ b/GPy/examples/regression.py @@ -25,80 +25,51 @@ def olympic_marathon_men(optimize=True, plot=True): return m -def coregionalization_toy2(optimize=True, plot=True): +def coregionalization_toy(optimize=True, plot=True): """ A simple demonstration of coregionalization on two sinusoidal functions. """ #build a design matrix with a column of integers indicating the output X1 = np.random.rand(50, 1) * 8 X2 = np.random.rand(30, 1) * 5 - index = np.vstack((np.zeros_like(X1), np.ones_like(X2))) - X = np.hstack((np.vstack((X1, X2)), index)) #build a suitable set of observed variables Y1 = np.sin(X1) + np.random.randn(*X1.shape) * 0.05 Y2 = np.sin(X2) + np.random.randn(*X2.shape) * 0.05 + 2. - Y = np.vstack((Y1, Y2)) - #build the kernel - k1 = GPy.kern.RBF(1) + GPy.kern.Bias(1) - k2 = GPy.kern.Coregionalize(2,1) - k = k1**k2 - m = GPy.models.GPRegression(X, Y, kernel=k) + m = GPy.models.GPCoregionalizedRegression(X_list=[X1,X2], Y_list=[Y1,Y2]) if optimize: m.optimize('bfgs', max_iters=100) if plot: - m.plot(fixed_inputs=[(1,0)]) - m.plot(fixed_inputs=[(1,1)], ax=pb.gca()) - + slices = GPy.util.multioutput.get_slices([X1,X2]) + m.plot(fixed_inputs=[(1,0)],which_data_rows=slices[0],Y_metadata={'output_index':0}) + m.plot(fixed_inputs=[(1,1)],which_data_rows=slices[1],Y_metadata={'output_index':1},ax=pb.gca()) return m -#FIXME: Needs recovering once likelihoods are consolidated -#def coregionalization_toy(optimize=True, plot=True): -# """ -# A simple demonstration of coregionalization on two sinusoidal functions. -# """ -# X1 = np.random.rand(50, 1) * 8 -# X2 = np.random.rand(30, 1) * 5 -# X = np.vstack((X1, X2)) -# Y1 = np.sin(X1) + np.random.randn(*X1.shape) * 0.05 -# Y2 = -np.sin(X2) + np.random.randn(*X2.shape) * 0.05 -# Y = np.vstack((Y1, Y2)) -# -# k1 = GPy.kern.RBF(1) -# m = GPy.models.GPMultioutputRegression(X_list=[X1,X2],Y_list=[Y1,Y2],kernel_list=[k1]) -# m.constrain_fixed('.*rbf_var', 1.) -# m.optimize(max_iters=100) -# -# fig, axes = pb.subplots(2,1) -# m.plot(fixed_inputs=[(1,0)],ax=axes[0]) -# m.plot(fixed_inputs=[(1,1)],ax=axes[1]) -# axes[0].set_title('Output 0') -# axes[1].set_title('Output 1') -# return m - def coregionalization_sparse(optimize=True, plot=True): """ A simple demonstration of coregionalization on two sinusoidal functions using sparse approximations. """ - #fetch the data from the non sparse examples - m = coregionalization_toy2(optimize=False, plot=False) - X, Y = m.X, m.Y + #build a design matrix with a column of integers indicating the output + X1 = np.random.rand(50, 1) * 8 + X2 = np.random.rand(30, 1) * 5 - k = GPy.kern.RBF(1)**GPy.kern.Coregionalize(2) + #build a suitable set of observed variables + Y1 = np.sin(X1) + np.random.randn(*X1.shape) * 0.05 + Y2 = np.sin(X2) + np.random.randn(*X2.shape) * 0.05 + 2. - #construct a model - m = GPy.models.SparseGPRegression(X,Y, num_inducing=25, kernel=k) - m.Z[:,1].fix() # don't optimize the inducing input indexes + m = GPy.models.SparseGPCoregionalizedRegression(X_list=[X1,X2], Y_list=[Y1,Y2]) if optimize: - m.optimize('bfgs', max_iters=100, messages=1) + m.optimize('bfgs', max_iters=100) if plot: - m.plot(fixed_inputs=[(1,0)]) - m.plot(fixed_inputs=[(1,1)], ax=pb.gca()) + slices = GPy.util.multioutput.get_slices([X1,X2]) + m.plot(fixed_inputs=[(1,0)],which_data_rows=slices[0],Y_metadata={'output_index':0}) + m.plot(fixed_inputs=[(1,1)],which_data_rows=slices[1],Y_metadata={'output_index':1},ax=pb.gca()) + pb.ylim(-3,) return m diff --git a/GPy/inference/latent_function_inference/expectation_propagation.py b/GPy/inference/latent_function_inference/expectation_propagation.py index 514a6dc7..ff60d2e3 100644 --- a/GPy/inference/latent_function_inference/expectation_propagation.py +++ b/GPy/inference/latent_function_inference/expectation_propagation.py @@ -11,9 +11,9 @@ class EP(object): :param epsilon: Convergence criterion, maximum squared difference allowed between mean updates to stop iterations (float) :type epsilon: float - :param eta: Power EP thing TODO: Ricardo: what, exactly? + :param eta: parameter for fractional EP updates. :type eta: float64 - :param delta: Power EP thing TODO: Ricardo: what, exactly? + :param delta: damping EP updates factor. :type delta: float64 """ self.epsilon, self.eta, self.delta = epsilon, eta, delta diff --git a/GPy/inference/latent_function_inference/var_dtc.py b/GPy/inference/latent_function_inference/var_dtc.py index e2aa95f5..669ccd08 100644 --- a/GPy/inference/latent_function_inference/var_dtc.py +++ b/GPy/inference/latent_function_inference/var_dtc.py @@ -176,7 +176,6 @@ class VarDTC(object): #construct a posterior object post = Posterior(woodbury_inv=woodbury_inv, woodbury_vector=woodbury_vector, K=Kmm, mean=None, cov=None, K_chol=Lm) - return post, log_marginal, grad_dict class VarDTCMissingData(object): @@ -365,7 +364,7 @@ class VarDTCMissingData(object): return post, log_marginal, grad_dict def _compute_dL_dpsi(num_inducing, num_data, output_dim, beta, Lm, VVT_factor, Cpsi1Vf, DBi_plus_BiPBi, psi1, het_noise, uncertain_inputs): - dL_dpsi0 = -0.5 * output_dim * (beta * np.ones([num_data, 1])).flatten() + dL_dpsi0 = -0.5 * output_dim * (beta[:,None] * np.ones([num_data, 1])).flatten() dL_dpsi1 = np.dot(VVT_factor, Cpsi1Vf.T) dL_dpsi2_beta = 0.5 * backsub_both_sides(Lm, output_dim * np.eye(num_inducing) - DBi_plus_BiPBi) if het_noise: diff --git a/GPy/likelihoods/gaussian.py b/GPy/likelihoods/gaussian.py index c7001278..6a5030bc 100644 --- a/GPy/likelihoods/gaussian.py +++ b/GPy/likelihoods/gaussian.py @@ -92,7 +92,7 @@ class Gaussian(Likelihood): def predictive_variance(self, mu, sigma, predictive_mean=None): return self.variance + sigma**2 - def predictive_quantiles(self, mu, var, quantiles, Y_metadata): + def predictive_quantiles(self, mu, var, quantiles, Y_metadata=None): return [stats.norm.ppf(q/100.)*np.sqrt(var) + mu for q in quantiles] def pdf_link(self, link_f, y, Y_metadata=None): diff --git a/GPy/likelihoods/likelihood.py b/GPy/likelihoods/likelihood.py index 6df573df..aabe93ef 100644 --- a/GPy/likelihoods/likelihood.py +++ b/GPy/likelihoods/likelihood.py @@ -411,10 +411,7 @@ class Likelihood(Parameterized): #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, 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] diff --git a/GPy/likelihoods/mixed_noise.py b/GPy/likelihoods/mixed_noise.py index 5f4d0705..c2435508 100644 --- a/GPy/likelihoods/mixed_noise.py +++ b/GPy/likelihoods/mixed_noise.py @@ -11,7 +11,7 @@ import itertools class MixedNoise(Likelihood): def __init__(self, likelihoods_list, name='mixed_noise'): - + #NOTE at the moment this likelihood only works for using a list of gaussians super(Likelihood, self).__init__(name=name) self.add_parameters(*likelihoods_list) @@ -24,11 +24,11 @@ class MixedNoise(Likelihood): variance = np.zeros(ind.size) for lik, j in zip(self.likelihoods_list, range(len(self.likelihoods_list))): variance[ind==j] = lik.variance - return variance[:,None] + return variance def betaY(self,Y,Y_metadata): #TODO not here. - return Y/self.gaussian_variance(Y_metadata=Y_metadata) + return Y/self.gaussian_variance(Y_metadata=Y_metadata)[:,None] def update_gradients(self, gradients): self.gradient = gradients @@ -39,24 +39,27 @@ class MixedNoise(Likelihood): 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'].flatten() - _variance = np.array([self.likelihoods_list[j].variance for j in ind ]) - if full_cov: - var += np.eye(var.shape[0])*_variance - else: - var += _variance - return mu, var + 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 else: - raise NotImplementedError + var += _variance + return mu, var - def predictive_variance(self, mu, sigma, **other_shit): - if isinstance(noise_index,int): - _variance = self.variance[noise_index] - else: - _variance = np.array([ self.variance[j] for j in noise_index ])[:,None] + def predictive_variance(self, mu, sigma, Y_metadata): + _variance = self.gaussian_variance(Y_metadata) return _variance + sigma**2 + def predictive_quantiles(self, mu, var, quantiles, Y_metadata): + ind = Y_metadata['output_index'].flatten() + outputs = np.unique(ind) + Q = np.zeros( (mu.size,len(quantiles)) ) + for j in outputs: + q = self.likelihoods_list[j].predictive_quantiles(mu[ind==j,:], + var[ind==j,:],quantiles,Y_metadata=None) + Q[ind==j,:] = np.hstack(q) + return [q[:,None] for q in Q.T] def samples(self, gp, Y_metadata): """ @@ -75,4 +78,3 @@ class MixedNoise(Likelihood): _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 - diff --git a/GPy/models/gp_coregionalized_regression.py b/GPy/models/gp_coregionalized_regression.py index 6d478fd9..5854fe63 100644 --- a/GPy/models/gp_coregionalized_regression.py +++ b/GPy/models/gp_coregionalized_regression.py @@ -36,7 +36,7 @@ class GPCoregionalizedRegression(GP): #Kernel if kernel is None: - 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) + kernel = util.multioutput.ICM(input_dim=X.shape[1]-1, num_outputs=Ny, kernel=kern.RBF(X.shape[1]-1), W_rank=1,name=kernel_name) #Likelihood likelihood = util.multioutput.build_likelihood(Y_list,self.output_index,likelihoods_list) diff --git a/GPy/models/sparse_gp_coregionalized_regression.py b/GPy/models/sparse_gp_coregionalized_regression.py index a97696d2..ea6ec9f1 100644 --- a/GPy/models/sparse_gp_coregionalized_regression.py +++ b/GPy/models/sparse_gp_coregionalized_regression.py @@ -43,14 +43,14 @@ class SparseGPCoregionalizedRegression(SparseGP): #Kernel if kernel is None: - 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) + kernel = util.multioutput.ICM(input_dim=X.shape[1]-1, num_outputs=Ny, kernel=kern.RBF(X.shape[1]-1), W_rank=1,name=kernel_name) #Likelihood likelihood = util.multioutput.build_likelihood(Y_list,self.output_index,likelihoods_list) #Inducing inputs list if len(Z_list): - assert len(Z_list) == self.output_dim, 'Number of outputs do not match length of inducing inputs list.' + assert len(Z_list) == Ny, 'Number of outputs do not match length of inducing inputs list.' else: if isinstance(num_inducing,np.int): num_inducing = [num_inducing] * Ny diff --git a/GPy/plotting/matplot_dep/models_plots.py b/GPy/plotting/matplot_dep/models_plots.py index cbb213b1..b626758f 100644 --- a/GPy/plotting/matplot_dep/models_plots.py +++ b/GPy/plotting/matplot_dep/models_plots.py @@ -123,6 +123,8 @@ def plot_fit(model, plot_limits=None, which_data_rows='all', #add inducing inputs (if a sparse model is used) if hasattr(model,"Z"): #Zu = model.Z[:,free_dims] * model._Xscale[:,free_dims] + model._Xoffset[:,free_dims] + if isinstance(model,SparseGPCoregionalizedRegression): + Z = Z[Z[:,-1] == Y_metadata['output_index'],:] Zu = Z[:,free_dims] z_height = ax.get_ylim()[0] plots['inducing_inputs'] = ax.plot(Zu, np.zeros_like(Zu) + z_height, 'r|', mew=1.5, markersize=12) diff --git a/GPy/util/multioutput.py b/GPy/util/multioutput.py index d9e8b704..b769d6a5 100644 --- a/GPy/util/multioutput.py +++ b/GPy/util/multioutput.py @@ -56,8 +56,6 @@ def ICM(input_dim, num_outputs, kernel, W_rank=1,W=None,kappa=None,name='X'): warnings.warn("kernel's input dimension overwritten to fit input_dim parameter.") 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() return K