diff --git a/GPy/examples/poisson.py b/GPy/examples/poisson.py index 71d80b30..e15f310d 100644 --- a/GPy/examples/poisson.py +++ b/GPy/examples/poisson.py @@ -25,16 +25,14 @@ seed=default_seed """ X = np.arange(0,100,5)[:,None] -F = np.round(np.sin(X/18.) + .1*X) -E = np.random.randint(-3,3,20)[:,None] +F = np.round(np.sin(X/18.) + .1*X) + np.arange(5,25)[:,None] +E = np.random.randint(-5,5,20)[:,None] Y = F + E -pb.plot(X,F,'k-') -pb.plot(X,Y,'ro') pb.figure() -likelihood = GPy.inference.likelihoods.poisson(Y,scale=6.) +likelihood = GPy.inference.likelihoods.poisson(Y,scale=1.) m = GPy.models.GP(X,likelihood=likelihood) -#m = GPy.models.GP(data['X'],Y=likelihood.Y) +#m = GPy.models.GP(X,Y=likelihood.Y) m.constrain_positive('var') m.constrain_positive('len') diff --git a/GPy/inference/EP.py b/GPy/inference/EP.py index 5d571888..5c473a8f 100644 --- a/GPy/inference/EP.py +++ b/GPy/inference/EP.py @@ -48,13 +48,13 @@ class EP: self.tau_tilde = np.zeros(self.N) self.v_tilde = np.zeros(self.N) - def restart_EP(self): - """ - Set the EP approximation to initial state - """ - self.tau_tilde = np.zeros(self.N) - self.v_tilde = np.zeros(self.N) - self.mu = np.zeros(self.N) + def _compute_GP_variables(self): + #Variables to be called from GP + mu_tilde = self.v_tilde/self.tau_tilde #When calling EP, this variable is used instead of Y in the GP model + sigma_sum = 1./self.tau_ + 1./self.tau_tilde + mu_diff_2 = (self.v_/self.tau_ - mu_tilde)**2 + Z_ep = np.sum(np.log(self.Z_hat)) + 0.5*np.sum(np.log(sigma_sum)) + 0.5*np.sum(mu_diff_2/sigma_sum) #Normalization constant + return self.tau_tilde[:,None], mu_tilde[:,None], Z_ep class Full(EP): def fit_EP(self): @@ -122,12 +122,7 @@ class Full(EP): self.np1.append(self.tau_tilde.copy()) self.np2.append(self.v_tilde.copy()) - #Variables to be called from GP - mu_tilde = self.v_tilde/self.tau_tilde #When calling EP, this variable is used instead of Y in the GP model - sigma_sum = 1./self.tau_ + 1./self.tau_tilde - mu_diff_2 = (self.v_/self.tau_ - mu_tilde)**2 - Z_ep = np.sum(np.log(self.Z_hat)) + 0.5*np.sum(np.log(sigma_sum)) + 0.5*np.sum(mu_diff_2/sigma_sum) #Normalization constant - return self.tau_tilde[:,None], mu_tilde[:,None], Z_ep + return self._compute_GP_variables() class DTC(EP): def fit_EP(self): @@ -212,7 +207,8 @@ class DTC(EP): epsilon_np2 = sum((self.v_tilde-self.np2[-1])**2)/self.N self.np1.append(self.tau_tilde.copy()) self.np2.append(self.v_tilde.copy()) - return self.tau_tilde[:,None], self.v_tilde[:,None], self.Z_hat[:,None], self.tau_[:,None], self.v_[:,None] + + return self._compute_GP_variables() class FITC(EP): def fit_EP(self): @@ -313,4 +309,5 @@ class FITC(EP): epsilon_np2 = sum((self.v_tilde-self.np2[-1])**2)/self.N self.np1.append(self.tau_tilde.copy()) self.np2.append(self.v_tilde.copy()) - return self.tau_tilde[:,None], self.v_tilde[:,None], self.Z_hat[:,None], self.tau_[:,None], self.v_[:,None] + + return self._compute_GP_variables() diff --git a/GPy/inference/likelihoods.py b/GPy/inference/likelihoods.py index 864afa57..b170dc3d 100644 --- a/GPy/inference/likelihoods.py +++ b/GPy/inference/likelihoods.py @@ -9,65 +9,18 @@ import pylab as pb from ..util.plot import gpplot class likelihood: - def __init__(self,Y,location=0,scale=1): - """ - Likelihood class for doing Expectation propagation + """ + Likelihood class for doing Expectation propagation - :param Y: observed output (Nx1 numpy.darray) - ..Note:: Y values allowed depend on the likelihood used - """ + :param Y: observed output (Nx1 numpy.darray) + ..Note:: Y values allowed depend on the likelihood used + """ + def __init__(self,Y,location=0,scale=1): self.Y = Y self.N = self.Y.shape[0] self.location = location self.scale = scale - def plot1D(self,X,mean,var,Z=None,mean_Z=None,var_Z=None,samples=0): - """ - Plot the predictive distribution of the GP model for 1-dimensional inputs - - :param X: The points at which to make a prediction - :param Mean: mean values at X - :param Var: variance values at X - :param Z: Set of points to be highlighted in the plot, i.e. inducing points - :param mean_Z: mean values at Z - :param var_Z: variance values at Z - :samples: Number of samples to plot - """ - assert X.shape[1] == 1, 'Number of dimensions must be 1' - gpplot(X,mean,var.flatten()) - if samples: #NOTE why don't we put samples as a parameter of gpplot - s = np.random.multivariate_normal(mean.flatten(),np.diag(var),samples) - pb.plot(X.flatten(),s.T, alpha = 0.4, c='#3465a4', linewidth = 0.8) - #pb.subplot(211) - #self.plot1Da(X,mean,var,Z,mean_Z,var_Z) - - - def plot1Da(self,X,mean,var,Z=None,mean_Z=None,var_Z=None): - """ - Plot the predictive distribution of the GP model for 1-dimensional inputs - - :param X_new: The points at which to make a prediction - :param Mean_new: mean values at X_new - :param Var_new: variance values at X_new - :param X_u: input (inducing) points used to train the model - :param Mean_u: mean values at X_u - :param Var_new: variance values at X_u - """ - assert X.shape[1] == 1, 'Number of dimensions must be 1' - gpplot(X,mean,var.flatten()) - pb.errorbar(Z.flatten(),mean_Z.flatten(),2*np.sqrt(var_Z.flatten()),fmt='r+') - pb.plot(Z,mean_Z,'ro') - - """ - def plot1Db(self,X_obs,X,phi,Z=None): - assert X_obs.shape[1] == 1, 'Number of dimensions must be 1' - gpplot(X,phi,np.zeros(X.shape[0])) - pb.plot(X_obs,(self.Y+1)/2,'kx',mew=1.5) - pb.ylim(-0.2,1.2) - if Z is not None: - pb.plot(Z,Z*0+.5,'r|',mew=1.5,markersize=12) - - """ def plot2D(self,X,X_new,F_new,U=None): """ Predictive distribution of the fitted GP model for 2-dimensional inputs @@ -106,6 +59,10 @@ class probit(likelihood): L(x) = \\Phi (Y_i*f_i) $$ """ + def __init__(self,Y,location=0,scale=1): + assert np.sum(np.abs(Y)-1) == 0, "Output values must be either -1 or 1" + likelihood.__init__(self,Y,location,scale) + def moments_match(self,i,tau_i,v_i): """ Moments match of the marginal approximation in EP algorithm @@ -146,6 +103,10 @@ class poisson(likelihood): L(x) = \exp(\lambda) * \lambda**Y_i / Y_i! $$ """ + def __init__(self,Y,location=0,scale=1): + assert len(Y[Y<0]) == 0, "Output cannot have negative values" + likelihood.__init__(self,Y,location,scale) + def moments_match(self,i,tau_i,v_i): """ Moments match of the marginal approximation in EP algorithm @@ -203,20 +164,19 @@ class poisson(likelihood): sigma2_hat = m2 - mu_hat**2 # Second central moment return float(Z_hat), float(mu_hat), float(sigma2_hat) - def plot1Db(self,X,X_new,F_new,F2_new=None,U=None): - pb.subplot(212) - #gpplot(X_new,F_new,np.sqrt(F2_new)) - pb.plot(X_new,F_new)#,np.sqrt(F2_new)) #FIXME - pb.plot(X,self.Y,'kx',mew=1.5) - if U is not None: - pb.plot(U,np.ones(U.shape[0])*self.Y.min()*.8,'r|',mew=1.5,markersize=12) def predictive_mean(self,mu,variance): return np.exp(mu*self.scale + self.location) - def predictive_variance(self,mu,variance): - return mu + def _log_likelihood_gradients(): raise NotImplementedError + def plot(self,X,phi,X_obs,Z=None): + assert X_obs.shape[1] == 1, 'Number of dimensions must be 1' + gpplot(X,phi,np.zeros(X.shape[0])) + pb.plot(X_obs,self.Y,'kx',mew=1.5) + if Z is not None: + pb.plot(Z,Z*0+pb.ylim()[0],'k|',mew=1.5,markersize=12) + class gaussian(likelihood): """ Gaussian likelihood diff --git a/GPy/models/sparse_GP.py b/GPy/models/sparse_GP.py index 655f6252..f5381eed 100644 --- a/GPy/models/sparse_GP.py +++ b/GPy/models/sparse_GP.py @@ -35,9 +35,13 @@ class sparse_GP(GP): :type beta: float :param normalize_(X|Y) : whether to normalize the data before computing (predictions will be in original scales) :type normalize_(X|Y): bool + :parm likelihood: a GPy likelihood, defaults to gaussian + :param epsilon_ep: convergence criterion for the Expectation Propagation algorithm, defaults to 0.1 + :param powerep: power-EP parameters [$\eta$,$\delta$], defaults to [1.,1.] + :type powerep: list """ - def __init__(self,X,Y=None,kernel=None, X_uncertainty=None, beta=100., Z=None,Zslices=None,M=10,normalize_X=False,normalize_Y=False,likelihood=None,method_ep='DTC',epsilon_ep=1e-3,epsilon_em=.1,power_ep=[1.,1.]): + def __init__(self,X,Y=None,kernel=None,X_uncertainty=None,beta=100.,Z=None,Zslices=None,M=10,normalize_X=False,normalize_Y=False,likelihood=None,method_ep='DTC',epsilon_ep=1e-3,power_ep=[1.,1.]): if Z is None: self.Z = np.random.permutation(X.copy())[:M] @@ -53,7 +57,7 @@ class sparse_GP(GP): self.has_uncertain_inputs=True self.X_uncertainty = X_uncertainty - GP.__init__(self, X=X, Y=Y, kernel=kernel, normalize_X=normalize_X, normalize_Y=normalize_Y,likelihood=likelihood,epsilon_ep=epsilon_ep,epsilon_em=epsilon_em,power_ep=power_ep) + GP.__init__(self, X=X, Y=Y, kernel=kernel, normalize_X=normalize_X, normalize_Y=normalize_Y,likelihood=likelihood,epsilon_ep=epsilon_ep,power_ep=power_ep) self.trYYT = np.sum(np.square(self.Y)) if not self.EP else None @@ -91,6 +95,9 @@ class sparse_GP(GP): else: if self.hetero_noise: print "rick's stuff here" + + + else: self.psi0 = self.kern.Kdiag(self.X,slices=self.Xslices).sum() self.psi1 = self.kern.K(self.Z,self.X)