mirror of
https://github.com/SheffieldML/GPy.git
synced 2026-05-30 14:35:15 +02:00
moving mcmc
This commit is contained in:
James Hensman
parent
70afcd3ddd
commit
d1200a1be7
8 changed files with 12 additions and 59 deletions
157
GPy/inference/mcmc/hmc.py
Normal file
157
GPy/inference/mcmc/hmc.py
Normal file
|
|
@ -0,0 +1,157 @@
|
|||
"""HMC implementation"""
|
||||
|
||||
import numpy as np
|
||||
|
||||
|
||||
class HMC:
|
||||
def __init__(self,model,M=None,stepsize=1e-1):
|
||||
self.model = model
|
||||
self.stepsize = stepsize
|
||||
self.p = np.empty_like(model.optimizer_array.copy())
|
||||
if M is None:
|
||||
self.M = np.eye(self.p.size)
|
||||
else:
|
||||
self.M = M
|
||||
self.Minv = np.linalg.inv(self.M)
|
||||
|
||||
def sample(self, m_iters=1000, hmc_iters=20):
|
||||
params = np.empty((m_iters,self.p.size))
|
||||
for i in xrange(m_iters):
|
||||
self.p[:] = np.random.multivariate_normal(np.zeros(self.p.size),self.M)
|
||||
H_old = self._computeH()
|
||||
theta_old = self.model.optimizer_array.copy()
|
||||
params[i] = self.model.unfixed_param_array
|
||||
#Matropolis
|
||||
self._update(hmc_iters)
|
||||
H_new = self._computeH()
|
||||
|
||||
if H_old>H_new:
|
||||
k = 1.
|
||||
else:
|
||||
k = np.exp(H_old-H_new)
|
||||
if np.random.rand()<k:
|
||||
params[i] = self.model.unfixed_param_array
|
||||
else:
|
||||
self.model.optimizer_array = theta_old
|
||||
return params
|
||||
|
||||
def _update(self, hmc_iters):
|
||||
for i in xrange(hmc_iters):
|
||||
self.p[:] += -self.stepsize/2.*self.model._transform_gradients(self.model.objective_function_gradients())
|
||||
self.model.optimizer_array = self.model.optimizer_array + self.stepsize*np.dot(self.Minv, self.p)
|
||||
self.p[:] += -self.stepsize/2.*self.model._transform_gradients(self.model.objective_function_gradients())
|
||||
|
||||
def _computeH(self,):
|
||||
return self.model.objective_function()+self.p.size*np.log(2*np.pi)/2.+np.log(np.linalg.det(self.M))/2.+np.dot(self.p, np.dot(self.Minv,self.p[:,None]))/2.
|
||||
|
||||
class HMC_shortcut:
|
||||
def __init__(self,model,M=None,stepsize_range=[1e-6, 1e-1],groupsize=5, Hstd_th=[1e-5, 3.]):
|
||||
self.model = model
|
||||
self.stepsize_range = np.log(stepsize_range)
|
||||
self.p = np.empty_like(model.optimizer_array.copy())
|
||||
self.groupsize = groupsize
|
||||
self.Hstd_th = Hstd_th
|
||||
if M is None:
|
||||
self.M = np.eye(self.p.size)
|
||||
else:
|
||||
self.M = M
|
||||
self.Minv = np.linalg.inv(self.M)
|
||||
|
||||
def sample(self, m_iters=1000, hmc_iters=20):
|
||||
params = np.empty((m_iters,self.p.size))
|
||||
for i in xrange(m_iters):
|
||||
# sample a stepsize from the uniform distribution
|
||||
stepsize = np.exp(np.random.rand()*(self.stepsize_range[1]-self.stepsize_range[0])+self.stepsize_range[0])
|
||||
self.p[:] = np.random.multivariate_normal(np.zeros(self.p.size),self.M)
|
||||
H_old = self._computeH()
|
||||
params[i] = self.model.unfixed_param_array
|
||||
theta_old = self.model.optimizer_array.copy()
|
||||
#Matropolis
|
||||
self._update(hmc_iters, stepsize)
|
||||
H_new = self._computeH()
|
||||
|
||||
if H_old>H_new:
|
||||
k = 1.
|
||||
else:
|
||||
k = np.exp(H_old-H_new)
|
||||
if np.random.rand()<k:
|
||||
params[i] = self.model.unfixed_param_array
|
||||
else:
|
||||
self.model.optimizer_array = theta_old
|
||||
return params
|
||||
|
||||
def _update(self, hmc_iters, stepsize):
|
||||
theta_buf = np.empty((2*hmc_iters+1,self.model.optimizer_array.size))
|
||||
p_buf = np.empty((2*hmc_iters+1,self.p.size))
|
||||
H_buf = np.empty((2*hmc_iters+1,))
|
||||
# Set initial position
|
||||
theta_buf[hmc_iters] = self.model.optimizer_array
|
||||
p_buf[hmc_iters] = self.p
|
||||
H_buf[hmc_iters] = self._computeH()
|
||||
|
||||
reversal = []
|
||||
pos = 1
|
||||
i=0
|
||||
while i<hmc_iters:
|
||||
self.p[:] += -stepsize/2.*self.model._transform_gradients(self.model.objective_function_gradients())
|
||||
self.model.optimizer_array = self.model.optimizer_array + stepsize*np.dot(self.Minv, self.p)
|
||||
self.p[:] += -stepsize/2.*self.model._transform_gradients(self.model.objective_function_gradients())
|
||||
|
||||
theta_buf[hmc_iters+pos] = self.model.optimizer_array
|
||||
p_buf[hmc_iters+pos] = self.p
|
||||
H_buf[hmc_iters+pos] = self._computeH()
|
||||
i+=1
|
||||
|
||||
if i<self.groupsize:
|
||||
pos += 1
|
||||
continue
|
||||
else:
|
||||
if len(reversal)==0:
|
||||
Hlist = range(hmc_iters+pos,hmc_iters+pos-self.groupsize,-1)
|
||||
if self._testH(H_buf[Hlist]):
|
||||
pos += 1
|
||||
else:
|
||||
# Reverse the trajectory for the 1st time
|
||||
reversal.append(pos)
|
||||
if hmc_iters-i>pos:
|
||||
pos = -1
|
||||
i += pos
|
||||
self.model.optimizer_array = theta_buf[hmc_iters]
|
||||
self.p[:] = -p_buf[hmc_iters]
|
||||
else:
|
||||
pos_new = pos-hmc_iters+i
|
||||
self.model.optimizer_array = theta_buf[hmc_iters+pos_new]
|
||||
self.p[:] = -p_buf[hmc_iters+pos_new]
|
||||
break
|
||||
else:
|
||||
Hlist = range(hmc_iters+pos,hmc_iters+pos+self.groupsize)
|
||||
# print Hlist
|
||||
# print self._testH(H_buf[Hlist])
|
||||
|
||||
if self._testH(H_buf[Hlist]):
|
||||
pos += -1
|
||||
else:
|
||||
# Reverse the trajectory for the 2nd time
|
||||
r = (hmc_iters - i)%((reversal[0]-pos)*2)
|
||||
if r>(reversal[0]-pos):
|
||||
pos_new = 2*reversal[0] - r - pos
|
||||
else:
|
||||
pos_new = pos + r
|
||||
self.model.optimizer_array = theta_buf[hmc_iters+pos_new]
|
||||
self.p[:] = p_buf[hmc_iters+pos_new] # the sign of momentum might be wrong!
|
||||
# print reversal[0],pos,pos_new
|
||||
# print H_buf
|
||||
break
|
||||
|
||||
def _testH(self, Hlist):
|
||||
Hstd = np.std(Hlist)
|
||||
# print Hlist
|
||||
# print Hstd
|
||||
if Hstd<self.Hstd_th[0] or Hstd>self.Hstd_th[1]:
|
||||
return False
|
||||
else:
|
||||
return True
|
||||
|
||||
def _computeH(self,):
|
||||
return self.model.objective_function()+self.p.size*np.log(2*np.pi)/2.+np.log(np.linalg.det(self.M))/2.+np.dot(self.p, np.dot(self.Minv,self.p[:,None]))/2.
|
||||
|
||||
81
GPy/inference/mcmc/samplers.py
Normal file
81
GPy/inference/mcmc/samplers.py
Normal file
|
|
@ -0,0 +1,81 @@
|
|||
# Copyright (c) 2012, GPy authors (see AUTHORS.txt).
|
||||
# Licensed under the BSD 3-clause license (see LICENSE.txt)
|
||||
|
||||
|
||||
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()
|
||||
self.D = current.size
|
||||
self.chains = []
|
||||
if cov is None:
|
||||
self.cov = model.Laplace_covariance()
|
||||
else:
|
||||
self.cov = cov
|
||||
self.scale = 2.4/np.sqrt(self.D)
|
||||
self.new_chain(current)
|
||||
|
||||
def new_chain(self, start=None):
|
||||
self.chains.append([])
|
||||
if start is None:
|
||||
self.model.randomize()
|
||||
else:
|
||||
self.model._set_params_transformed(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()
|
||||
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()
|
||||
|
||||
if fprop>fcurrent:#sample accepted, going 'uphill'
|
||||
accepted[it] = True
|
||||
current = prop
|
||||
fcurrent = fprop
|
||||
else:
|
||||
u = np.random.rand()
|
||||
if np.exp(fprop-fcurrent)>u:#sample accepted downhill
|
||||
accepted[it] = True
|
||||
current = prop
|
||||
fcurrent = fprop
|
||||
|
||||
#store current value
|
||||
if (it > Nburn) & ((it%Nthin)==0):
|
||||
self.chains[-1].append(current)
|
||||
|
||||
#tuning!
|
||||
if it & ((it%tune_interval)==0) & tune & ((it<Nburn) | (tune_throughout)):
|
||||
pc = np.mean(accepted[it-tune_interval:it])
|
||||
self.cov = np.cov(np.vstack(self.chains[-1][-tune_interval:]).T)
|
||||
if pc > .25:
|
||||
self.scale *= 1.1
|
||||
if pc < .15:
|
||||
self.scale /= 1.1
|
||||
|
||||
def predict(self,function,args):
|
||||
"""Make a prediction for the function, to which we will pass the additional arguments"""
|
||||
param = self.model._get_params()
|
||||
fs = []
|
||||
for p in self.chain:
|
||||
self.model._set_params(p)
|
||||
fs.append(function(*args))
|
||||
self.model._set_params(param)# reset model to starting state
|
||||
return fs
|
||||
Loading…
Add table
Add a link
Reference in a new issue