mirror of
https://github.com/SheffieldML/GPy.git
synced 2026-05-07 02:52:40 +02:00
Merge branch 'devel' into gpy_one_fixes
This commit is contained in:
Max Zwiessele
commit
f4f3d8546e
32 changed files with 36241 additions and 6 deletions
|
|
@ -29,3 +29,11 @@ from .src.splitKern import SplitKern,DEtime
|
|||
from .src.splitKern import DEtime as DiffGenomeKern
|
||||
from .src.spline import Spline
|
||||
from .src.basis_funcs import LogisticBasisFuncKernel, LinearSlopeBasisFuncKernel, BasisFuncKernel, ChangePointBasisFuncKernel, DomainKernel
|
||||
|
||||
from .src.sde_matern import sde_Matern32
|
||||
from .src.sde_matern import sde_Matern52
|
||||
from .src.sde_linear import sde_Linear
|
||||
from .src.sde_standard_periodic import sde_StdPeriodic
|
||||
from .src.sde_static import sde_White, sde_Bias
|
||||
from .src.sde_stationary import sde_RBF,sde_Exponential,sde_RatQuad
|
||||
from .src.sde_brownian import sde_Brownian
|
||||
|
|
|
|||
57
GPy/kern/_src/sde_brownian.py
Normal file
57
GPy/kern/_src/sde_brownian.py
Normal file
|
|
@ -0,0 +1,57 @@
|
|||
# -*- coding: utf-8 -*-
|
||||
"""
|
||||
Classes in this module enhance Brownian motion covariance function with the
|
||||
Stochastic Differential Equation (SDE) functionality.
|
||||
"""
|
||||
|
||||
from .brownian import Brownian
|
||||
|
||||
import numpy as np
|
||||
|
||||
class sde_Brownian(Brownian):
|
||||
"""
|
||||
|
||||
Class provide extra functionality to transfer this covariance function into
|
||||
SDE form.
|
||||
|
||||
Linear kernel:
|
||||
|
||||
.. math::
|
||||
|
||||
k(x,y) = \sigma^2 min(x,y)
|
||||
|
||||
"""
|
||||
|
||||
def sde_update_gradient_full(self, gradients):
|
||||
"""
|
||||
Update gradient in the order in which parameters are represented in the
|
||||
kernel
|
||||
"""
|
||||
|
||||
self.variance.gradient = gradients[0]
|
||||
|
||||
def sde(self):
|
||||
"""
|
||||
Return the state space representation of the covariance.
|
||||
"""
|
||||
|
||||
variance = float(self.variance.values) # this is initial variancve in Bayesian linear regression
|
||||
|
||||
F = np.array( ((0,1.0),(0,0) ))
|
||||
L = np.array( ((1.0,),(0,)) )
|
||||
Qc = np.array( ((variance,),) )
|
||||
H = np.array( ((1.0,0),) )
|
||||
|
||||
Pinf = np.array( ( (0, -0.5*variance ), (-0.5*variance, 0) ) )
|
||||
#P0 = Pinf.copy()
|
||||
P0 = np.zeros((2,2))
|
||||
#Pinf = np.array( ( (t0, 1.0), (1.0, 1.0/t0) ) ) * variance
|
||||
dF = np.zeros((2,2,1))
|
||||
dQc = np.ones( (1,1,1) )
|
||||
|
||||
dPinf = np.zeros((2,2,1))
|
||||
dPinf[:,:,0] = np.array( ( (0, -0.5), (-0.5, 0) ) )
|
||||
#dP0 = dPinf.copy()
|
||||
dP0 = np.zeros((2,2,1))
|
||||
|
||||
return (F, L, Qc, H, Pinf, P0, dF, dQc, dPinf, dP0)
|
||||
64
GPy/kern/_src/sde_linear.py
Normal file
64
GPy/kern/_src/sde_linear.py
Normal file
|
|
@ -0,0 +1,64 @@
|
|||
# -*- coding: utf-8 -*-
|
||||
"""
|
||||
Classes in this module enhance Linear covariance function with the
|
||||
Stochastic Differential Equation (SDE) functionality.
|
||||
"""
|
||||
from .linear import Linear
|
||||
|
||||
import numpy as np
|
||||
|
||||
class sde_Linear(Linear):
|
||||
"""
|
||||
|
||||
Class provide extra functionality to transfer this covariance function into
|
||||
SDE form.
|
||||
|
||||
Linear kernel:
|
||||
|
||||
.. math::
|
||||
|
||||
k(x,y) = \sum_{i=1}^{input dim} \sigma^2_i x_iy_i
|
||||
|
||||
"""
|
||||
def __init__(self, input_dim, X, variances=None, ARD=False, active_dims=None, name='linear'):
|
||||
"""
|
||||
Modify the init method, because one extra parameter is required. X - points
|
||||
on the X axis.
|
||||
"""
|
||||
|
||||
super(sde_Linear, self).__init__(input_dim, variances, ARD, active_dims, name)
|
||||
|
||||
self.t0 = np.min(X)
|
||||
|
||||
|
||||
def sde_update_gradient_full(self, gradients):
|
||||
"""
|
||||
Update gradient in the order in which parameters are represented in the
|
||||
kernel
|
||||
"""
|
||||
|
||||
self.variances.gradient = gradients[0]
|
||||
|
||||
def sde(self):
|
||||
"""
|
||||
Return the state space representation of the covariance.
|
||||
"""
|
||||
|
||||
variance = float(self.variances.values) # this is initial variancve in Bayesian linear regression
|
||||
t0 = float(self.t0)
|
||||
|
||||
F = np.array( ((0,1.0),(0,0) ))
|
||||
L = np.array( ((0,),(1.0,)) )
|
||||
Qc = np.zeros((1,1))
|
||||
H = np.array( ((1.0,0),) )
|
||||
|
||||
Pinf = np.zeros((2,2))
|
||||
P0 = np.array( ( (t0**2, t0), (t0, 1) ) ) * variance
|
||||
dF = np.zeros((2,2,1))
|
||||
dQc = np.zeros( (1,1,1) )
|
||||
|
||||
dPinf = np.zeros((2,2,1))
|
||||
dP0 = np.zeros((2,2,1))
|
||||
dP0[:,:,0] = P0 / variance
|
||||
|
||||
return (F, L, Qc, H, Pinf, P0, dF, dQc, dPinf, dP0)
|
||||
135
GPy/kern/_src/sde_matern.py
Normal file
135
GPy/kern/_src/sde_matern.py
Normal file
|
|
@ -0,0 +1,135 @@
|
|||
# -*- coding: utf-8 -*-
|
||||
"""
|
||||
Classes in this module enhance Matern covariance functions with the
|
||||
Stochastic Differential Equation (SDE) functionality.
|
||||
"""
|
||||
from .stationary import Matern32
|
||||
from .stationary import Matern52
|
||||
import numpy as np
|
||||
|
||||
class sde_Matern32(Matern32):
|
||||
"""
|
||||
|
||||
Class provide extra functionality to transfer this covariance function into
|
||||
SDE forrm.
|
||||
|
||||
Matern 3/2 kernel:
|
||||
|
||||
.. math::
|
||||
|
||||
k(r) = \sigma^2 (1 + \sqrt{3} r) \exp(- \sqrt{3} r) \\ \\ \\ \\ \text{ where } r = \sqrt{\sum_{i=1}^{input dim} \frac{(x_i-y_i)^2}{\ell_i^2} }
|
||||
|
||||
"""
|
||||
def sde_update_gradient_full(self, gradients):
|
||||
"""
|
||||
Update gradient in the order in which parameters are represented in the
|
||||
kernel
|
||||
"""
|
||||
|
||||
self.variance.gradient = gradients[0]
|
||||
self.lengthscale.gradient = gradients[1]
|
||||
|
||||
def sde(self):
|
||||
"""
|
||||
Return the state space representation of the covariance.
|
||||
"""
|
||||
|
||||
variance = float(self.variance.values)
|
||||
lengthscale = float(self.lengthscale.values)
|
||||
|
||||
foo = np.sqrt(3.)/lengthscale
|
||||
F = np.array(((0, 1.0), (-foo**2, -2*foo)))
|
||||
L = np.array(( (0,), (1.0,) ))
|
||||
Qc = np.array(((12.*np.sqrt(3) / lengthscale**3 * variance,),))
|
||||
H = np.array(((1.0, 0),))
|
||||
Pinf = np.array(((variance, 0.0), (0.0, 3.*variance/(lengthscale**2))))
|
||||
P0 = Pinf.copy()
|
||||
|
||||
# Allocate space for the derivatives
|
||||
dF = np.empty([F.shape[0],F.shape[1],2])
|
||||
dQc = np.empty([Qc.shape[0],Qc.shape[1],2])
|
||||
dPinf = np.empty([Pinf.shape[0],Pinf.shape[1],2])
|
||||
# The partial derivatives
|
||||
dFvariance = np.zeros((2,2))
|
||||
dFlengthscale = np.array(((0,0), (6./lengthscale**3,2*np.sqrt(3)/lengthscale**2)))
|
||||
dQcvariance = np.array((12.*np.sqrt(3)/lengthscale**3))
|
||||
dQclengthscale = np.array((-3*12*np.sqrt(3)/lengthscale**4*variance))
|
||||
dPinfvariance = np.array(((1,0),(0,3./lengthscale**2)))
|
||||
dPinflengthscale = np.array(((0,0), (0,-6*variance/lengthscale**3)))
|
||||
# Combine the derivatives
|
||||
dF[:,:,0] = dFvariance
|
||||
dF[:,:,1] = dFlengthscale
|
||||
dQc[:,:,0] = dQcvariance
|
||||
dQc[:,:,1] = dQclengthscale
|
||||
dPinf[:,:,0] = dPinfvariance
|
||||
dPinf[:,:,1] = dPinflengthscale
|
||||
dP0 = dPinf.copy()
|
||||
|
||||
return (F, L, Qc, H, Pinf, P0, dF, dQc, dPinf, dP0)
|
||||
|
||||
class sde_Matern52(Matern52):
|
||||
"""
|
||||
|
||||
Class provide extra functionality to transfer this covariance function into
|
||||
SDE forrm.
|
||||
|
||||
Matern 5/2 kernel:
|
||||
|
||||
.. math::
|
||||
|
||||
k(r) = \sigma^2 (1 + \sqrt{5} r + \frac{5}{3}r^2) \exp(- \sqrt{5} r) \\ \\ \\ \\ \text{ where } r = \sqrt{\sum_{i=1}^{input dim} \frac{(x_i-y_i)^2}{\ell_i^2} }
|
||||
|
||||
"""
|
||||
def sde_update_gradient_full(self, gradients):
|
||||
"""
|
||||
Update gradient in the order in which parameters are represented in the
|
||||
kernel
|
||||
"""
|
||||
|
||||
self.variance.gradient = gradients[0]
|
||||
self.lengthscale.gradient = gradients[1]
|
||||
|
||||
def sde(self):
|
||||
"""
|
||||
Return the state space representation of the covariance.
|
||||
"""
|
||||
|
||||
variance = float(self.variance.values)
|
||||
lengthscale = float(self.lengthscale.values)
|
||||
|
||||
lamda = np.sqrt(5.0)/lengthscale
|
||||
kappa = 5.0/3.0*variance/lengthscale**2
|
||||
|
||||
F = np.array(((0, 1,0), (0, 0, 1), (-lamda**3, -3.0*lamda**2, -3*lamda)))
|
||||
L = np.array(((0,),(0,),(1,)))
|
||||
Qc = np.array((((variance*400.0*np.sqrt(5.0)/3.0/lengthscale**5),),))
|
||||
H = np.array(((1,0,0),))
|
||||
|
||||
Pinf = np.array(((variance,0,-kappa), (0, kappa, 0), (-kappa, 0, 25.0*variance/lengthscale**4)))
|
||||
P0 = Pinf.copy()
|
||||
# Allocate space for the derivatives
|
||||
dF = np.empty((3,3,2))
|
||||
dQc = np.empty((1,1,2))
|
||||
dPinf = np.empty((3,3,2))
|
||||
|
||||
# The partial derivatives
|
||||
dFvariance = np.zeros((3,3))
|
||||
dFlengthscale = np.array(((0,0,0),(0,0,0),(15.0*np.sqrt(5.0)/lengthscale**4,
|
||||
30.0/lengthscale**3, 3*np.sqrt(5.0)/lengthscale**2)))
|
||||
dQcvariance = np.array((((400*np.sqrt(5)/3/lengthscale**5,),)))
|
||||
dQclengthscale = np.array((((-variance*2000*np.sqrt(5)/3/lengthscale**6,),)))
|
||||
|
||||
dPinf_variance = Pinf/variance
|
||||
kappa2 = -2.0*kappa/lengthscale
|
||||
dPinf_lengthscale = np.array(((0,0,-kappa2),(0,kappa2,0),(-kappa2,
|
||||
0,-100*variance/lengthscale**5)))
|
||||
# Combine the derivatives
|
||||
dF[:,:,0] = dFvariance
|
||||
dF[:,:,1] = dFlengthscale
|
||||
dQc[:,:,0] = dQcvariance
|
||||
dQc[:,:,1] = dQclengthscale
|
||||
dPinf[:,:,0] = dPinf_variance
|
||||
dPinf[:,:,1] = dPinf_lengthscale
|
||||
dP0 = dPinf.copy()
|
||||
|
||||
return (F, L, Qc, H, Pinf, P0, dF, dQc, dPinf, dP0)
|
||||
178
GPy/kern/_src/sde_standard_periodic.py
Normal file
178
GPy/kern/_src/sde_standard_periodic.py
Normal file
|
|
@ -0,0 +1,178 @@
|
|||
# -*- coding: utf-8 -*-
|
||||
"""
|
||||
Classes in this module enhance Matern covariance functions with the
|
||||
Stochastic Differential Equation (SDE) functionality.
|
||||
"""
|
||||
from .standard_periodic import StdPeriodic
|
||||
|
||||
import numpy as np
|
||||
import scipy as sp
|
||||
|
||||
from scipy import special as special
|
||||
|
||||
class sde_StdPeriodic(StdPeriodic):
|
||||
"""
|
||||
|
||||
Class provide extra functionality to transfer this covariance function into
|
||||
SDE form.
|
||||
|
||||
Standard Periodic kernel:
|
||||
|
||||
.. math::
|
||||
|
||||
k(x,y) = \theta_1 \exp \left[ - \frac{1}{2} {}\sum_{i=1}^{input\_dim}
|
||||
\left( \frac{\sin(\frac{\pi}{\lambda_i} (x_i - y_i) )}{l_i} \right)^2 \right] }
|
||||
|
||||
"""
|
||||
def sde_update_gradient_full(self, gradients):
|
||||
"""
|
||||
Update gradient in the order in which parameters are represented in the
|
||||
kernel
|
||||
"""
|
||||
|
||||
self.variance.gradient = gradients[0]
|
||||
self.wavelengths.gradient = gradients[1]
|
||||
self.lengthscales.gradient = gradients[2]
|
||||
|
||||
def sde(self):
|
||||
"""
|
||||
Return the state space representation of the covariance.
|
||||
|
||||
|
||||
! Note: one must constrain lengthscale not to drop below 0.25.
|
||||
After this bessel functions of the first kind grows to very high.
|
||||
|
||||
! Note: one must keep wevelength also not very low. Because then
|
||||
the gradients wrt wavelength become ustable.
|
||||
However this might depend on the data. For test example with
|
||||
300 data points the low limit is 0.15.
|
||||
"""
|
||||
|
||||
# Params to use: (in that order)
|
||||
#self.variance
|
||||
#self.wavelengths
|
||||
#self.lengthscales
|
||||
N = 7 # approximation order
|
||||
|
||||
|
||||
w0 = 2*np.pi/self.wavelengths # frequency
|
||||
lengthscales = 2*self.lengthscales
|
||||
|
||||
[q2,dq2l] = seriescoeff(N,lengthscales,self.variance)
|
||||
# lengthscale is multiplied by 2 because of slightly different
|
||||
# formula for periodic covariance function.
|
||||
# For the same reason:
|
||||
|
||||
dq2l = 2*dq2l
|
||||
|
||||
if np.any( np.isfinite(q2) == False):
|
||||
raise ValueError("SDE periodic covariance error 1")
|
||||
|
||||
if np.any( np.isfinite(dq2l) == False):
|
||||
raise ValueError("SDE periodic covariance error 2")
|
||||
|
||||
F = np.kron(np.diag(range(0,N+1)),np.array( ((0, -w0), (w0, 0)) ) )
|
||||
L = np.eye(2*(N+1))
|
||||
Qc = np.zeros((2*(N+1), 2*(N+1)))
|
||||
P_inf = np.kron(np.diag(q2),np.eye(2))
|
||||
H = np.kron(np.ones((1,N+1)),np.array((1,0)) )
|
||||
P0 = P_inf.copy()
|
||||
|
||||
# Derivatives
|
||||
dF = np.empty((F.shape[0], F.shape[1], 3))
|
||||
dQc = np.empty((Qc.shape[0], Qc.shape[1], 3))
|
||||
dP_inf = np.empty((P_inf.shape[0], P_inf.shape[1], 3))
|
||||
|
||||
# Derivatives wrt self.variance
|
||||
dF[:,:,0] = np.zeros(F.shape)
|
||||
dQc[:,:,0] = np.zeros(Qc.shape)
|
||||
dP_inf[:,:,0] = P_inf / self.variance
|
||||
|
||||
# Derivatives self.wavelengths
|
||||
dF[:,:,1] = np.kron(np.diag(range(0,N+1)),np.array( ((0, w0), (-w0, 0)) ) / self.wavelengths );
|
||||
dQc[:,:,1] = np.zeros(Qc.shape)
|
||||
dP_inf[:,:,1] = np.zeros(P_inf.shape)
|
||||
|
||||
# Derivatives self.lengthscales
|
||||
dF[:,:,2] = np.zeros(F.shape)
|
||||
dQc[:,:,2] = np.zeros(Qc.shape)
|
||||
dP_inf[:,:,2] = np.kron(np.diag(dq2l),np.eye(2))
|
||||
dP0 = dP_inf.copy()
|
||||
|
||||
return (F, L, Qc, H, P_inf, P0, dF, dQc, dP_inf, dP0)
|
||||
|
||||
|
||||
|
||||
|
||||
def seriescoeff(m=6,lengthScale=1.0,magnSigma2=1.0, true_covariance=False):
|
||||
"""
|
||||
Calculate the coefficients q_j^2 for the covariance function
|
||||
approximation:
|
||||
|
||||
k(\tau) = \sum_{j=0}^{+\infty} q_j^2 \cos(j\omega_0 \tau)
|
||||
|
||||
Reference is:
|
||||
|
||||
[1] Arno Solin and Simo Särkkä (2014). Explicit link between periodic
|
||||
covariance functions and state space models. In Proceedings of the
|
||||
Seventeenth International Conference on Artifcial Intelligence and
|
||||
Statistics (AISTATS 2014). JMLR: W&CP, volume 33.
|
||||
|
||||
Note! Only the infinite approximation (through Bessel function)
|
||||
is currently implemented.
|
||||
|
||||
Input:
|
||||
----------------
|
||||
|
||||
m: int
|
||||
Degree of approximation. Default 6.
|
||||
lengthScale: float
|
||||
Length scale parameter in the kerenl
|
||||
magnSigma2:float
|
||||
Multiplier in front of the kernel.
|
||||
|
||||
|
||||
Output:
|
||||
-----------------
|
||||
|
||||
coeffs: array(m+1)
|
||||
Covariance series coefficients
|
||||
|
||||
coeffs_dl: array(m+1)
|
||||
Derivatives of the coefficients with respect to lengthscale.
|
||||
|
||||
"""
|
||||
|
||||
if true_covariance:
|
||||
|
||||
bb = lambda j,m: (1.0 + np.array((j != 0), dtype=np.float64) ) / (2**(j)) *\
|
||||
sp.special.binom(j, sp.floor( (j-m)/2.0 * np.array(m<=j, dtype=np.float64) ))*\
|
||||
np.array(m<=j, dtype=np.float64) *np.array(sp.mod(j-m,2)==0, dtype=np.float64)
|
||||
|
||||
M,J = np.meshgrid(range(0,m+1),range(0,m+1))
|
||||
|
||||
coeffs = bb(J,M) / sp.misc.factorial(J) * sp.exp( -lengthScale**(-2) ) *\
|
||||
(lengthScale**(-2))**J *magnSigma2
|
||||
|
||||
coeffs_dl = np.sum( coeffs*lengthScale**(-3)*(2.0-2.0*J*lengthScale**2),0)
|
||||
|
||||
coeffs = np.sum(coeffs,0)
|
||||
|
||||
else:
|
||||
coeffs = 2*magnSigma2*sp.exp( -lengthScale**(-2) ) * special.iv(range(0,m+1),1.0/lengthScale**(2))
|
||||
if np.any( np.isfinite(coeffs) == False):
|
||||
raise ValueError("sde_standard_periodic: Coefficients are not finite!")
|
||||
#import pdb; pdb.set_trace()
|
||||
coeffs[0] = 0.5*coeffs[0]
|
||||
|
||||
# Derivatives wrt (lengthScale)
|
||||
coeffs_dl = np.zeros(m+1)
|
||||
coeffs_dl[1:] = magnSigma2*lengthScale**(-3) * sp.exp(-lengthScale**(-2))*\
|
||||
(-4*special.iv(range(0,m),lengthScale**(-2)) + 4*(1+np.arange(1,m+1)*lengthScale**(2))*special.iv(range(1,m+1),lengthScale**(-2)) )
|
||||
|
||||
# The first element
|
||||
coeffs_dl[0] = magnSigma2*lengthScale**(-3) * np.exp(-lengthScale**(-2))*\
|
||||
(2*special.iv(0,lengthScale**(-2)) - 2*special.iv(1,lengthScale**(-2)) )
|
||||
|
||||
|
||||
return coeffs, coeffs_dl
|
||||
101
GPy/kern/_src/sde_static.py
Normal file
101
GPy/kern/_src/sde_static.py
Normal file
|
|
@ -0,0 +1,101 @@
|
|||
# -*- coding: utf-8 -*-
|
||||
"""
|
||||
Classes in this module enhance Static covariance functions with the
|
||||
Stochastic Differential Equation (SDE) functionality.
|
||||
"""
|
||||
from .static import White
|
||||
from .static import Bias
|
||||
|
||||
import numpy as np
|
||||
|
||||
class sde_White(White):
|
||||
"""
|
||||
|
||||
Class provide extra functionality to transfer this covariance function into
|
||||
SDE forrm.
|
||||
|
||||
White kernel:
|
||||
|
||||
.. math::
|
||||
|
||||
k(x,y) = \alpha*\delta(x-y)
|
||||
|
||||
"""
|
||||
|
||||
def sde_update_gradient_full(self, gradients):
|
||||
"""
|
||||
Update gradient in the order in which parameters are represented in the
|
||||
kernel
|
||||
"""
|
||||
|
||||
self.variance.gradient = gradients[0]
|
||||
|
||||
def sde(self):
|
||||
"""
|
||||
Return the state space representation of the covariance.
|
||||
"""
|
||||
|
||||
variance = float(self.variance.values)
|
||||
|
||||
F = np.array( ((-np.inf,),) )
|
||||
L = np.array( ((1.0,),) )
|
||||
Qc = np.array( ((variance,),) )
|
||||
H = np.array( ((1.0,),) )
|
||||
|
||||
Pinf = np.array( ((variance,),) )
|
||||
P0 = Pinf.copy()
|
||||
|
||||
dF = np.zeros((1,1,1))
|
||||
dQc = np.zeros((1,1,1))
|
||||
dQc[:,:,0] = np.array( ((1.0,),) )
|
||||
|
||||
dPinf = np.zeros((1,1,1))
|
||||
dPinf[:,:,0] = np.array( ((1.0,),) )
|
||||
dP0 = dPinf.copy()
|
||||
|
||||
return (F, L, Qc, H, Pinf, P0, dF, dQc, dPinf, dP0)
|
||||
|
||||
|
||||
class sde_Bias(Bias):
|
||||
"""
|
||||
|
||||
Class provide extra functionality to transfer this covariance function into
|
||||
SDE forrm.
|
||||
|
||||
Bias kernel:
|
||||
|
||||
.. math::
|
||||
|
||||
k(x,y) = \alpha
|
||||
|
||||
"""
|
||||
def sde_update_gradient_full(self, gradients):
|
||||
"""
|
||||
Update gradient in the order in which parameters are represented in the
|
||||
kernel
|
||||
"""
|
||||
|
||||
self.variance.gradient = gradients[0]
|
||||
|
||||
def sde(self):
|
||||
"""
|
||||
Return the state space representation of the covariance.
|
||||
"""
|
||||
variance = float(self.variance.values)
|
||||
|
||||
F = np.array( ((0.0,),))
|
||||
L = np.array( ((1.0,),))
|
||||
Qc = np.zeros((1,1))
|
||||
H = np.array( ((1.0,),))
|
||||
|
||||
Pinf = np.zeros((1,1))
|
||||
P0 = np.array( ((variance,),) )
|
||||
|
||||
dF = np.zeros((1,1,1))
|
||||
dQc = np.zeros((1,1,1))
|
||||
|
||||
dPinf = np.zeros((1,1,1))
|
||||
dP0 = np.zeros((1,1,1))
|
||||
dP0[:,:,0] = np.array( ((1.0,),) )
|
||||
|
||||
return (F, L, Qc, H, Pinf, P0, dF, dQc, dPinf, dP0)
|
||||
190
GPy/kern/_src/sde_stationary.py
Normal file
190
GPy/kern/_src/sde_stationary.py
Normal file
|
|
@ -0,0 +1,190 @@
|
|||
# -*- coding: utf-8 -*-
|
||||
"""
|
||||
Classes in this module enhance several stationary covariance functions with the
|
||||
Stochastic Differential Equation (SDE) functionality.
|
||||
"""
|
||||
from .rbf import RBF
|
||||
from .stationary import Exponential
|
||||
from .stationary import RatQuad
|
||||
|
||||
import numpy as np
|
||||
import scipy as sp
|
||||
|
||||
class sde_RBF(RBF):
|
||||
"""
|
||||
|
||||
Class provide extra functionality to transfer this covariance function into
|
||||
SDE form.
|
||||
|
||||
Radial Basis Function kernel:
|
||||
|
||||
.. math::
|
||||
|
||||
k(r) = \sigma^2 \exp \\bigg(- \\frac{1}{2} r^2 \\bigg) \\ \\ \\ \\ \text{ where } r = \sqrt{\sum_{i=1}^{input dim} \frac{(x_i-y_i)^2}{\ell_i^2} }
|
||||
|
||||
"""
|
||||
def sde_update_gradient_full(self, gradients):
|
||||
"""
|
||||
Update gradient in the order in which parameters are represented in the
|
||||
kernel
|
||||
"""
|
||||
|
||||
self.variance.gradient = gradients[0]
|
||||
self.lengthscale.gradient = gradients[1]
|
||||
|
||||
def sde(self):
|
||||
"""
|
||||
Return the state space representation of the covariance.
|
||||
"""
|
||||
|
||||
N = 10# approximation order ( number of terms in exponent series expansion)
|
||||
roots_rounding_decimals = 6
|
||||
|
||||
fn = np.math.factorial(N)
|
||||
|
||||
kappa = 1.0/2.0/self.lengthscale**2
|
||||
|
||||
Qc = np.array((self.variance*np.sqrt(np.pi/kappa)*fn*(4*kappa)**N,),)
|
||||
|
||||
pp = np.zeros((2*N+1,)) # array of polynomial coefficients from higher power to lower
|
||||
|
||||
for n in range(0, N+1): # (2N+1) - number of polynomial coefficients
|
||||
pp[2*(N-n)] = fn*(4.0*kappa)**(N-n)/np.math.factorial(n)*(-1)**n
|
||||
|
||||
pp = sp.poly1d(pp)
|
||||
roots = sp.roots(pp)
|
||||
|
||||
neg_real_part_roots = roots[np.round(np.real(roots) ,roots_rounding_decimals) < 0]
|
||||
aa = sp.poly1d(neg_real_part_roots, r=True).coeffs
|
||||
|
||||
F = np.diag(np.ones((N-1,)),1)
|
||||
F[-1,:] = -aa[-1:0:-1]
|
||||
|
||||
L= np.zeros((N,1))
|
||||
L[N-1,0] = 1
|
||||
|
||||
H = np.zeros((1,N))
|
||||
H[0,0] = 1
|
||||
|
||||
# Infinite covariance:
|
||||
Pinf = sp.linalg.solve_lyapunov(F, -np.dot(L,np.dot( Qc[0,0],L.T)))
|
||||
Pinf = 0.5*(Pinf + Pinf.T)
|
||||
# Allocating space for derivatives
|
||||
dF = np.empty([F.shape[0],F.shape[1],2])
|
||||
dQc = np.empty([Qc.shape[0],Qc.shape[1],2])
|
||||
dPinf = np.empty([Pinf.shape[0],Pinf.shape[1],2])
|
||||
|
||||
# Derivatives:
|
||||
dFvariance = np.zeros(F.shape)
|
||||
dFlengthscale = np.zeros(F.shape)
|
||||
dFlengthscale[-1,:] = -aa[-1:0:-1]/self.lengthscale * np.arange(-N,0,1)
|
||||
|
||||
dQcvariance = Qc/self.variance
|
||||
dQclengthscale = np.array(((self.variance*np.sqrt(2*np.pi)*fn*2**N*self.lengthscale**(-2*N)*(1-2*N,),)))
|
||||
|
||||
dPinf_variance = Pinf/self.variance
|
||||
|
||||
lp = Pinf.shape[0]
|
||||
coeff = np.arange(1,lp+1).reshape(lp,1) + np.arange(1,lp+1).reshape(1,lp) - 2
|
||||
coeff[np.mod(coeff,2) != 0] = 0
|
||||
dPinf_lengthscale = -1/self.lengthscale*Pinf*coeff
|
||||
|
||||
dF[:,:,0] = dFvariance
|
||||
dF[:,:,1] = dFlengthscale
|
||||
dQc[:,:,0] = dQcvariance
|
||||
dQc[:,:,1] = dQclengthscale
|
||||
dPinf[:,:,0] = dPinf_variance
|
||||
dPinf[:,:,1] = dPinf_lengthscale
|
||||
|
||||
P0 = Pinf.copy()
|
||||
dP0 = dPinf.copy()
|
||||
|
||||
# Benefits of this are not very sound. Helps only in one case:
|
||||
# SVD Kalman + RBF kernel
|
||||
import GPy.models.state_space_main as ssm
|
||||
(F, L, Qc, H, Pinf, P0, dF, dQc, dPinf,dP0, T) = ssm.balance_ss_model(F, L, Qc, H, Pinf, P0, dF, dQc, dPinf, dP0 )
|
||||
|
||||
return (F, L, Qc, H, Pinf, P0, dF, dQc, dPinf, dP0)
|
||||
|
||||
class sde_Exponential(Exponential):
|
||||
"""
|
||||
|
||||
Class provide extra functionality to transfer this covariance function into
|
||||
SDE form.
|
||||
|
||||
Exponential kernel:
|
||||
|
||||
.. math::
|
||||
|
||||
k(r) = \sigma^2 \exp \\bigg(- \\frac{1}{2} r \\bigg) \\ \\ \\ \\ \text{ where } r = \sqrt{\sum_{i=1}^{input dim} \frac{(x_i-y_i)^2}{\ell_i^2} }
|
||||
|
||||
"""
|
||||
|
||||
def sde_update_gradient_full(self, gradients):
|
||||
"""
|
||||
Update gradient in the order in which parameters are represented in the
|
||||
kernel
|
||||
"""
|
||||
|
||||
self.variance.gradient = gradients[0]
|
||||
self.lengthscale.gradient = gradients[1]
|
||||
|
||||
def sde(self):
|
||||
"""
|
||||
Return the state space representation of the covariance.
|
||||
"""
|
||||
variance = float(self.variance.values)
|
||||
lengthscale = float(self.lengthscale)
|
||||
|
||||
F = np.array(((-1.0/lengthscale,),))
|
||||
L = np.array(((1.0,),))
|
||||
Qc = np.array( ((2.0*variance/lengthscale,),) )
|
||||
H = np.array(((1.0,),))
|
||||
Pinf = np.array(((variance,),))
|
||||
P0 = Pinf.copy()
|
||||
|
||||
dF = np.zeros((1,1,2));
|
||||
dQc = np.zeros((1,1,2));
|
||||
dPinf = np.zeros((1,1,2));
|
||||
|
||||
dF[:,:,0] = 0.0
|
||||
dF[:,:,1] = 1.0/lengthscale**2
|
||||
|
||||
dQc[:,:,0] = 2.0/lengthscale
|
||||
dQc[:,:,1] = -2.0*variance/lengthscale**2
|
||||
|
||||
dPinf[:,:,0] = 1.0
|
||||
dPinf[:,:,1] = 0.0
|
||||
|
||||
dP0 = dPinf.copy()
|
||||
|
||||
return (F, L, Qc, H, Pinf, P0, dF, dQc, dPinf, dP0)
|
||||
|
||||
class sde_RatQuad(RatQuad):
|
||||
"""
|
||||
|
||||
Class provide extra functionality to transfer this covariance function into
|
||||
SDE form.
|
||||
|
||||
Rational Quadratic kernel:
|
||||
|
||||
.. math::
|
||||
|
||||
k(r) = \sigma^2 \\bigg( 1 + \\frac{r^2}{2} \\bigg)^{- \alpha} \\ \\ \\ \\ \text{ where } r = \sqrt{\sum_{i=1}^{input dim} \frac{(x_i-y_i)^2}{\ell_i^2} }
|
||||
|
||||
"""
|
||||
|
||||
def sde(self):
|
||||
"""
|
||||
Return the state space representation of the covariance.
|
||||
"""
|
||||
|
||||
assert False, 'Not Implemented'
|
||||
|
||||
# Params to use:
|
||||
|
||||
# self.lengthscale
|
||||
# self.variance
|
||||
#self.power
|
||||
|
||||
#return (F, L, Qc, H, Pinf, dF, dQc, dPinf)
|
||||
|
|
@ -263,4 +263,94 @@ class Add(CombinationKernel):
|
|||
i_s[k._all_dims_active] += k.input_sensitivity(summarize)
|
||||
return i_s
|
||||
else:
|
||||
|
||||
return super(Add, self).input_sensitivity(summarize)
|
||||
|
||||
def sde_update_gradient_full(self, gradients):
|
||||
"""
|
||||
Update gradient in the order in which parameters are represented in the
|
||||
kernel
|
||||
"""
|
||||
part_start_param_index = 0
|
||||
for p in self.parts:
|
||||
if not p.is_fixed:
|
||||
part_param_num = len(p.param_array) # number of parameters in the part
|
||||
p.sde_update_gradient_full(gradients[part_start_param_index:(part_start_param_index+part_param_num)])
|
||||
part_start_param_index += part_param_num
|
||||
|
||||
def sde(self):
|
||||
"""
|
||||
Support adding kernels for sde representation
|
||||
"""
|
||||
|
||||
import scipy.linalg as la
|
||||
|
||||
F = None
|
||||
L = None
|
||||
Qc = None
|
||||
H = None
|
||||
Pinf = None
|
||||
P0 = None
|
||||
dF = None
|
||||
dQc = None
|
||||
dPinf = None
|
||||
dP0 = None
|
||||
n = 0
|
||||
nq = 0
|
||||
nd = 0
|
||||
|
||||
# Assign models
|
||||
for p in self.parts:
|
||||
(Ft,Lt,Qct,Ht,Pinft,P0t,dFt,dQct,dPinft,dP0t) = p.sde()
|
||||
F = la.block_diag(F,Ft) if (F is not None) else Ft
|
||||
L = la.block_diag(L,Lt) if (L is not None) else Lt
|
||||
Qc = la.block_diag(Qc,Qct) if (Qc is not None) else Qct
|
||||
H = np.hstack((H,Ht)) if (H is not None) else Ht
|
||||
|
||||
Pinf = la.block_diag(Pinf,Pinft) if (Pinf is not None) else Pinft
|
||||
P0 = la.block_diag(P0,P0t) if (P0 is not None) else P0t
|
||||
|
||||
if dF is not None:
|
||||
dF = np.pad(dF,((0,dFt.shape[0]),(0,dFt.shape[1]),(0,dFt.shape[2])),
|
||||
'constant', constant_values=0)
|
||||
dF[-dFt.shape[0]:,-dFt.shape[1]:,-dFt.shape[2]:] = dFt
|
||||
else:
|
||||
dF = dFt
|
||||
|
||||
if dQc is not None:
|
||||
dQc = np.pad(dQc,((0,dQct.shape[0]),(0,dQct.shape[1]),(0,dQct.shape[2])),
|
||||
'constant', constant_values=0)
|
||||
dQc[-dQct.shape[0]:,-dQct.shape[1]:,-dQct.shape[2]:] = dQct
|
||||
else:
|
||||
dQc = dQct
|
||||
|
||||
if dPinf is not None:
|
||||
dPinf = np.pad(dPinf,((0,dPinft.shape[0]),(0,dPinft.shape[1]),(0,dPinft.shape[2])),
|
||||
'constant', constant_values=0)
|
||||
dPinf[-dPinft.shape[0]:,-dPinft.shape[1]:,-dPinft.shape[2]:] = dPinft
|
||||
else:
|
||||
dPinf = dPinft
|
||||
|
||||
if dP0 is not None:
|
||||
dP0 = np.pad(dP0,((0,dP0t.shape[0]),(0,dP0t.shape[1]),(0,dP0t.shape[2])),
|
||||
'constant', constant_values=0)
|
||||
dP0[-dP0t.shape[0]:,-dP0t.shape[1]:,-dP0t.shape[2]:] = dP0t
|
||||
else:
|
||||
dP0 = dP0t
|
||||
|
||||
n += Ft.shape[0]
|
||||
nq += Qct.shape[0]
|
||||
nd += dFt.shape[2]
|
||||
|
||||
assert (F.shape[0] == n and F.shape[1]==n), "SDE add: Check of F Dimensions failed"
|
||||
assert (L.shape[0] == n and L.shape[1]==nq), "SDE add: Check of L Dimensions failed"
|
||||
assert (Qc.shape[0] == nq and Qc.shape[1]==nq), "SDE add: Check of Qc Dimensions failed"
|
||||
assert (H.shape[0] == 1 and H.shape[1]==n), "SDE add: Check of H Dimensions failed"
|
||||
assert (Pinf.shape[0] == n and Pinf.shape[1]==n), "SDE add: Check of Pinf Dimensions failed"
|
||||
assert (P0.shape[0] == n and P0.shape[1]==n), "SDE add: Check of P0 Dimensions failed"
|
||||
assert (dF.shape[0] == n and dF.shape[1]==n and dF.shape[2]==nd), "SDE add: Check of dF Dimensions failed"
|
||||
assert (dQc.shape[0] == nq and dQc.shape[1]==nq and dQc.shape[2]==nd), "SDE add: Check of dQc Dimensions failed"
|
||||
assert (dPinf.shape[0] == n and dPinf.shape[1]==n and dPinf.shape[2]==nd), "SDE add: Check of dPinf Dimensions failed"
|
||||
assert (dP0.shape[0] == n and dP0.shape[1]==n and dP0.shape[2]==nd), "SDE add: Check of dP0 Dimensions failed"
|
||||
|
||||
return (F,L,Qc,H,Pinf,P0,dF,dQc,dPinf,dP0)
|
||||
|
|
|
|||
|
|
@ -305,7 +305,6 @@ class Kern(Parameterized):
|
|||
def _check_active_dims(self, X):
|
||||
assert X.shape[1] >= len(self._all_dims_active), "At least {} dimensional X needed, X.shape={!s}".format(len(self._all_dims_active), X.shape)
|
||||
|
||||
|
||||
class CombinationKernel(Kern):
|
||||
"""
|
||||
Abstract super class for combination kernels.
|
||||
|
|
|
|||
|
|
@ -105,3 +105,114 @@ class Prod(CombinationKernel):
|
|||
return i_s
|
||||
else:
|
||||
return super(Prod, self).input_sensitivity(summarize)
|
||||
|
||||
def sde_update_gradient_full(self, gradients):
|
||||
"""
|
||||
Update gradient in the order in which parameters are represented in the
|
||||
kernel
|
||||
"""
|
||||
part_start_param_index = 0
|
||||
for p in self.parts:
|
||||
if not p.is_fixed:
|
||||
part_param_num = len(p.param_array) # number of parameters in the part
|
||||
p.sde_update_gradient_full(gradients[part_start_param_index:(part_start_param_index+part_param_num)])
|
||||
part_start_param_index += part_param_num
|
||||
|
||||
def sde(self):
|
||||
"""
|
||||
"""
|
||||
F = np.array((0,), ndmin=2)
|
||||
L = np.array((1,), ndmin=2)
|
||||
Qc = np.array((1,), ndmin=2)
|
||||
H = np.array((1,), ndmin=2)
|
||||
Pinf = np.array((1,), ndmin=2)
|
||||
P0 = np.array((1,), ndmin=2)
|
||||
dF = None
|
||||
dQc = None
|
||||
dPinf = None
|
||||
dP0 = None
|
||||
|
||||
# Assign models
|
||||
for p in self.parts:
|
||||
(Ft,Lt,Qct,Ht,P_inft, P0t, dFt,dQct,dP_inft,dP0t) = p.sde()
|
||||
|
||||
# check derivative dimensions ->
|
||||
number_of_parameters = len(p.param_array)
|
||||
assert dFt.shape[2] == number_of_parameters, "Dynamic matrix derivative shape is wrong"
|
||||
assert dQct.shape[2] == number_of_parameters, "Diffusion matrix derivative shape is wrong"
|
||||
assert dP_inft.shape[2] == number_of_parameters, "Infinite covariance matrix derivative shape is wrong"
|
||||
# check derivative dimensions <-
|
||||
|
||||
# exception for periodic kernel
|
||||
if (p.name == 'std_periodic'):
|
||||
Qct = P_inft
|
||||
dQct = dP_inft
|
||||
|
||||
dF = dkron(F,dF,Ft,dFt,'sum')
|
||||
dQc = dkron(Qc,dQc,Qct,dQct,'prod')
|
||||
dPinf = dkron(Pinf,dPinf,P_inft,dP_inft,'prod')
|
||||
dP0 = dkron(P0,dP0,P0t,dP0t,'prod')
|
||||
|
||||
F = np.kron(F,np.eye(Ft.shape[0])) + np.kron(np.eye(F.shape[0]),Ft)
|
||||
L = np.kron(L,Lt)
|
||||
Qc = np.kron(Qc,Qct)
|
||||
Pinf = np.kron(Pinf,P_inft)
|
||||
P0 = np.kron(P0,P_inft)
|
||||
H = np.kron(H,Ht)
|
||||
|
||||
return (F,L,Qc,H,Pinf,P0,dF,dQc,dPinf,dP0)
|
||||
|
||||
def dkron(A,dA,B,dB, operation='prod'):
|
||||
"""
|
||||
Function computes the derivative of Kronecker product A*B
|
||||
(or Kronecker sum A+B).
|
||||
|
||||
Input:
|
||||
-----------------------
|
||||
|
||||
A: 2D matrix
|
||||
Some matrix
|
||||
dA: 3D (or 2D matrix)
|
||||
Derivarives of A
|
||||
B: 2D matrix
|
||||
Some matrix
|
||||
dB: 3D (or 2D matrix)
|
||||
Derivarives of B
|
||||
|
||||
operation: str 'prod' or 'sum'
|
||||
Which operation is considered. If the operation is 'sum' it is assumed
|
||||
that A and are square matrices.s
|
||||
|
||||
Output:
|
||||
dC: 3D matrix
|
||||
Derivative of Kronecker product A*B (or Kronecker sum A+B)
|
||||
"""
|
||||
|
||||
if dA is None:
|
||||
dA_param_num = 0
|
||||
dA = np.zeros((A.shape[0], A.shape[1],1))
|
||||
else:
|
||||
dA_param_num = dA.shape[2]
|
||||
|
||||
if dB is None:
|
||||
dB_param_num = 0
|
||||
dB = np.zeros((B.shape[0], B.shape[1],1))
|
||||
else:
|
||||
dB_param_num = dB.shape[2]
|
||||
|
||||
# Space allocation for derivative matrix
|
||||
dC = np.zeros((A.shape[0]*B.shape[0], A.shape[1]*B.shape[1], dA_param_num + dB_param_num))
|
||||
|
||||
for k in range(dA_param_num):
|
||||
if operation == 'prod':
|
||||
dC[:,:,k] = np.kron(dA[:,:,k],B);
|
||||
else:
|
||||
dC[:,:,k] = np.kron(dA[:,:,k],np.eye( B.shape[0] ))
|
||||
|
||||
for k in range(dB_param_num):
|
||||
if operation == 'prod':
|
||||
dC[:,:,dA_param_num+k] = np.kron(A,dB[:,:,k])
|
||||
else:
|
||||
dC[:,:,dA_param_num+k] = np.kron(np.eye( A.shape[0] ),dB[:,:,k])
|
||||
|
||||
return dC
|
||||
|
|
|
|||
59
GPy/kern/src/sde_brownian.py
Normal file
59
GPy/kern/src/sde_brownian.py
Normal file
|
|
@ -0,0 +1,59 @@
|
|||
# -*- coding: utf-8 -*-
|
||||
# Copyright (c) 2015, Alex Grigorevskiy, Arno Solin
|
||||
# Licensed under the BSD 3-clause license (see LICENSE.txt)
|
||||
"""
|
||||
Classes in this module enhance Brownian motion covariance function with the
|
||||
Stochastic Differential Equation (SDE) functionality.
|
||||
"""
|
||||
|
||||
from .brownian import Brownian
|
||||
|
||||
import numpy as np
|
||||
|
||||
class sde_Brownian(Brownian):
|
||||
"""
|
||||
|
||||
Class provide extra functionality to transfer this covariance function into
|
||||
SDE form.
|
||||
|
||||
Linear kernel:
|
||||
|
||||
.. math::
|
||||
|
||||
k(x,y) = \sigma^2 min(x,y)
|
||||
|
||||
"""
|
||||
|
||||
def sde_update_gradient_full(self, gradients):
|
||||
"""
|
||||
Update gradient in the order in which parameters are represented in the
|
||||
kernel
|
||||
"""
|
||||
|
||||
self.variance.gradient = gradients[0]
|
||||
|
||||
def sde(self):
|
||||
"""
|
||||
Return the state space representation of the covariance.
|
||||
"""
|
||||
|
||||
variance = float(self.variance.values) # this is initial variancve in Bayesian linear regression
|
||||
|
||||
F = np.array( ((0,1.0),(0,0) ))
|
||||
L = np.array( ((1.0,),(0,)) )
|
||||
Qc = np.array( ((variance,),) )
|
||||
H = np.array( ((1.0,0),) )
|
||||
|
||||
Pinf = np.array( ( (0, -0.5*variance ), (-0.5*variance, 0) ) )
|
||||
#P0 = Pinf.copy()
|
||||
P0 = np.zeros((2,2))
|
||||
#Pinf = np.array( ( (t0, 1.0), (1.0, 1.0/t0) ) ) * variance
|
||||
dF = np.zeros((2,2,1))
|
||||
dQc = np.ones( (1,1,1) )
|
||||
|
||||
dPinf = np.zeros((2,2,1))
|
||||
dPinf[:,:,0] = np.array( ( (0, -0.5), (-0.5, 0) ) )
|
||||
#dP0 = dPinf.copy()
|
||||
dP0 = np.zeros((2,2,1))
|
||||
|
||||
return (F, L, Qc, H, Pinf, P0, dF, dQc, dPinf, dP0)
|
||||
66
GPy/kern/src/sde_linear.py
Normal file
66
GPy/kern/src/sde_linear.py
Normal file
|
|
@ -0,0 +1,66 @@
|
|||
# -*- coding: utf-8 -*-
|
||||
# Copyright (c) 2015, Alex Grigorevskiy, Arno Solin
|
||||
# Licensed under the BSD 3-clause license (see LICENSE.txt)
|
||||
"""
|
||||
Classes in this module enhance Linear covariance function with the
|
||||
Stochastic Differential Equation (SDE) functionality.
|
||||
"""
|
||||
from .linear import Linear
|
||||
|
||||
import numpy as np
|
||||
|
||||
class sde_Linear(Linear):
|
||||
"""
|
||||
|
||||
Class provide extra functionality to transfer this covariance function into
|
||||
SDE form.
|
||||
|
||||
Linear kernel:
|
||||
|
||||
.. math::
|
||||
|
||||
k(x,y) = \sum_{i=1}^{input dim} \sigma^2_i x_iy_i
|
||||
|
||||
"""
|
||||
def __init__(self, input_dim, X, variances=None, ARD=False, active_dims=None, name='linear'):
|
||||
"""
|
||||
Modify the init method, because one extra parameter is required. X - points
|
||||
on the X axis.
|
||||
"""
|
||||
|
||||
super(sde_Linear, self).__init__(input_dim, variances, ARD, active_dims, name)
|
||||
|
||||
self.t0 = np.min(X)
|
||||
|
||||
|
||||
def sde_update_gradient_full(self, gradients):
|
||||
"""
|
||||
Update gradient in the order in which parameters are represented in the
|
||||
kernel
|
||||
"""
|
||||
|
||||
self.variances.gradient = gradients[0]
|
||||
|
||||
def sde(self):
|
||||
"""
|
||||
Return the state space representation of the covariance.
|
||||
"""
|
||||
|
||||
variance = float(self.variances.values) # this is initial variancve in Bayesian linear regression
|
||||
t0 = float(self.t0)
|
||||
|
||||
F = np.array( ((0,1.0),(0,0) ))
|
||||
L = np.array( ((0,),(1.0,)) )
|
||||
Qc = np.zeros((1,1))
|
||||
H = np.array( ((1.0,0),) )
|
||||
|
||||
Pinf = np.zeros((2,2))
|
||||
P0 = np.array( ( (t0**2, t0), (t0, 1) ) ) * variance
|
||||
dF = np.zeros((2,2,1))
|
||||
dQc = np.zeros( (1,1,1) )
|
||||
|
||||
dPinf = np.zeros((2,2,1))
|
||||
dP0 = np.zeros((2,2,1))
|
||||
dP0[:,:,0] = P0 / variance
|
||||
|
||||
return (F, L, Qc, H, Pinf, P0, dF, dQc, dPinf, dP0)
|
||||
137
GPy/kern/src/sde_matern.py
Normal file
137
GPy/kern/src/sde_matern.py
Normal file
|
|
@ -0,0 +1,137 @@
|
|||
# -*- coding: utf-8 -*-
|
||||
# Copyright (c) 2015, Alex Grigorevskiy, Arno Solin
|
||||
# Licensed under the BSD 3-clause license (see LICENSE.txt)
|
||||
"""
|
||||
Classes in this module enhance Matern covariance functions with the
|
||||
Stochastic Differential Equation (SDE) functionality.
|
||||
"""
|
||||
from .stationary import Matern32
|
||||
from .stationary import Matern52
|
||||
import numpy as np
|
||||
|
||||
class sde_Matern32(Matern32):
|
||||
"""
|
||||
|
||||
Class provide extra functionality to transfer this covariance function into
|
||||
SDE forrm.
|
||||
|
||||
Matern 3/2 kernel:
|
||||
|
||||
.. math::
|
||||
|
||||
k(r) = \sigma^2 (1 + \sqrt{3} r) \exp(- \sqrt{3} r) \\ \\ \\ \\ \text{ where } r = \sqrt{\sum_{i=1}^{input dim} \frac{(x_i-y_i)^2}{\ell_i^2} }
|
||||
|
||||
"""
|
||||
def sde_update_gradient_full(self, gradients):
|
||||
"""
|
||||
Update gradient in the order in which parameters are represented in the
|
||||
kernel
|
||||
"""
|
||||
|
||||
self.variance.gradient = gradients[0]
|
||||
self.lengthscale.gradient = gradients[1]
|
||||
|
||||
def sde(self):
|
||||
"""
|
||||
Return the state space representation of the covariance.
|
||||
"""
|
||||
|
||||
variance = float(self.variance.values)
|
||||
lengthscale = float(self.lengthscale.values)
|
||||
|
||||
foo = np.sqrt(3.)/lengthscale
|
||||
F = np.array(((0, 1.0), (-foo**2, -2*foo)))
|
||||
L = np.array(( (0,), (1.0,) ))
|
||||
Qc = np.array(((12.*np.sqrt(3) / lengthscale**3 * variance,),))
|
||||
H = np.array(((1.0, 0),))
|
||||
Pinf = np.array(((variance, 0.0), (0.0, 3.*variance/(lengthscale**2))))
|
||||
P0 = Pinf.copy()
|
||||
|
||||
# Allocate space for the derivatives
|
||||
dF = np.empty([F.shape[0],F.shape[1],2])
|
||||
dQc = np.empty([Qc.shape[0],Qc.shape[1],2])
|
||||
dPinf = np.empty([Pinf.shape[0],Pinf.shape[1],2])
|
||||
# The partial derivatives
|
||||
dFvariance = np.zeros((2,2))
|
||||
dFlengthscale = np.array(((0,0), (6./lengthscale**3,2*np.sqrt(3)/lengthscale**2)))
|
||||
dQcvariance = np.array((12.*np.sqrt(3)/lengthscale**3))
|
||||
dQclengthscale = np.array((-3*12*np.sqrt(3)/lengthscale**4*variance))
|
||||
dPinfvariance = np.array(((1,0),(0,3./lengthscale**2)))
|
||||
dPinflengthscale = np.array(((0,0), (0,-6*variance/lengthscale**3)))
|
||||
# Combine the derivatives
|
||||
dF[:,:,0] = dFvariance
|
||||
dF[:,:,1] = dFlengthscale
|
||||
dQc[:,:,0] = dQcvariance
|
||||
dQc[:,:,1] = dQclengthscale
|
||||
dPinf[:,:,0] = dPinfvariance
|
||||
dPinf[:,:,1] = dPinflengthscale
|
||||
dP0 = dPinf.copy()
|
||||
|
||||
return (F, L, Qc, H, Pinf, P0, dF, dQc, dPinf, dP0)
|
||||
|
||||
class sde_Matern52(Matern52):
|
||||
"""
|
||||
|
||||
Class provide extra functionality to transfer this covariance function into
|
||||
SDE forrm.
|
||||
|
||||
Matern 5/2 kernel:
|
||||
|
||||
.. math::
|
||||
|
||||
k(r) = \sigma^2 (1 + \sqrt{5} r + \frac{5}{3}r^2) \exp(- \sqrt{5} r) \\ \\ \\ \\ \text{ where } r = \sqrt{\sum_{i=1}^{input dim} \frac{(x_i-y_i)^2}{\ell_i^2} }
|
||||
|
||||
"""
|
||||
def sde_update_gradient_full(self, gradients):
|
||||
"""
|
||||
Update gradient in the order in which parameters are represented in the
|
||||
kernel
|
||||
"""
|
||||
|
||||
self.variance.gradient = gradients[0]
|
||||
self.lengthscale.gradient = gradients[1]
|
||||
|
||||
def sde(self):
|
||||
"""
|
||||
Return the state space representation of the covariance.
|
||||
"""
|
||||
|
||||
variance = float(self.variance.values)
|
||||
lengthscale = float(self.lengthscale.values)
|
||||
|
||||
lamda = np.sqrt(5.0)/lengthscale
|
||||
kappa = 5.0/3.0*variance/lengthscale**2
|
||||
|
||||
F = np.array(((0, 1,0), (0, 0, 1), (-lamda**3, -3.0*lamda**2, -3*lamda)))
|
||||
L = np.array(((0,),(0,),(1,)))
|
||||
Qc = np.array((((variance*400.0*np.sqrt(5.0)/3.0/lengthscale**5),),))
|
||||
H = np.array(((1,0,0),))
|
||||
|
||||
Pinf = np.array(((variance,0,-kappa), (0, kappa, 0), (-kappa, 0, 25.0*variance/lengthscale**4)))
|
||||
P0 = Pinf.copy()
|
||||
# Allocate space for the derivatives
|
||||
dF = np.empty((3,3,2))
|
||||
dQc = np.empty((1,1,2))
|
||||
dPinf = np.empty((3,3,2))
|
||||
|
||||
# The partial derivatives
|
||||
dFvariance = np.zeros((3,3))
|
||||
dFlengthscale = np.array(((0,0,0),(0,0,0),(15.0*np.sqrt(5.0)/lengthscale**4,
|
||||
30.0/lengthscale**3, 3*np.sqrt(5.0)/lengthscale**2)))
|
||||
dQcvariance = np.array((((400*np.sqrt(5)/3/lengthscale**5,),)))
|
||||
dQclengthscale = np.array((((-variance*2000*np.sqrt(5)/3/lengthscale**6,),)))
|
||||
|
||||
dPinf_variance = Pinf/variance
|
||||
kappa2 = -2.0*kappa/lengthscale
|
||||
dPinf_lengthscale = np.array(((0,0,-kappa2),(0,kappa2,0),(-kappa2,
|
||||
0,-100*variance/lengthscale**5)))
|
||||
# Combine the derivatives
|
||||
dF[:,:,0] = dFvariance
|
||||
dF[:,:,1] = dFlengthscale
|
||||
dQc[:,:,0] = dQcvariance
|
||||
dQc[:,:,1] = dQclengthscale
|
||||
dPinf[:,:,0] = dPinf_variance
|
||||
dPinf[:,:,1] = dPinf_lengthscale
|
||||
dP0 = dPinf.copy()
|
||||
|
||||
return (F, L, Qc, H, Pinf, P0, dF, dQc, dPinf, dP0)
|
||||
180
GPy/kern/src/sde_standard_periodic.py
Normal file
180
GPy/kern/src/sde_standard_periodic.py
Normal file
|
|
@ -0,0 +1,180 @@
|
|||
# -*- coding: utf-8 -*-
|
||||
# Copyright (c) 2015, Alex Grigorevskiy, Arno Solin
|
||||
# Licensed under the BSD 3-clause license (see LICENSE.txt)
|
||||
"""
|
||||
Classes in this module enhance Matern covariance functions with the
|
||||
Stochastic Differential Equation (SDE) functionality.
|
||||
"""
|
||||
from .standard_periodic import StdPeriodic
|
||||
|
||||
import numpy as np
|
||||
import scipy as sp
|
||||
|
||||
from scipy import special as special
|
||||
|
||||
class sde_StdPeriodic(StdPeriodic):
|
||||
"""
|
||||
|
||||
Class provide extra functionality to transfer this covariance function into
|
||||
SDE form.
|
||||
|
||||
Standard Periodic kernel:
|
||||
|
||||
.. math::
|
||||
|
||||
k(x,y) = \theta_1 \exp \left[ - \frac{1}{2} {}\sum_{i=1}^{input\_dim}
|
||||
\left( \frac{\sin(\frac{\pi}{\lambda_i} (x_i - y_i) )}{l_i} \right)^2 \right] }
|
||||
|
||||
"""
|
||||
def sde_update_gradient_full(self, gradients):
|
||||
"""
|
||||
Update gradient in the order in which parameters are represented in the
|
||||
kernel
|
||||
"""
|
||||
|
||||
self.variance.gradient = gradients[0]
|
||||
self.period.gradient = gradients[1]
|
||||
self.lengthscale.gradient = gradients[2]
|
||||
|
||||
def sde(self):
|
||||
"""
|
||||
Return the state space representation of the covariance.
|
||||
|
||||
|
||||
! Note: one must constrain lengthscale not to drop below 0.25.
|
||||
After this bessel functions of the first kind grows to very high.
|
||||
|
||||
! Note: one must keep wevelength also not very low. Because then
|
||||
the gradients wrt wavelength become ustable.
|
||||
However this might depend on the data. For test example with
|
||||
300 data points the low limit is 0.15.
|
||||
"""
|
||||
|
||||
# Params to use: (in that order)
|
||||
#self.variance
|
||||
#self.period
|
||||
#self.lengthscale
|
||||
N = 7 # approximation order
|
||||
|
||||
|
||||
w0 = 2*np.pi/self.period # frequency
|
||||
lengthscale = 2*self.lengthscale
|
||||
|
||||
[q2,dq2l] = seriescoeff(N,lengthscale,self.variance)
|
||||
# lengthscale is multiplied by 2 because of slightly different
|
||||
# formula for periodic covariance function.
|
||||
# For the same reason:
|
||||
|
||||
dq2l = 2*dq2l
|
||||
|
||||
if np.any( np.isfinite(q2) == False):
|
||||
raise ValueError("SDE periodic covariance error 1")
|
||||
|
||||
if np.any( np.isfinite(dq2l) == False):
|
||||
raise ValueError("SDE periodic covariance error 2")
|
||||
|
||||
F = np.kron(np.diag(range(0,N+1)),np.array( ((0, -w0), (w0, 0)) ) )
|
||||
L = np.eye(2*(N+1))
|
||||
Qc = np.zeros((2*(N+1), 2*(N+1)))
|
||||
P_inf = np.kron(np.diag(q2),np.eye(2))
|
||||
H = np.kron(np.ones((1,N+1)),np.array((1,0)) )
|
||||
P0 = P_inf.copy()
|
||||
|
||||
# Derivatives
|
||||
dF = np.empty((F.shape[0], F.shape[1], 3))
|
||||
dQc = np.empty((Qc.shape[0], Qc.shape[1], 3))
|
||||
dP_inf = np.empty((P_inf.shape[0], P_inf.shape[1], 3))
|
||||
|
||||
# Derivatives wrt self.variance
|
||||
dF[:,:,0] = np.zeros(F.shape)
|
||||
dQc[:,:,0] = np.zeros(Qc.shape)
|
||||
dP_inf[:,:,0] = P_inf / self.variance
|
||||
|
||||
# Derivatives self.period
|
||||
dF[:,:,1] = np.kron(np.diag(range(0,N+1)),np.array( ((0, w0), (-w0, 0)) ) / self.period );
|
||||
dQc[:,:,1] = np.zeros(Qc.shape)
|
||||
dP_inf[:,:,1] = np.zeros(P_inf.shape)
|
||||
|
||||
# Derivatives self.lengthscales
|
||||
dF[:,:,2] = np.zeros(F.shape)
|
||||
dQc[:,:,2] = np.zeros(Qc.shape)
|
||||
dP_inf[:,:,2] = np.kron(np.diag(dq2l),np.eye(2))
|
||||
dP0 = dP_inf.copy()
|
||||
|
||||
return (F, L, Qc, H, P_inf, P0, dF, dQc, dP_inf, dP0)
|
||||
|
||||
|
||||
|
||||
|
||||
def seriescoeff(m=6,lengthScale=1.0,magnSigma2=1.0, true_covariance=False):
|
||||
"""
|
||||
Calculate the coefficients q_j^2 for the covariance function
|
||||
approximation:
|
||||
|
||||
k(\tau) = \sum_{j=0}^{+\infty} q_j^2 \cos(j\omega_0 \tau)
|
||||
|
||||
Reference is:
|
||||
|
||||
[1] Arno Solin and Simo Särkkä (2014). Explicit link between periodic
|
||||
covariance functions and state space models. In Proceedings of the
|
||||
Seventeenth International Conference on Artifcial Intelligence and
|
||||
Statistics (AISTATS 2014). JMLR: W&CP, volume 33.
|
||||
|
||||
Note! Only the infinite approximation (through Bessel function)
|
||||
is currently implemented.
|
||||
|
||||
Input:
|
||||
----------------
|
||||
|
||||
m: int
|
||||
Degree of approximation. Default 6.
|
||||
lengthScale: float
|
||||
Length scale parameter in the kerenl
|
||||
magnSigma2:float
|
||||
Multiplier in front of the kernel.
|
||||
|
||||
|
||||
Output:
|
||||
-----------------
|
||||
|
||||
coeffs: array(m+1)
|
||||
Covariance series coefficients
|
||||
|
||||
coeffs_dl: array(m+1)
|
||||
Derivatives of the coefficients with respect to lengthscale.
|
||||
|
||||
"""
|
||||
|
||||
if true_covariance:
|
||||
|
||||
bb = lambda j,m: (1.0 + np.array((j != 0), dtype=np.float64) ) / (2**(j)) *\
|
||||
sp.special.binom(j, sp.floor( (j-m)/2.0 * np.array(m<=j, dtype=np.float64) ))*\
|
||||
np.array(m<=j, dtype=np.float64) *np.array(sp.mod(j-m,2)==0, dtype=np.float64)
|
||||
|
||||
M,J = np.meshgrid(range(0,m+1),range(0,m+1))
|
||||
|
||||
coeffs = bb(J,M) / sp.misc.factorial(J) * sp.exp( -lengthScale**(-2) ) *\
|
||||
(lengthScale**(-2))**J *magnSigma2
|
||||
|
||||
coeffs_dl = np.sum( coeffs*lengthScale**(-3)*(2.0-2.0*J*lengthScale**2),0)
|
||||
|
||||
coeffs = np.sum(coeffs,0)
|
||||
|
||||
else:
|
||||
coeffs = 2*magnSigma2*sp.exp( -lengthScale**(-2) ) * special.iv(range(0,m+1),1.0/lengthScale**(2))
|
||||
if np.any( np.isfinite(coeffs) == False):
|
||||
raise ValueError("sde_standard_periodic: Coefficients are not finite!")
|
||||
#import pdb; pdb.set_trace()
|
||||
coeffs[0] = 0.5*coeffs[0]
|
||||
|
||||
# Derivatives wrt (lengthScale)
|
||||
coeffs_dl = np.zeros(m+1)
|
||||
coeffs_dl[1:] = magnSigma2*lengthScale**(-3) * sp.exp(-lengthScale**(-2))*\
|
||||
(-4*special.iv(range(0,m),lengthScale**(-2)) + 4*(1+np.arange(1,m+1)*lengthScale**(2))*special.iv(range(1,m+1),lengthScale**(-2)) )
|
||||
|
||||
# The first element
|
||||
coeffs_dl[0] = magnSigma2*lengthScale**(-3) * np.exp(-lengthScale**(-2))*\
|
||||
(2*special.iv(0,lengthScale**(-2)) - 2*special.iv(1,lengthScale**(-2)) )
|
||||
|
||||
|
||||
return coeffs, coeffs_dl
|
||||
103
GPy/kern/src/sde_static.py
Normal file
103
GPy/kern/src/sde_static.py
Normal file
|
|
@ -0,0 +1,103 @@
|
|||
# -*- coding: utf-8 -*-
|
||||
# Copyright (c) 2015, Alex Grigorevskiy, Arno Solin
|
||||
# Licensed under the BSD 3-clause license (see LICENSE.txt)
|
||||
"""
|
||||
Classes in this module enhance Static covariance functions with the
|
||||
Stochastic Differential Equation (SDE) functionality.
|
||||
"""
|
||||
from .static import White
|
||||
from .static import Bias
|
||||
|
||||
import numpy as np
|
||||
|
||||
class sde_White(White):
|
||||
"""
|
||||
|
||||
Class provide extra functionality to transfer this covariance function into
|
||||
SDE forrm.
|
||||
|
||||
White kernel:
|
||||
|
||||
.. math::
|
||||
|
||||
k(x,y) = \alpha*\delta(x-y)
|
||||
|
||||
"""
|
||||
|
||||
def sde_update_gradient_full(self, gradients):
|
||||
"""
|
||||
Update gradient in the order in which parameters are represented in the
|
||||
kernel
|
||||
"""
|
||||
|
||||
self.variance.gradient = gradients[0]
|
||||
|
||||
def sde(self):
|
||||
"""
|
||||
Return the state space representation of the covariance.
|
||||
"""
|
||||
|
||||
variance = float(self.variance.values)
|
||||
|
||||
F = np.array( ((-np.inf,),) )
|
||||
L = np.array( ((1.0,),) )
|
||||
Qc = np.array( ((variance,),) )
|
||||
H = np.array( ((1.0,),) )
|
||||
|
||||
Pinf = np.array( ((variance,),) )
|
||||
P0 = Pinf.copy()
|
||||
|
||||
dF = np.zeros((1,1,1))
|
||||
dQc = np.zeros((1,1,1))
|
||||
dQc[:,:,0] = np.array( ((1.0,),) )
|
||||
|
||||
dPinf = np.zeros((1,1,1))
|
||||
dPinf[:,:,0] = np.array( ((1.0,),) )
|
||||
dP0 = dPinf.copy()
|
||||
|
||||
return (F, L, Qc, H, Pinf, P0, dF, dQc, dPinf, dP0)
|
||||
|
||||
|
||||
class sde_Bias(Bias):
|
||||
"""
|
||||
|
||||
Class provide extra functionality to transfer this covariance function into
|
||||
SDE forrm.
|
||||
|
||||
Bias kernel:
|
||||
|
||||
.. math::
|
||||
|
||||
k(x,y) = \alpha
|
||||
|
||||
"""
|
||||
def sde_update_gradient_full(self, gradients):
|
||||
"""
|
||||
Update gradient in the order in which parameters are represented in the
|
||||
kernel
|
||||
"""
|
||||
|
||||
self.variance.gradient = gradients[0]
|
||||
|
||||
def sde(self):
|
||||
"""
|
||||
Return the state space representation of the covariance.
|
||||
"""
|
||||
variance = float(self.variance.values)
|
||||
|
||||
F = np.array( ((0.0,),))
|
||||
L = np.array( ((1.0,),))
|
||||
Qc = np.zeros((1,1))
|
||||
H = np.array( ((1.0,),))
|
||||
|
||||
Pinf = np.zeros((1,1))
|
||||
P0 = np.array( ((variance,),) )
|
||||
|
||||
dF = np.zeros((1,1,1))
|
||||
dQc = np.zeros((1,1,1))
|
||||
|
||||
dPinf = np.zeros((1,1,1))
|
||||
dP0 = np.zeros((1,1,1))
|
||||
dP0[:,:,0] = np.array( ((1.0,),) )
|
||||
|
||||
return (F, L, Qc, H, Pinf, P0, dF, dQc, dPinf, dP0)
|
||||
192
GPy/kern/src/sde_stationary.py
Normal file
192
GPy/kern/src/sde_stationary.py
Normal file
|
|
@ -0,0 +1,192 @@
|
|||
# -*- coding: utf-8 -*-
|
||||
# Copyright (c) 2015, Alex Grigorevskiy, Arno Solin
|
||||
# Licensed under the BSD 3-clause license (see LICENSE.txt)
|
||||
"""
|
||||
Classes in this module enhance several stationary covariance functions with the
|
||||
Stochastic Differential Equation (SDE) functionality.
|
||||
"""
|
||||
from .rbf import RBF
|
||||
from .stationary import Exponential
|
||||
from .stationary import RatQuad
|
||||
|
||||
import numpy as np
|
||||
import scipy as sp
|
||||
|
||||
class sde_RBF(RBF):
|
||||
"""
|
||||
|
||||
Class provide extra functionality to transfer this covariance function into
|
||||
SDE form.
|
||||
|
||||
Radial Basis Function kernel:
|
||||
|
||||
.. math::
|
||||
|
||||
k(r) = \sigma^2 \exp \\bigg(- \\frac{1}{2} r^2 \\bigg) \\ \\ \\ \\ \text{ where } r = \sqrt{\sum_{i=1}^{input dim} \frac{(x_i-y_i)^2}{\ell_i^2} }
|
||||
|
||||
"""
|
||||
def sde_update_gradient_full(self, gradients):
|
||||
"""
|
||||
Update gradient in the order in which parameters are represented in the
|
||||
kernel
|
||||
"""
|
||||
|
||||
self.variance.gradient = gradients[0]
|
||||
self.lengthscale.gradient = gradients[1]
|
||||
|
||||
def sde(self):
|
||||
"""
|
||||
Return the state space representation of the covariance.
|
||||
"""
|
||||
|
||||
N = 10# approximation order ( number of terms in exponent series expansion)
|
||||
roots_rounding_decimals = 6
|
||||
|
||||
fn = np.math.factorial(N)
|
||||
|
||||
kappa = 1.0/2.0/self.lengthscale**2
|
||||
|
||||
Qc = np.array((self.variance*np.sqrt(np.pi/kappa)*fn*(4*kappa)**N,),)
|
||||
|
||||
pp = np.zeros((2*N+1,)) # array of polynomial coefficients from higher power to lower
|
||||
|
||||
for n in range(0, N+1): # (2N+1) - number of polynomial coefficients
|
||||
pp[2*(N-n)] = fn*(4.0*kappa)**(N-n)/np.math.factorial(n)*(-1)**n
|
||||
|
||||
pp = sp.poly1d(pp)
|
||||
roots = sp.roots(pp)
|
||||
|
||||
neg_real_part_roots = roots[np.round(np.real(roots) ,roots_rounding_decimals) < 0]
|
||||
aa = sp.poly1d(neg_real_part_roots, r=True).coeffs
|
||||
|
||||
F = np.diag(np.ones((N-1,)),1)
|
||||
F[-1,:] = -aa[-1:0:-1]
|
||||
|
||||
L= np.zeros((N,1))
|
||||
L[N-1,0] = 1
|
||||
|
||||
H = np.zeros((1,N))
|
||||
H[0,0] = 1
|
||||
|
||||
# Infinite covariance:
|
||||
Pinf = sp.linalg.solve_lyapunov(F, -np.dot(L,np.dot( Qc[0,0],L.T)))
|
||||
Pinf = 0.5*(Pinf + Pinf.T)
|
||||
# Allocating space for derivatives
|
||||
dF = np.empty([F.shape[0],F.shape[1],2])
|
||||
dQc = np.empty([Qc.shape[0],Qc.shape[1],2])
|
||||
dPinf = np.empty([Pinf.shape[0],Pinf.shape[1],2])
|
||||
|
||||
# Derivatives:
|
||||
dFvariance = np.zeros(F.shape)
|
||||
dFlengthscale = np.zeros(F.shape)
|
||||
dFlengthscale[-1,:] = -aa[-1:0:-1]/self.lengthscale * np.arange(-N,0,1)
|
||||
|
||||
dQcvariance = Qc/self.variance
|
||||
dQclengthscale = np.array(((self.variance*np.sqrt(2*np.pi)*fn*2**N*self.lengthscale**(-2*N)*(1-2*N,),)))
|
||||
|
||||
dPinf_variance = Pinf/self.variance
|
||||
|
||||
lp = Pinf.shape[0]
|
||||
coeff = np.arange(1,lp+1).reshape(lp,1) + np.arange(1,lp+1).reshape(1,lp) - 2
|
||||
coeff[np.mod(coeff,2) != 0] = 0
|
||||
dPinf_lengthscale = -1/self.lengthscale*Pinf*coeff
|
||||
|
||||
dF[:,:,0] = dFvariance
|
||||
dF[:,:,1] = dFlengthscale
|
||||
dQc[:,:,0] = dQcvariance
|
||||
dQc[:,:,1] = dQclengthscale
|
||||
dPinf[:,:,0] = dPinf_variance
|
||||
dPinf[:,:,1] = dPinf_lengthscale
|
||||
|
||||
P0 = Pinf.copy()
|
||||
dP0 = dPinf.copy()
|
||||
|
||||
# Benefits of this are not very sound. Helps only in one case:
|
||||
# SVD Kalman + RBF kernel
|
||||
import GPy.models.state_space_main as ssm
|
||||
(F, L, Qc, H, Pinf, P0, dF, dQc, dPinf,dP0, T) = ssm.balance_ss_model(F, L, Qc, H, Pinf, P0, dF, dQc, dPinf, dP0 )
|
||||
|
||||
return (F, L, Qc, H, Pinf, P0, dF, dQc, dPinf, dP0)
|
||||
|
||||
class sde_Exponential(Exponential):
|
||||
"""
|
||||
|
||||
Class provide extra functionality to transfer this covariance function into
|
||||
SDE form.
|
||||
|
||||
Exponential kernel:
|
||||
|
||||
.. math::
|
||||
|
||||
k(r) = \sigma^2 \exp \\bigg(- \\frac{1}{2} r \\bigg) \\ \\ \\ \\ \text{ where } r = \sqrt{\sum_{i=1}^{input dim} \frac{(x_i-y_i)^2}{\ell_i^2} }
|
||||
|
||||
"""
|
||||
|
||||
def sde_update_gradient_full(self, gradients):
|
||||
"""
|
||||
Update gradient in the order in which parameters are represented in the
|
||||
kernel
|
||||
"""
|
||||
|
||||
self.variance.gradient = gradients[0]
|
||||
self.lengthscale.gradient = gradients[1]
|
||||
|
||||
def sde(self):
|
||||
"""
|
||||
Return the state space representation of the covariance.
|
||||
"""
|
||||
variance = float(self.variance.values)
|
||||
lengthscale = float(self.lengthscale)
|
||||
|
||||
F = np.array(((-1.0/lengthscale,),))
|
||||
L = np.array(((1.0,),))
|
||||
Qc = np.array( ((2.0*variance/lengthscale,),) )
|
||||
H = np.array(((1.0,),))
|
||||
Pinf = np.array(((variance,),))
|
||||
P0 = Pinf.copy()
|
||||
|
||||
dF = np.zeros((1,1,2));
|
||||
dQc = np.zeros((1,1,2));
|
||||
dPinf = np.zeros((1,1,2));
|
||||
|
||||
dF[:,:,0] = 0.0
|
||||
dF[:,:,1] = 1.0/lengthscale**2
|
||||
|
||||
dQc[:,:,0] = 2.0/lengthscale
|
||||
dQc[:,:,1] = -2.0*variance/lengthscale**2
|
||||
|
||||
dPinf[:,:,0] = 1.0
|
||||
dPinf[:,:,1] = 0.0
|
||||
|
||||
dP0 = dPinf.copy()
|
||||
|
||||
return (F, L, Qc, H, Pinf, P0, dF, dQc, dPinf, dP0)
|
||||
|
||||
class sde_RatQuad(RatQuad):
|
||||
"""
|
||||
|
||||
Class provide extra functionality to transfer this covariance function into
|
||||
SDE form.
|
||||
|
||||
Rational Quadratic kernel:
|
||||
|
||||
.. math::
|
||||
|
||||
k(r) = \sigma^2 \\bigg( 1 + \\frac{r^2}{2} \\bigg)^{- \alpha} \\ \\ \\ \\ \text{ where } r = \sqrt{\sum_{i=1}^{input dim} \frac{(x_i-y_i)^2}{\ell_i^2} }
|
||||
|
||||
"""
|
||||
|
||||
def sde(self):
|
||||
"""
|
||||
Return the state space representation of the covariance.
|
||||
"""
|
||||
|
||||
assert False, 'Not Implemented'
|
||||
|
||||
# Params to use:
|
||||
|
||||
# self.lengthscale
|
||||
# self.variance
|
||||
#self.power
|
||||
|
||||
#return (F, L, Qc, H, Pinf, dF, dQc, dPinf)
|
||||
|
|
@ -320,6 +320,18 @@ class Exponential(Stationary):
|
|||
def dK_dr(self, r):
|
||||
return -0.5*self.K_of_r(r)
|
||||
|
||||
# def sde(self):
|
||||
# """
|
||||
# Return the state space representation of the covariance.
|
||||
# """
|
||||
# F = np.array([[-1/self.lengthscale]])
|
||||
# L = np.array([[1]])
|
||||
# Qc = np.array([[2*self.variance/self.lengthscale]])
|
||||
# H = np.array([[1]])
|
||||
# Pinf = np.array([[self.variance]])
|
||||
# # TODO: return the derivatives as well
|
||||
#
|
||||
# return (F, L, Qc, H, Pinf)
|
||||
|
||||
|
||||
|
||||
|
|
@ -388,6 +400,41 @@ class Matern32(Stationary):
|
|||
F1lower = np.array([f(lower) for f in F1])[:, None]
|
||||
return(self.lengthscale ** 3 / (12.*np.sqrt(3) * self.variance) * G + 1. / self.variance * np.dot(Flower, Flower.T) + self.lengthscale ** 2 / (3.*self.variance) * np.dot(F1lower, F1lower.T))
|
||||
|
||||
def sde(self):
|
||||
"""
|
||||
Return the state space representation of the covariance.
|
||||
"""
|
||||
variance = float(self.variance.values)
|
||||
lengthscale = float(self.lengthscale.values)
|
||||
foo = np.sqrt(3.)/lengthscale
|
||||
F = np.array([[0, 1], [-foo**2, -2*foo]])
|
||||
L = np.array([[0], [1]])
|
||||
Qc = np.array([[12.*np.sqrt(3) / lengthscale**3 * variance]])
|
||||
H = np.array([[1, 0]])
|
||||
Pinf = np.array([[variance, 0],
|
||||
[0, 3.*variance/(lengthscale**2)]])
|
||||
# Allocate space for the derivatives
|
||||
dF = np.empty([F.shape[0],F.shape[1],2])
|
||||
dQc = np.empty([Qc.shape[0],Qc.shape[1],2])
|
||||
dPinf = np.empty([Pinf.shape[0],Pinf.shape[1],2])
|
||||
# The partial derivatives
|
||||
dFvariance = np.zeros([2,2])
|
||||
dFlengthscale = np.array([[0,0],
|
||||
[6./lengthscale**3,2*np.sqrt(3)/lengthscale**2]])
|
||||
dQcvariance = np.array([12.*np.sqrt(3)/lengthscale**3])
|
||||
dQclengthscale = np.array([-3*12*np.sqrt(3)/lengthscale**4*variance])
|
||||
dPinfvariance = np.array([[1,0],[0,3./lengthscale**2]])
|
||||
dPinflengthscale = np.array([[0,0],
|
||||
[0,-6*variance/lengthscale**3]])
|
||||
# Combine the derivatives
|
||||
dF[:,:,0] = dFvariance
|
||||
dF[:,:,1] = dFlengthscale
|
||||
dQc[:,:,0] = dQcvariance
|
||||
dQc[:,:,1] = dQclengthscale
|
||||
dPinf[:,:,0] = dPinfvariance
|
||||
dPinf[:,:,1] = dPinflengthscale
|
||||
|
||||
return (F, L, Qc, H, Pinf, dF, dQc, dPinf)
|
||||
|
||||
class Matern52(Stationary):
|
||||
"""
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue