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

Sample data

In [21]:
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 0x1138e28d0>]

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 [10]:
y_masked = np.ma.masked_equal(np.append(y, np.zeros(10)), value=0)
In [11]:
y = pm.MvNormalCov('y', mu, Sigma, value=y_masked, observed=True)
In [12]:
M = pm.MCMC(locals())
M.use_step_method(pm.AdaptiveMetropolis, y)
In [13]:
M.sample(5000, 4000)
 [-----------------100%-----------------] 5000 of 5000 complete in 615.7 sec
/Users/fonnescj/GitHub/pymc/pymc/StepMethods.py:1272: UserWarning: 
Covariance was not positive definite and proposal_sd cannot be computed by 
Cholesky decomposition. The next jumps will be based on the last 
valid covariance matrix. This situation may have arisen because no 
jumps were accepted during the last `interval`. One solution is to 
increase the interval, or specify an initial covariance matrix with 
a smaller variance. For this simulation, each time a similar error 
occurs, proposal_sd will be reduced by a factor .9 to reduce the 
jumps and increase the likelihood of accepted jumps.
  warnings.warn(adjustmentwarning)
In [14]:
pm.Matplot.summary_plot(y)
Could not calculate Gelman-Rubin statistics. Requires multiple chains of equal length.

Plot predictions (blue points)

In [22]:
plt.plot(x, _y, 'r.')
plt.plot(np.linspace(-6, 6, 10), y.value[-10:], 'bo')
Out[22]:
[<matplotlib.lines.Line2D at 0x113ab2470>]
In [23]:
pm.Matplot.plot(rho2)
Plotting ρ2
In [24]:
pm.Matplot.plot(eta2)
Plotting η2
In [25]:
pm.Matplot.plot(sigma2)
Plotting σ2