In [1]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import pymc as pm

Sample data

In [2]:
x = (-5, -4.9, -4.8, -4.7, -4.6, -4.5, -4.4, -4.3, -4.2, -4.1, -4, 
-3.9, -3.8, -3.7, -3.6, -3.5, -3.4, -3.3, -3.2, -3.1, -3, -2.9, 
-2.8, -2.7, -2.6, -2.5, -2.4, -2.3, -2.2, -2.1, -2, -1.9, -1.8, 
-1.7, -1.6, -1.5, -1.4, -1.3, -1.2, -1.1, -1, -0.9, -0.8, -0.7, 
-0.6, -0.5, -0.4, -0.3, -0.2, -0.1, 0, 0.1, 0.2, 0.3, 0.4, 0.5, 
0.6, 0.7, 0.8, 0.9, 1, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 
1.9, 2, 2.1, 2.2, 2.3, 2.4, 2.5, 2.6, 2.7, 2.8, 2.9, 3, 3.1, 
3.2, 3.3, 3.4, 3.5, 3.6, 3.7, 3.8, 3.9, 4, 4.1, 4.2, 4.3, 4.4, 
4.5, 4.6, 4.7, 4.8, 4.9, 5)

y = (1.04442478194401, 0.948306088493654, 0.357037759697332, 0.492336514646604, 
0.520651364364746, 0.112629866592809, 0.470995468454158, -0.168442254267804, 
0.0720344402575861, -0.188108980535916, -0.0160163306512027, 
-0.0388792158617705, -0.0600673630622568, 0.113568725264636, 
0.447160403837629, 0.664421188556779, -0.139510743820276, 0.458823971660986, 
0.141214654640904, -0.286957663528091, -0.466537724021695, -0.308185884317105, 
-1.57664872694079, -1.44463024170082, -1.51206214603847, -1.49393593601901, 
-2.02292464164487, -1.57047488853653, -1.22973445533419, -1.51502367058357, 
-1.41493587255224, -1.10140254663611, -0.591866485375275, -1.08781838696462, 
-0.800375653733931, -1.00764767602679, -0.0471028950122742, -0.536820626879737, 
-0.151688056391446, -0.176771681318393, -0.240094952335518, -1.16827876746502, 
-0.493597351974992, -0.831683011472805, -0.152347043914137, 0.0190364158178343, 
-1.09355955218051, -0.328157917911376, -0.585575679802941, -0.472837120425201, 
-0.503633622750049, -0.0124446353828312, -0.465529814250314, 
-0.101621725887347, -0.26988462590405, 0.398726664193302, 0.113805181040188, 
0.331353802465398, 0.383592361618461, 0.431647298655434, 0.580036473774238, 
0.830404669466897, 1.17919105883462, 0.871037583886711, 1.12290553424174, 
0.752564860804382, 0.76897960270623, 1.14738839410786, 0.773151715269892, 
0.700611498974798, 0.0412951045437818, 0.303526087747629, -0.139399513324585, 
-0.862987735433697, -1.23399179134008, -1.58924289116396, -1.35105117911049, 
-0.990144529089174, -1.91175364127672, -1.31836236129543, -1.65955735224704, 
-1.83516148300526, -2.03817062501248, -1.66764011409214, -0.552154350554687, 
-0.547807883952654, -0.905389222477036, -0.737156477425302, -0.40211249920415, 
0.129669958952991, 0.271142753510592, 0.176311762529962, 0.283580281859344, 
0.635808289696458, 1.69976647982837, 1.10748978734239, 0.365412229181044, 
0.788821368082444, 0.879731888124867, 1.02180766619069, 0.551526067300283
)

N = len(x)
In [3]:
plt.plot(x, y, 'r.')
Out[3]:
[<matplotlib.lines.Line2D at 0x118442588>]

Kernel function

In [4]:
squared_exponential = lambda x, y, η2, ρ2: η2 * np.exp(-ρ2 * (x - y)**2)

Parameters of kernel function

In [5]:
eta2 = pm.HalfCauchy('eta2', 0, 5, value=1)
rho2 = pm.HalfCauchy('rho2', 0, 5, value=1)
sigma2 = pm.HalfCauchy('sigma2', 0, 5, value=1)

Points on which to evaluate GP

In [6]:
x_pred = np.append(x, np.linspace(-6, 6, 10))

Covariance

In [7]:
@pm.deterministic
def Sigma(eta2=eta2, rho2=rho2, sigma2=sigma2):
    S = np.array([[squared_exponential(xi, yi, eta2, rho2) for xi in x_pred] for yi in x_pred])
    # Add variance to 
    S[np.diag_indices_from(S)] += sigma2
    return S

Default mean

In [8]:
mu = np.zeros(len(x_pred))

Use masked array to make predictions on linear mesh along with observed values

In [9]:
y_masked = np.ma.masked_equal(np.append(y, np.zeros(10)), value=0)
In [10]:
y_masked
Out[10]:
masked_array(data = [1.04442478194401 0.948306088493654 0.357037759697332 0.492336514646604
 0.520651364364746 0.112629866592809 0.470995468454158 -0.168442254267804
 0.0720344402575861 -0.188108980535916 -0.0160163306512027
 -0.0388792158617705 -0.0600673630622568 0.113568725264636
 0.447160403837629 0.664421188556779 -0.139510743820276 0.458823971660986
 0.141214654640904 -0.286957663528091 -0.466537724021695 -0.308185884317105
 -1.57664872694079 -1.44463024170082 -1.51206214603847 -1.49393593601901
 -2.02292464164487 -1.57047488853653 -1.22973445533419 -1.51502367058357
 -1.41493587255224 -1.10140254663611 -0.591866485375275 -1.08781838696462
 -0.800375653733931 -1.00764767602679 -0.0471028950122742
 -0.536820626879737 -0.151688056391446 -0.176771681318393
 -0.240094952335518 -1.16827876746502 -0.493597351974992 -0.831683011472805
 -0.152347043914137 0.0190364158178343 -1.09355955218051 -0.328157917911376
 -0.585575679802941 -0.472837120425201 -0.503633622750049
 -0.0124446353828312 -0.465529814250314 -0.101621725887347
 -0.26988462590405 0.398726664193302 0.113805181040188 0.331353802465398
 0.383592361618461 0.431647298655434 0.580036473774238 0.830404669466897
 1.17919105883462 0.871037583886711 1.12290553424174 0.752564860804382
 0.76897960270623 1.14738839410786 0.773151715269892 0.700611498974798
 0.0412951045437818 0.303526087747629 -0.139399513324585 -0.862987735433697
 -1.23399179134008 -1.58924289116396 -1.35105117911049 -0.990144529089174
 -1.91175364127672 -1.31836236129543 -1.65955735224704 -1.83516148300526
 -2.03817062501248 -1.66764011409214 -0.552154350554687 -0.547807883952654
 -0.905389222477036 -0.737156477425302 -0.40211249920415 0.129669958952991
 0.271142753510592 0.176311762529962 0.283580281859344 0.635808289696458
 1.69976647982837 1.10748978734239 0.365412229181044 0.788821368082444
 0.879731888124867 1.02180766619069 0.551526067300283 -- -- -- -- -- -- --
 -- -- --],
             mask = [False False False False False False False False False False False False
 False False False False False False False False False False False False
 False False False False False False False False False False False False
 False False False False False False False False False False False False
 False False False False False False False False False False False False
 False False False False False False False False False False False False
 False False False False False False False False False False False False
 False False False False False False False False False False False False
 False False False False False  True  True  True  True  True  True  True
  True  True  True],
       fill_value = 0.0)
In [11]:
gp = pm.MvNormalCov('gp', mu, Sigma, value=y_masked, observed=True)
In [12]:
M = pm.MCMC(locals())
#M.use_step_method(pm.AdaptiveMetropolis, y)
In [13]:
M.sample(2000, 1000)
 [-----------------100%-----------------] 2000 of 2000 complete in 255.6 sec
In [14]:
pm.Matplot.summary_plot(gp)
Could not calculate Gelman-Rubin statistics. Requires multiple chains of equal length.

Plot predictions (blue points)

In [17]:
plt.plot(x, y, 'r.')
plt.plot(np.linspace(-6, 6, 10), gp.value[-10:], 'bo')
Out[17]:
[<matplotlib.lines.Line2D at 0x11ee49898>]
In [18]:
pm.Matplot.plot(rho2)
Plotting rho2
In [19]:
pm.Matplot.plot(eta2)
Plotting eta2
In [20]:
pm.Matplot.plot(sigma2)
Plotting sigma2