%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import pymc as pm
Sample data
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)
plt.plot(x, y, 'r.')
[<matplotlib.lines.Line2D at 0x1138e28d0>]
Kernel function
squared_exponential = lambda x, y, η2, ρ2: η2 * np.exp(-ρ2 * (x - y)**2)
Parameters of kernel function
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
x_pred = np.append(x, np.linspace(-6, 6, 10))
Covariance
@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
mu = np.zeros(len(x_pred))
Use masked array to make predictions on linear mesh along with observed values
y_masked = np.ma.masked_equal(np.append(y, np.zeros(10)), value=0)
y = pm.MvNormalCov('y', mu, Sigma, value=y_masked, observed=True)
M = pm.MCMC(locals())
M.use_step_method(pm.AdaptiveMetropolis, y)
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)
pm.Matplot.summary_plot(y)
Could not calculate Gelman-Rubin statistics. Requires multiple chains of equal length.
Plot predictions (blue points)
plt.plot(x, _y, 'r.')
plt.plot(np.linspace(-6, 6, 10), y.value[-10:], 'bo')
[<matplotlib.lines.Line2D at 0x113ab2470>]
pm.Matplot.plot(rho2)
Plotting ρ2
pm.Matplot.plot(eta2)
Plotting η2
pm.Matplot.plot(sigma2)
Plotting σ2