From d3730bb5185fd70366f56200bdd2cc6b302fd795 Mon Sep 17 00:00:00 2001 From: Ilias Bilionis Date: Tue, 4 Aug 2015 01:18:32 -0400 Subject: [PATCH] Fixed MCMC sampler. --- GPy/core/parameterization/parameter_core.py | 10 ++++++ GPy/inference/mcmc/__init__.py | 1 + GPy/inference/mcmc/samplers.py | 32 ++++++++---------- ib_tests/test_regression.py | 36 +++++++++++++++++++++ 4 files changed, 61 insertions(+), 18 deletions(-) create mode 100644 ib_tests/test_regression.py diff --git a/GPy/core/parameterization/parameter_core.py b/GPy/core/parameterization/parameter_core.py index 6add95b0..a0d0981d 100644 --- a/GPy/core/parameterization/parameter_core.py +++ b/GPy/core/parameterization/parameter_core.py @@ -684,6 +684,16 @@ class OptimizationHandlable(Indexable): if self._has_fixes(): return g[self._fixes_] return g + def _log_det_jacobian(self): + """ + Return the logarithm of the Jacobian needed for a proper change of + variables. + """ + J = np.ones(self.param_array.shape) + [np.put(J, i, c.jacobianfactor(self.param_array[i])) + for c, i in self.constraints.iteritems() if c != __fixed__] + return np.log(J).sum() + @property def num_params(self): """ diff --git a/GPy/inference/mcmc/__init__.py b/GPy/inference/mcmc/__init__.py index 956448d4..b30b6ff0 100644 --- a/GPy/inference/mcmc/__init__.py +++ b/GPy/inference/mcmc/__init__.py @@ -1 +1,2 @@ from hmc import HMC +from samplers import * diff --git a/GPy/inference/mcmc/samplers.py b/GPy/inference/mcmc/samplers.py index 444d99d7..c2367359 100644 --- a/GPy/inference/mcmc/samplers.py +++ b/GPy/inference/mcmc/samplers.py @@ -4,23 +4,18 @@ import numpy as np from scipy import linalg, optimize -import Tango import sys -import re -import numdifftools as ndt -import pdb -import cPickle class Metropolis_Hastings: def __init__(self,model,cov=None): """Metropolis Hastings, with tunings according to Gelman et al. """ self.model = model - current = self.model._get_params_transformed() + current = self.model.optimizer_array self.D = current.size self.chains = [] if cov is None: - self.cov = model.Laplace_covariance() + self.cov = np.eye(self.D) else: self.cov = cov self.scale = 2.4/np.sqrt(self.D) @@ -31,20 +26,20 @@ class Metropolis_Hastings: if start is None: self.model.randomize() else: - self.model._set_params_transformed(start) + self.model.optimizer_array = start - - - def sample(self, Ntotal, Nburn, Nthin, tune=True, tune_throughout=False, tune_interval=400): - current = self.model._get_params_transformed() - fcurrent = self.model.log_likelihood() + self.model.log_prior() + def sample(self, Ntotal=10000, Nburn=1000, Nthin=10, tune=True, tune_throughout=False, tune_interval=400): + current = self.model.optimizer_array + fcurrent = self.model.log_likelihood() + self.model.log_prior() + \ + self.model._log_det_jacobian() accepted = np.zeros(Ntotal,dtype=np.bool) for it in range(Ntotal): print "sample %d of %d\r"%(it,Ntotal), sys.stdout.flush() prop = np.random.multivariate_normal(current, self.cov*self.scale*self.scale) - self.model._set_params_transformed(prop) - fprop = self.model.log_likelihood() + self.model.log_prior() + self.model.optimizer_array = prop + fprop = self.model.log_likelihood() + self.model.log_prior() + \ + self.model._log_det_jacobian() if fprop>fcurrent:#sample accepted, going 'uphill' accepted[it] = True @@ -72,10 +67,11 @@ class Metropolis_Hastings: def predict(self,function,args): """Make a prediction for the function, to which we will pass the additional arguments""" - param = self.model._get_params() + param = self.model.param_array fs = [] for p in self.chain: - self.model._set_params(p) + self.model.param_array = p fs.append(function(*args)) - self.model._set_params(param)# reset model to starting state + # reset model to starting state + self.model.param_array = param return fs diff --git a/ib_tests/test_regression.py b/ib_tests/test_regression.py new file mode 100644 index 00000000..cd9e2943 --- /dev/null +++ b/ib_tests/test_regression.py @@ -0,0 +1,36 @@ +""" +Test the regression we get with the new transformations. + +Author: + Ilias Bilionis + +Date: + 3/8/2015 + +""" + + +import sys +import os +# Make sure we load the GP that is here +sys.path.insert(0, os.path.abspath(os.path.join(os.path.dirname(__file__), '..'))) +import GPy +import matplotlib.pyplot as plt +import numpy as np +import triangle + + +if __name__ == '__main__': + m = GPy.examples.regression.olympic_marathon_men(optimize=True) + plt.show(block=True) + print m + mcmc = GPy.inference.mcmc.samplers.Metropolis_Hastings(m) + mcmc.sample(Ntotal=100000, Nburn=10000, Nthin=100, tune_interval=1000, tune_throughout=True) + samples = np.array(mcmc.chains[-1]) + fig = triangle.corner(samples) + m.plot() + fig = plt.figure() + for i in xrange(samples.shape[1]): + ax = fig.add_subplot(samples.shape[1], 1, i + 1) + ax.plot(samples[:, i], linewidth=1.5) + plt.show(block=True)