GPMCMC

Put the following in PyMCModel.py:

import pymc as pm
import pymc.gp as gp
from pymc.gp.cov_funs import matern
import numpy as np
import matplotlib.pyplot as pl

from numpy.random import normal

x = np.arange(-1.,1.,.1)

# Prior parameters of C
diff_degree = pm.Uniform('diff_degree', 1., 3)
amp = pm.Lognormal('amp', mu=.4, tau=1.)
scale = pm.Lognormal('scale', mu=.5, tau=1.)

# The covariance dtrm C is valued as a Covariance object.
@pm.deterministic
def C(eval_fun = gp.matern.euclidean, diff_degree=diff_degree, amp=amp, scale=scale):
    return gp.NearlyFullRankCovariance(eval_fun, diff_degree=diff_degree, amp=amp, scale=scale)


# Prior parameters of M
a = pm.Normal('a', mu=1., tau=1.)
b = pm.Normal('b', mu=.5, tau=1.)
c = pm.Normal('c', mu=2., tau=1.)

# The mean M is valued as a Mean object.
def linfun(x, a, b, c):
    # return a * x ** 2 + b * x + c
    return 0.*x + c
@pm.deterministic
def M(eval_fun = linfun, a=a, b=b, c=c):
    return gp.Mean(eval_fun, a=a, b=b, c=c)

# The GP submodel
fmesh = np.linspace(-np.pi/3.3,np.pi/3.3,4)
sm = gp.GPSubmodel('sm',M,C,fmesh)

# Observation precision
V = .0001

# The data d is just array-valued. It's normally distributed about GP.f(obs_x).
init_val = np.random.normal(size=len(fmesh))
d = pm.Normal('d',mu=sm.f_eval, tau=1./V, value=init_val, observed=True)

Then execute the following:

import PyMCmodel
from pymc import *
from pylab import *
from PyMCmodel import *

x = linspace(-1,1,400)

GPSampler = MCMC(PyMCmodel)
# Uncomment this to use the GPNormal step method instead of the default GPMetropolis
GPSampler.use_step_method(gp.GPEvaluationGibbs, GPSampler.sm, GPSampler.V, GPSampler.d)

GPSampler.isample(iter=5000,burn=1000,thin=100)

# Uncomment this for a medium run.
# GPSampler.isample(iter=500,burn=0,thin=10)

# Uncomment this for a short run.
# GPSampler.isample(iter=50,burn=0,thin=1)


N_samps = len(GPSampler.sm.f.trace())

close('all')

mid_traces = []
subplot(1,2,1)
for i in range(0,N_samps):
    f=GPSampler.sm.f.trace()[i](x)
    plot(x,f)
    mid_traces.append(f[len(f)/2])
    plot(fmesh,GPSampler.d.value,'k.',markersize=16)
axis([x.min(),x.max(),-5.,10.])
title('Some samples of f')

subplot(1,2,2)

plot(mid_traces)
title('The trace of f at midpoint')

# Plot posterior of C and tau
figure()
subplot(2,2,1)
try:
    plot(GPSampler.diff_degree.trace())
    title("Degree of differentiability of f")
except:
    pass

subplot(2,2,2)
plot(GPSampler.amp.trace())
title("Pointwise prior variance of f")

subplot(2,2,3)
plot(GPSampler.scale.trace())
title("X-axis scaling of f")

# subplot(2,2,4)
# plot(GPSampler.V.trace())
# title('Observation precision')


# Plot posterior of M
figure()
subplot(1,3,1)
plot(GPSampler.a.trace())
title("Quadratic coefficient of M")

subplot(1,3,2)
plot(GPSampler.b.trace())
title("Linear coefficient of M")

subplot(1,3,3)
plot(GPSampler.c.trace())
title("Constant term of M")