Merge branch 'master' of github.com:SheffieldML/GPy

GPy/kern/periodic_Matern32.py

@@ -1,6 +1,11 @@
+# Copyright (c) 2012, GPy authors (see AUTHORS.txt).
+# Licensed under the BSD 3-clause license (see LICENSE.txt)
+
+
 from kernpart import kernpart
 import numpy as np
 from GPy.util.linalg import mdot, pdinv
+from GPy.util.decorators import silence_errors
 
 class periodic_Matern32(kernpart):
     """
@@ -39,12 +44,16 @@ class periodic_Matern32(kernpart):
         def f(x):
             return alpha*np.cos(omega*x+phase)
         return f
+
+    @silence_errors
     def _cos_factorization(self,alpha,omega,phase):
         r1 = np.sum(alpha*np.cos(phase),axis=1)[:,None]
         r2 = np.sum(alpha*np.sin(phase),axis=1)[:,None]
         r =  np.sqrt(r1**2 + r2**2)
         psi = np.where(r1 != 0, (np.arctan(r2/r1) + (r1<0.)*np.pi),np.arcsin(r2))
         return r,omega[:,0:1], psi
+
+    @silence_errors
     def _int_computation(self,r1,omega1,phi1,r2,omega2,phi2):
         Gint1 = 1./(omega1+omega2.T)*( np.sin((omega1+omega2.T)*self.upper+phi1+phi2.T) - np.sin((omega1+omega2.T)*self.lower+phi1+phi2.T)) + 1./(omega1-omega2.T)*( np.sin((omega1-omega2.T)*self.upper+phi1-phi2.T) - np.sin((omega1-omega2.T)*self.lower+phi1-phi2.T) )
         Gint2 = 1./(omega1+omega2.T)*( np.sin((omega1+omega2.T)*self.upper+phi1+phi2.T) - np.sin((omega1+omega2.T)*self.lower+phi1+phi2.T)) +  np.cos(phi1-phi2.T)*(self.upper-self.lower)
@@ -55,6 +64,7 @@ class periodic_Matern32(kernpart):
     def _get_params(self):
         """return the value of the parameters."""
         return np.hstack((self.variance,self.lengthscale,self.period))
+    
     def _set_params(self,x):
         """set the value of the parameters."""
         assert x.size==3
@@ -101,6 +111,7 @@ class periodic_Matern32(kernpart):
         FX  = self._cos(self.basis_alpha[None,:],self.basis_omega[None,:],self.basis_phi[None,:])(X)
         np.add(target,np.diag(mdot(FX,self.Gi,FX.T)),target)
 
+    @silence_errors
     def dK_dtheta(self,dL_dK,X,X2,target):
         """derivative of the covariance matrix with respect to the parameters (shape is NxMxNparam)"""
         if X2 is None: X2 = X
@@ -172,6 +183,7 @@ class periodic_Matern32(kernpart):
         #np.add(target[:,:,2],dK_dper, target[:,:,2])
         target[2] += np.sum(dK_dper*dL_dK)
 
+    @silence_errors
     def dKdiag_dtheta(self,dL_dKdiag,X,target):
         """derivative of the diagonal covariance matrix with respect to the parameters"""
         FX  = self._cos(self.basis_alpha[None,:],self.basis_omega[None,:],self.basis_phi[None,:])(X)

GPy/kern/periodic_Matern52.py

@@ -1,6 +1,11 @@
+# Copyright (c) 2012, GPy authors (see AUTHORS.txt).
+# Licensed under the BSD 3-clause license (see LICENSE.txt)
+
+
 from kernpart import kernpart
 import numpy as np
 from GPy.util.linalg import mdot, pdinv
+from GPy.util.decorators import silence_errors
 
 class periodic_Matern52(kernpart):
     """
@@ -40,6 +45,7 @@ class periodic_Matern52(kernpart):
             return alpha*np.cos(omega*x+phase)
         return f
 
+    @silence_errors
     def _cos_factorization(self,alpha,omega,phase):
         r1 = np.sum(alpha*np.cos(phase),axis=1)[:,None]
         r2 = np.sum(alpha*np.sin(phase),axis=1)[:,None]
@@ -57,6 +63,7 @@ class periodic_Matern52(kernpart):
     def _get_params(self):
         """return the value of the parameters."""
         return np.hstack((self.variance,self.lengthscale,self.period))
+    
     def _set_params(self,x):
         """set the value of the parameters."""
         assert x.size==3
@@ -105,6 +112,7 @@ class periodic_Matern52(kernpart):
         FX  = self._cos(self.basis_alpha[None,:],self.basis_omega[None,:],self.basis_phi[None,:])(X)
         np.add(target,np.diag(mdot(FX,self.Gi,FX.T)),target)
 
+    @silence_errors
     def dK_dtheta(self,dL_dK,X,X2,target):
         """derivative of the covariance matrix with respect to the parameters (shape is NxMxNparam)"""
         if X2 is None: X2 = X
@@ -184,6 +192,7 @@ class periodic_Matern52(kernpart):
         #np.add(target[:,:,2],dK_dper, target[:,:,2])
         target[2] += np.sum(dK_dper*dL_dK)
 
+    @silence_errors
     def dKdiag_dtheta(self,dL_dKdiag,X,target):
         """derivative of the diagonal of the covariance matrix with respect to the parameters"""
         FX  = self._cos(self.basis_alpha[None,:],self.basis_omega[None,:],self.basis_phi[None,:])(X)

GPy/kern/periodic_exponential.py

@@ -1,6 +1,11 @@
+# Copyright (c) 2012, GPy authors (see AUTHORS.txt).
+# Licensed under the BSD 3-clause license (see LICENSE.txt)
+
+
 from kernpart import kernpart
 import numpy as np
 from GPy.util.linalg import mdot, pdinv
+from GPy.util.decorators import silence_errors
 
 class periodic_exponential(kernpart):
     """
@@ -40,6 +45,7 @@ class periodic_exponential(kernpart):
             return alpha*np.cos(omega*x+phase)
         return f
 
+    @silence_errors
     def _cos_factorization(self,alpha,omega,phase):
         r1 = np.sum(alpha*np.cos(phase),axis=1)[:,None]
         r2 = np.sum(alpha*np.sin(phase),axis=1)[:,None]
@@ -57,6 +63,7 @@ class periodic_exponential(kernpart):
     def _get_params(self):
         """return the value of the parameters."""
         return np.hstack((self.variance,self.lengthscale,self.period))
+    
     def _set_params(self,x):
         """set the value of the parameters."""
         assert x.size==3
@@ -101,6 +108,7 @@ class periodic_exponential(kernpart):
         FX  = self._cos(self.basis_alpha[None,:],self.basis_omega[None,:],self.basis_phi[None,:])(X)
         np.add(target,np.diag(mdot(FX,self.Gi,FX.T)),target)
 
+    @silence_errors
     def dK_dtheta(self,dL_dK,X,X2,target):
         """derivative of the covariance matrix with respect to the parameters (shape is NxMxNparam)"""
         if X2 is None: X2 = X
@@ -166,6 +174,7 @@ class periodic_exponential(kernpart):
         target[1] += np.sum(dK_dlen*dL_dK)
         target[2] += np.sum(dK_dper*dL_dK)
 
+    @silence_errors
     def dKdiag_dtheta(self,dL_dKdiag,X,target):
         """derivative of the diagonal of the covariance matrix with respect to the parameters"""
         FX  = self._cos(self.basis_alpha[None,:],self.basis_omega[None,:],self.basis_phi[None,:])(X)
@@ -225,4 +234,3 @@ class periodic_exponential(kernpart):
         target[0] += np.sum(np.diag(dK_dvar)*dL_dKdiag)
         target[1] += np.sum(np.diag(dK_dlen)*dL_dKdiag)
         target[2] += np.sum(np.diag(dK_dper)*dL_dKdiag)
-

GPy/models/sparse_GP.py

@@ -153,14 +153,26 @@ class sparse_GP(GP):
             #self.partial_for_likelihood += -np.diag(np.dot((self.C - 0.5 * mdot(self.C,self.psi2_beta_scaled,self.C) ) , self.psi1VVpsi1 ))*self.likelihood.precision #dD
         else:
             #likelihood is not heterscedatic
-            beta = self.likelihood.precision
+            self.partial_for_likelihood =   - 0.5 * self.N*self.D*self.likelihood.precision + 0.5 * np.sum(np.square(self.likelihood.Y))*self.likelihood.precision**2
-            dbeta =   0.5 * self.N*self.D/beta - 0.5 * np.sum(np.square(self.likelihood.Y))
+            self.partial_for_likelihood += 0.5 * self.D * (self.psi0.sum()*self.likelihood.precision**2 - np.trace(self.A)*self.likelihood.precision*sf2)
-            dbeta += - 0.5 * self.D * (self.psi0.sum() - np.trace(self.A)/beta*sf2)
+            self.partial_for_likelihood += 0.5 * self.D * trace_dot(self.Bi,self.A)*self.likelihood.precision
-            dbeta += - 0.5 * self.D * trace_dot(self.Bi,self.A)/beta
+            self.partial_for_likelihood += self.likelihood.precision*(0.5*trace_dot(self.psi2_beta_scaled,self.E*sf2) - np.trace(self.Cpsi1VVpsi1))
-            dbeta += np.trace(self.Cpsi1VVpsi1)/beta - 0.5 * trace_dot(np.dot(self.C,self.psi2_beta_scaled) , self.Cpsi1VVpsi1 )/beta
-            self.partial_for_likelihood = -dbeta*self.likelihood.precision**2
 
 
+
+    def log_likelihood(self):
+        """ Compute the (lower bound on the) log marginal likelihood """
+        sf2 = self.scale_factor**2
+        if self.likelihood.is_heteroscedastic:
+            A = -0.5*self.N*self.D*np.log(2.*np.pi) +0.5*np.sum(np.log(self.likelihood.precision)) -0.5*np.sum(self.V*self.likelihood.Y)
+            B = -0.5*self.D*(np.sum(self.likelihood.precision.flatten()*self.psi0) - np.trace(self.A)*sf2)
+        else:
+            A = -0.5*self.N*self.D*(np.log(2.*np.pi) + np.log(self.likelihood._variance)) -0.5*self.likelihood.precision*self.likelihood.trYYT
+            B = -0.5*self.D*(np.sum(self.likelihood.precision*self.psi0) - np.trace(self.A)*sf2)
+        C = -0.5*self.D * (self.B_logdet + self.M*np.log(sf2))
+        D = 0.5*np.trace(self.Cpsi1VVpsi1)
+        return A+B+C+D
+
     def _set_params(self, p):
         self.Z = p[:self.M*self.Q].reshape(self.M, self.Q)
         self.kern._set_params(p[self.Z.size:self.Z.size+self.kern.Nparam])
@@ -188,18 +200,6 @@ class sparse_GP(GP):
             #self.likelihood.fit_FITC(self.Kmm,self.psi1,self.psi0)
             self._set_params(self._get_params()) # update the GP
 
-    def log_likelihood(self):
-        """ Compute the (lower bound on the) log marginal likelihood """
-        sf2 = self.scale_factor**2
-        if self.likelihood.is_heteroscedastic:
-            A = -0.5*self.N*self.D*np.log(2.*np.pi) +0.5*np.sum(np.log(self.likelihood.precision)) -0.5*np.sum(self.V*self.likelihood.Y)
-            B = -0.5*self.D*(np.sum(self.likelihood.precision.flatten()*self.psi0) - np.trace(self.A)*sf2)
-        else:
-            A = -0.5*self.N*self.D*(np.log(2.*np.pi) - np.log(self.likelihood.precision)) -0.5*self.likelihood.precision*self.likelihood.trYYT
-            B = -0.5*self.D*(np.sum(self.likelihood.precision*self.psi0) - np.trace(self.A)*sf2)
-        C = -0.5*self.D * (self.B_logdet + self.M*np.log(sf2))
-        D = 0.5*np.trace(self.Cpsi1VVpsi1)
-        return A+B+C+D
 
     def _log_likelihood_gradients(self):
         return np.hstack((self.dL_dZ().flatten(), self.dL_dtheta(), self.likelihood._gradients(partial=self.partial_for_likelihood)))