I am trying to draw samples from a GP prior with Gibbs covariance - following the example in the API on the website if I just were to replace the tanh_func with a function of my own choosing (in this case I am sampling the warping function from another GP) – it is very clear that the resulting covariance K is not positive semi-definite.
(I think it might be because if you zoom into the plot for the warping function there are these feeble discontinuities – I am not sure why they arise as a monotonic transformation of a perfectly smooth function shouldn’t result in any discontinuities. )
import numpy as np
import pymc3 as pm
import matplotlib.pyplot as plt
import theano.tensor as tt
X = np.linspace(-10, 10, 1000)[:,None]
y = (X**3 - X**2 + X)
lengthscale = 2
eta = 2.0
cov = eta**2 * pm.gp.cov.ExpQuad(1, lengthscale)
K_base = cov(X).eval()
plt.matshow(K_base)
###### Samples from Gibbs covariance
def warp_gp_func(x, K):
sample = pm.MvNormal.dist(mu=np.zeros(K.shape[0]), cov=K + np.eye(1000)*1e-5, shape=(len(X),)).random(size=1).T*tt.ones(shape=(1000,1))
return sample.exp()
def tanh_func(x, ls1, ls2, w, x0):
"""
ls1: left saturation value
ls2: right saturation value
w: transition width
x0: transition location.
"""
return (ls1 + ls2) / 2.0 - (ls1 - ls2) / 2.0 * tt.tanh((x - x0) / w)
ls1 = 0.05
ls2 = 0.6
w = 0.3
x0 = 1.0
#cov = pm.gp.cov.Gibbs(1, tanh_func, args=(ls1, ls2, w, x0))
cov = pm.gp.cov.Gibbs(1, warp_gp_func, args=(K_base))
#th = tanh_func(X, ls1, ls2, w, x0).eval()
wf = warp_gp_func(X, K_base).eval()
plt.figure(figsize=(14, 4))
plt.plot(X, wf)
plt.ylabel("lengthscale")
plt.xlabel("X")
plt.title("Lengthscale as a function of X")
K = cov(X).eval()
plt.matshow(K)
### Below block throws an error
plt.figure(figsize=(14, 4))
plt.plot(X, pm.MvNormal.dist(mu=np.zeros(K.shape[0]), cov=K + np.eye(1000)*1e-6, shape=(1000,)).random(size=3).T)
plt.title("Samples from the GP prior")
plt.ylabel("y")
plt.xlabel("X");
