From 82d5479438f300f46783c3b90c1981f4dfac4f7b Mon Sep 17 00:00:00 2001 From: Ricardo Andrade Date: Wed, 23 Jan 2013 11:22:20 +0000 Subject: [PATCH 1/2] Poisson likelihood. --- GPy/inference/likelihoods.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/GPy/inference/likelihoods.py b/GPy/inference/likelihoods.py index 53dd4248..098e525a 100644 --- a/GPy/inference/likelihoods.py +++ b/GPy/inference/likelihoods.py @@ -102,3 +102,5 @@ class probit(likelihood): def _log_likelihood_gradients(): raise NotImplementedError +class poisson(likelihood): + pass From f949e4a2c8a1463b035398d10035bc6c49a3b9b2 Mon Sep 17 00:00:00 2001 From: Ricardo Andrade Date: Wed, 23 Jan 2013 11:25:41 +0000 Subject: [PATCH 2/2] Poisson and Gaussian likelihood. --- GPy/inference/likelihoods.py | 98 +++++++++++++++++++++++++++++++++++- 1 file changed, 97 insertions(+), 1 deletion(-) diff --git a/GPy/inference/likelihoods.py b/GPy/inference/likelihoods.py index 098e525a..0eeeb982 100644 --- a/GPy/inference/likelihoods.py +++ b/GPy/inference/likelihoods.py @@ -103,4 +103,100 @@ class probit(likelihood): raise NotImplementedError class poisson(likelihood): - pass + """ + Poisson likelihood + Y is expected to take values in {0,1,2,...} + ----- + $$ + L(x) = \exp(\lambda) * \lambda**Y_i / Y_i! + $$ + """ + def moments_match(self,i,tau_i,v_i): + """ + Moments match of the marginal approximation in EP algorithm + + :param i: number of observation (int) + :param tau_i: precision of the cavity distribution (float) + :param v_i: mean/variance of the cavity distribution (float) + """ + mu = v_i/tau_i + sigma = np.sqrt(1./tau_i) + def poisson_norm(f): + pdf_norm_f = stats.norm.pdf(f,loc=mu,scale=sigma) + rate = np.exp( (f*self.scale)+self.location) + poisson = stats.poisson.pmf(float(self.Y[i]),rate) + return pdf_norm_f*poisson + + def log_pnm(f): + return -(-.5*(f-mu)**2/sigma**2 - np.exp( (f*self.scale)+self.location) + ( (f*self.scale)+self.location)*self.Y[i]) + + golden_A = -1 if self.Y[i] == 0 else np.array([np.log(self.Y[i]),mu]).min() + golden_B = np.array([np.log(self.Y[i]),mu]).max() + golden_A = (golden_A - self.location)/self.scale + golden_B = (golden_B - self.location)/self.scale + opt = sp.optimize.golden(log_pnm,brack=(golden_A,golden_B)) + width = 3./np.log(max(self.Y[i],2)) + + # Simpson's approximamtion + A = opt - width + B = opt + width + K = 10*int(np.log(max(self.Y[i],150))) + h = (B-A)/K + grid_x = np.hstack([np.linspace(opt-width,opt,K/2+1)[1:-1], np.linspace(opt,opt+width,K/2+1)]) + x = np.hstack([A,B,grid_x[range(1,K,2)],grid_x[range(2,K-1,2)]]) + zeroth = np.hstack([poisson_norm(A),poisson_norm(B),[4*poisson_norm(f) for f in grid_x[range(1,K,2)]],[2*poisson_norm(f) for f in grid_x[range(2,K-1,2)]]]) + first = zeroth*x + second = first*x + Z_hat = sum(zeroth)*h/3 + mu_hat = sum(first)*h/(3*Z_hat) + m2 = sum(second)*h/(3*Z_hat) + sigma2_hat = m2 - mu_hat**2 + 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 + +class gaussian(likelihood): + """ + Gaussian likelihood + Y is expected to take values in (-inf,inf) + """ + def moments_match(self,i,tau_i,v_i): + """ + Moments match of the marginal approximation in EP algorithm + + :param i: number of observation (int) + :param tau_i: precision of the cavity distribution (float) + :param v_i: mean/variance of the cavity distribution (float) + """ + mu = v_i/tau_i + sigma = np.sqrt(1./tau_i) + s = 1. if self.Y[i] == 0 else 1./self.Y[i] + sigma2_hat = 1./(1./sigma**2 + 1./s**2) + mu_hat = sigma2_hat*(mu/sigma**2 + self.Y[i]/s**2) + Z_hat = 1./np.sqrt(2*np.pi) * 1./np.sqrt(sigma**2+s**2) * np.exp(-.5*(mu-self.Y[i])**2/(sigma**2 + s**2)) + return Z_hat, mu_hat, sigma2_hat + + def plot1Db(self,X,X_new,F_new,U=None): + assert X.shape[1] == 1, 'Number of dimensions must be 1' + gpplot(X_new,F_new,np.zeros(X_new.shape[0])) + 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,Sigma): + return mu #return stats.poisson.pmf(h) + + def _log_likelihood_gradients(): + raise NotImplementedError