bug fix for mcmc sampler and add test case

This commit is contained in:
Zhenwen Dai 2016-03-10 18:37:53 +00:00
parent ba74e29aee
commit f2b813551a
2 changed files with 18 additions and 5 deletions

View file

@ -37,16 +37,14 @@ class Metropolis_Hastings(object):
def sample(self, Ntotal=10000, Nburn=1000, Nthin=10, tune=True, tune_throughout=False, tune_interval=400): def sample(self, Ntotal=10000, Nburn=1000, Nthin=10, tune=True, tune_throughout=False, tune_interval=400):
current = self.model.optimizer_array current = self.model.optimizer_array
fcurrent = self.model.log_likelihood() + self.model.log_prior() + \ fcurrent = self.model.log_likelihood() + self.model.log_prior()
self.model._log_det_jacobian()
accepted = np.zeros(Ntotal,dtype=np.bool) accepted = np.zeros(Ntotal,dtype=np.bool)
for it in range(Ntotal): for it in range(Ntotal):
print("sample %d of %d\r"%(it,Ntotal),end="\t") print("sample %d of %d\r"%(it+1,Ntotal),end="")
sys.stdout.flush() sys.stdout.flush()
prop = np.random.multivariate_normal(current, self.cov*self.scale*self.scale) prop = np.random.multivariate_normal(current, self.cov*self.scale*self.scale)
self.model.optimizer_array = prop self.model.optimizer_array = prop
fprop = self.model.log_likelihood() + self.model.log_prior() + \ fprop = self.model.log_likelihood() + self.model.log_prior()
self.model._log_det_jacobian()
if fprop>fcurrent:#sample accepted, going 'uphill' if fprop>fcurrent:#sample accepted, going 'uphill'
accepted[it] = True accepted[it] = True

View file

@ -65,5 +65,20 @@ class HMCSamplerTest(unittest.TestCase):
hmc = GPy.inference.mcmc.HMC(m,stepsize=1e-2) hmc = GPy.inference.mcmc.HMC(m,stepsize=1e-2)
s = hmc.sample(num_samples=3) s = hmc.sample(num_samples=3)
class MCMCSamplerTest(unittest.TestCase):
def test_sampling(self):
np.random.seed(1)
x = np.linspace(0.,2*np.pi,100)[:,None]
y = -np.cos(x)+np.random.randn(*x.shape)*0.3+1
m = GPy.models.GPRegression(x,y)
m.kern.lengthscale.set_prior(GPy.priors.Gamma.from_EV(1.,10.))
m.kern.variance.set_prior(GPy.priors.Gamma.from_EV(1.,10.))
m.likelihood.variance.set_prior(GPy.priors.Gamma.from_EV(1.,10.))
mcmc = GPy.inference.mcmc.Metropolis_Hastings(m)
mcmc.sample(Ntotal=100, Nburn=10)
if __name__ == "__main__": if __name__ == "__main__":
unittest.main() unittest.main()