Merge remote-tracking branch 'rick70x7/master'

GPy/inference/likelihoods.py

@@ -102,3 +102,101 @@ class probit(likelihood):
     def _log_likelihood_gradients():
         raise NotImplementedError
 
+class poisson(likelihood):
+    """
+    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