gp-bug-example.py

#!/usr/bin/env python

import pymc
from pymc import gp
from pymc.gp.cov_funs import matern,gaussian
from pylab import *

# Load some data generated from a GP with mean=0, scale=1, amp=1
xdata,ydata = loadtxt('train.txt', unpack=1)

# Try to define non-informative prior for GP mean (beta).
# Uninformative gives results I would expect, but Normal with large
# variance gives unusually large posterior variance
# beta = pymc.Uninformative('beta', value=0)
beta = pymc.Normal('beta', 0, 0.000001, value=0)

# Other priors
scale = pymc.Exponential('scale', 1e-3, value=2)
amp = pymc.Exponential('amp', 1e-9, value=2)

def const_mean(x, beta):
    return x*0.0 + beta

@pymc.deterministic
def M(eval_fun = const_mean, beta=beta):
    return gp.Mean(eval_fun, beta=beta)

@pymc.deterministic
def C(eval_fun = gaussian.euclidean, amp=amp, scale=scale):
    return gp.Covariance(eval_fun, amp=amp, scale=scale)

sm = gp.GPSubmodel('sm', M,C, xdata, init_vals=ydata, obs_on_mesh=True)

last_gen = set([sm.f])

print pymc.utils.crawl_dataless(set(last_gen),[last_gen])



GPSampler = pymc.MCMC()
# GPSampler.assign_step_methods()
# print GPSampler.step_method_dict[beta]
# 
# GPSampler.isample(iter=5000, burn=1000, thin=10)
# 
# hist( GPSampler.beta.trace(), 25 )
# show()