syrte
4/7/2020 - 8:23 AM

## Custom GPy noise model

Custom GPy noise model

``````import GPy
import numpy as np

class HeteroNoise(GPy.kern.Kern):
"""
heteroscedastic white noise as function of X
"""

def __init__(self, input_dim, variance_function=1., active_dims=None):
"""
interpolating noise function:
variance_function=lambda X:np.interp(X[:, 0], x, yerr)
"""
super(HeteroNoise, self).__init__(input_dim, active_dims, name='noise')

if callable(variance_function):
self.variance_function = variance_function
if np.isscalar(variance_function):
self.variance_function = lambda X: np.full_like(X[:, 0], variance_function)
else:
raise TypeError("'variance_function' should be a function!")

def K(self, X, X2=None):
if X2 is None:
return np.diag(self.Kdiag(X))
else:
return np.zeros((X.shape[0], X2.shape[0]))

def Kdiag(self, X):
return self.variance_function(X)
``````
``````import GPy
import numpy as np

class Noise(GPy.kern.Kern):
def __init__(self, input_dim, variance=1., active_dims=None):
super(Noise, self).__init__(input_dim, active_dims, name='noise')
assert input_dim == 1, "For this kernel we assume input_dim=1"

from paramz.transformations import Logexp
self.variance = GPy.Param('variance', variance, Logexp())

def K(self, X, X2=None):
if X2 is None:
return np.diag(self.Kdiag(X)) * self.variance
else:
return np.zeros((X.shape[0], X2.shape[0]))

def Kdiag(self, X):
return np.interp(X[:, 0], gmag, cvar) * self.variance

if X2 is None:
dKdvar = np.diag(self.Kdiag(X))
else:

return np.sqrt((r / epsilon)**2 + 1)

def rbf_norm(x1, x2):
"""
x1, x2: shape (nsample, nfeature)

>>> rbf_norm(rand(10, 2), rand(5, 2)).shape
# (10, 5)
"""
return np.sqrt(((x1[..., None, :] - x2)**2).sum(axis=-1))

def rbf_eval(xi, x, w, epsilon=None, kern=rbf_kern_multiquadric):
"""
xi: shape (nsample, nfeature)
x: shape (nnode, nfeature)
nodes
w: shape (nnode,)
weights
"""
if epsilon is None:
epsilon = (x.max() - x.min()) / len(x)
r = rbf_norm(xi, x)
A = kern(r, epsilon)
return (A * w).sum(axis=-1)

class RbfNoise(GPy.kern.Kern):
def __init__(self, input_dim, x, w, p, rbf_kern=rbf_kern_multiquadric, active_dims=None):
"""
x: rbf nodes
w: weights of nodes
p: fiducial
"""
super(RbfNoise, self).__init__(input_dim, active_dims, name='rbfnoise')

self.epsilon = (x.max() - x.min()) / len(x)
self.rbf_kern = rbf_kern
self.x = x
self.w = GPy.Param('w', w)
self.p = GPy.Param('p', p)

if not np.isscalar(p):
raise NotImplementedError

def K(self, X, X2=None):
if X2 is None:
return np.diag(self.Kdiag(X))
else:
return np.zeros((X.shape[0], X2.shape[0]))

def lnKdiag(self, X):
lnKdiag = rbf_eval(X, self.x, self.w, epsilon=self.epsilon, kern=self.rbf_kern) + self.p
return lnKdiag

def Kdiag(self, X):
lnKdiag = self.lnKdiag(X)
return np.exp(lnKdiag)