New rational quadratic kernel

GPy/kern/__init__.py

@@ -2,5 +2,5 @@
 # Licensed under the BSD 3-clause license (see LICENSE.txt)
 
 
-from constructors import rbf, Matern32, Matern52, exponential, linear, white, bias, finite_dimensional, spline, Brownian, rbf_sympy, sympykern, periodic_exponential, periodic_Matern32, periodic_Matern52, prod, prod_orthogonal, symmetric, coregionalise
+from constructors import rbf, Matern32, Matern52, exponential, linear, white, bias, finite_dimensional, spline, Brownian, rbf_sympy, sympykern, periodic_exponential, periodic_Matern32, periodic_Matern52, prod, prod_orthogonal, symmetric, coregionalise, rational_quadratic
 from kern import kern

GPy/kern/constructors.py

@@ -22,6 +22,7 @@ from prod import prod as prodpart
 from prod_orthogonal import prod_orthogonal as prod_orthogonalpart
 from symmetric import symmetric as symmetric_part
 from coregionalise import coregionalise as coregionalise_part
+from rational_quadratic import rational_quadratic as rational_quadraticpart
 #TODO these s=constructors are not as clean as we'd like. Tidy the code up
 #using meta-classes to make the objects construct properly wthout them.
 
@@ -280,3 +281,18 @@ def coregionalise(Nout,R=1, W=None, kappa=None):
     return kern(1,[p])
 
 
+def rational_quadratic(D,variance=1., lengthscale=1., power=1.):
+    """
+     Construct rational quadratic kernel.
+
+    :param D: the number of input dimensions
+    :type D: int (D=1 is the only value currently supported)
+    :param variance: the variance :math:`\sigma^2`
+    :type variance: float
+    :param lengthscale: the lengthscale :math:`\ell`
+    :type lengthscale: float
+    :rtype: kern object
+
+    """
+    part = rational_quadraticpart(D,variance, lengthscale, power)
+    return kern(D, [part])

GPy/kern/rational_quadratic.py

@@ -0,0 +1,79 @@
+# 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
+
+class rational_quadratic(kernpart):
+    """
+    rational quadratic kernel
+
+    .. math::
+
+       k(r) = \sigma^2 \left(1 + \frac{r^2}{2 \ell^2})^{- \alpha} \ \ \ \ \  \\text{ where  } r^2 = (x-y)^2
+
+    :param D: the number of input dimensions
+    :type D: int (D=1 is the only value currently supported)
+    :param variance: the variance :math:`\sigma^2`
+    :type variance: float
+    :param lengthscale: the lengthscale :math:`\ell`
+    :type lengthscale: float
+    :rtype: kernpart object
+
+    """
+    def __init__(self,D,variance=1.,lengthscale=1.,power=1.):
+        assert D == 1, "For this kernel we assume D=1"
+        self.D = D
+        self.Nparam = 3
+        self.name = 'rat_quad'
+        self.variance = variance
+        self.lengthscale = lengthscale
+        self.power = power
+
+    def _get_params(self):
+        return np.hstack((self.variance,self.lengthscale,self.power))
+
+    def _set_params(self,x):
+        self.variance = x[0]
+        self.lengthscale = x[1]
+        self.power = x[2]
+
+    def _get_param_names(self):
+        return ['variance','lengthscale','power']
+
+    def K(self,X,X2,target):
+        if X2 is None: X2 = X
+        dist2 = np.square((X-X2.T)/self.lengthscale)
+        target += self.variance*(1 + dist2/2.)**(-self.power)
+
+    def Kdiag(self,X,target):
+        target += self.variance
+
+    def dK_dtheta(self,dL_dK,X,X2,target):
+        if X2 is None: X2 = X
+        dist2 = np.square((X-X2.T)/self.lengthscale)
+
+        dvar = (1 + dist2/2.)**(-self.power)
+        dl = self.power * self.variance * dist2 * self.lengthscale**(-3) * (1 + dist2/2./self.power)**(-self.power-1)
+        dp = - self.variance * np.log(1 + dist2/2.) * (1 + dist2/2.)**(-self.power)
+
+        target[0] += np.sum(dvar*dL_dK)
+        target[1] += np.sum(dl*dL_dK)
+        target[2] += np.sum(dp*dL_dK)
+
+    def dKdiag_dtheta(self,dL_dKdiag,X,target):
+        target[0] += np.sum(dL_dKdiag)
+        # here self.lengthscale and self.power have no influence on Kdiag so target[1:] are unchanged
+
+    def dK_dX(self,dL_dK,X,X2,target):
+        """derivative of the covariance matrix with respect to X."""
+        if X2 is None: X2 = X
+        dist2 = np.square((X-X2.T)/self.lengthscale)
+
+        dX = -self.variance*self.power * (X-X2.T)/self.lengthscale**2 *  (1 + dist2/2./self.power)**(-self.power-1)
+        target += np.sum(dL_dK*dX)
+
+    def dKdiag_dX(self,dL_dKdiag,X,target):
+        pass
+