Mercurial > pylearn

comparison pylearn/sampling/tests/test_hmc.py @ 1526:5804e44d7a1b

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
pep8
author Frederic Bastien <nouiz@nouiz.org>
date Fri, 09 Nov 2012 16:58:17 -0500
parents 1ee532a6f33b
children 211b217f9691
comparison
equal deleted inserted replaced
1525:9c24a2bdbe90 1526:5804e44d7a1b
2 import theano 2 import theano
3 from theano import tensor 3 from theano import tensor
4 4
5 from pylearn.sampling.hmc import HMC_sampler 5 from pylearn.sampling.hmc import HMC_sampler
6 6
7
7 def _sampler_on_2d_gaussian(sampler_cls, burnin, n_samples): 8 def _sampler_on_2d_gaussian(sampler_cls, burnin, n_samples):
8 batchsize=3 9 batchsize = 3
9 10
10 rng = numpy.random.RandomState(234) 11 rng = numpy.random.RandomState(234)
11 12
12 # 13 #
13 # Define a covariance and mu for a gaussian 14 # Define a covariance and mu for a gaussian
14 # 15 #
15 tmp = rng.randn(2,2).astype(theano.config.floatX) 16 tmp = rng.randn(2, 2).astype(theano.config.floatX)
16 tmp[0] += tmp[1] #induce some covariance 17 tmp[0] += tmp[1] # induce some covariance
17 cov = numpy.dot(tmp, tmp.T) 18 cov = numpy.dot(tmp, tmp.T)
18 cov_inv = numpy.linalg.inv(cov).astype(theano.config.floatX) 19 cov_inv = numpy.linalg.inv(cov).astype(theano.config.floatX)
19 mu = numpy.asarray([5, 9.5], dtype=theano.config.floatX) 20 mu = numpy.asarray([5, 9.5], dtype=theano.config.floatX)
20 21
21 def gaussian_energy(xlist): 22 def gaussian_energy(xlist):
22 x = xlist 23 x = xlist
23 return 0.5 * (tensor.dot((x-mu),cov_inv)*(x-mu)).sum(axis=1) 24 return 0.5 * (tensor.dot((x - mu), cov_inv) * (x - mu)).sum(axis=1)
24 25
25 26 position = theano.shared(rng.randn(batchsize,
26 position = theano.shared(rng.randn(batchsize, 2).astype(theano.config.floatX)) 27 2).astype(theano.config.floatX))
27 sampler = sampler_cls(position, gaussian_energy) 28 sampler = sampler_cls(position, gaussian_energy)
28 29
29 print 'initial position', position.get_value(borrow=True) 30 print 'initial position', position.get_value(borrow=True)
30 print 'initial stepsize', sampler.stepsize.get_value(borrow=True) 31 print 'initial stepsize', sampler.stepsize.get_value(borrow=True)
31 32
32 # DRAW SAMPLES 33 # DRAW SAMPLES
33 34 samples = [sampler.draw() for r in xrange(burnin)] # burn-in
34 samples = [sampler.draw() for r in xrange(burnin)] #burn-in
35 samples = numpy.asarray([sampler.draw() for r in xrange(n_samples)]) 35 samples = numpy.asarray([sampler.draw() for r in xrange(n_samples)])
36 36
37 assert sampler.avg_acceptance_rate.get_value() > 0 37 assert sampler.avg_acceptance_rate.get_value() > 0
38 assert sampler.avg_acceptance_rate.get_value() < 1 38 assert sampler.avg_acceptance_rate.get_value() < 1
39 39
43 43
44 print 'target mean:', mu 44 print 'target mean:', mu
45 print 'empirical mean: ', samples.mean(axis=0) 45 print 'empirical mean: ', samples.mean(axis=0)
46 #assert numpy.all(abs(mu - samples.mean(axis=0)) < 1) 46 #assert numpy.all(abs(mu - samples.mean(axis=0)) < 1)
47 47
48
49 print 'final stepsize', sampler.stepsize.get_value() 48 print 'final stepsize', sampler.stepsize.get_value()
50 print 'final acceptance_rate', sampler.avg_acceptance_rate.get_value() 49 print 'final acceptance_rate', sampler.avg_acceptance_rate.get_value()
51 50
52 print 'target cov', cov 51 print 'target cov', cov
53 s = samples[:,0,:] 52 s = samples[:, 0, :]
54 empirical_cov = numpy.cov(samples[:,0,:].T) 53 empirical_cov = numpy.cov(samples[:, 0, :].T)
55 print '' 54 print ''
56 print 'cov/empirical_cov', cov/empirical_cov 55 print 'cov/empirical_cov', cov / empirical_cov
57 empirical_cov = numpy.cov(samples[:,1,:].T) 56 empirical_cov = numpy.cov(samples[:, 1, :].T)
58 print 'cov/empirical_cov', cov/empirical_cov 57 print 'cov/empirical_cov', cov / empirical_cov
59 empirical_cov = numpy.cov(samples[:,2,:].T) 58 empirical_cov = numpy.cov(samples[:, 2, :].T)
60 print 'cov/empirical_cov', cov/empirical_cov 59 print 'cov/empirical_cov', cov / empirical_cov
61 return sampler 60 return sampler
61
62 62
63 def test_hmc(): 63 def test_hmc():
64 print ('HMC') 64 print ('HMC')
65 sampler = _sampler_on_2d_gaussian(HMC_sampler.new_from_shared_positions, burnin=3000/20, n_samples=90000/20) 65 sampler = _sampler_on_2d_gaussian(HMC_sampler.new_from_shared_positions,
66 burnin=3000 / 20, n_samples=90000 / 20)
66 assert abs(sampler.avg_acceptance_rate.get_value() - sampler.target_acceptance_rate) < .1 67 assert abs(sampler.avg_acceptance_rate.get_value() - sampler.target_acceptance_rate) < .1
67 assert sampler.stepsize.get_value() >= sampler.stepsize_min 68 assert sampler.stepsize.get_value() >= sampler.stepsize_min
68 assert sampler.stepsize.get_value() <= sampler.stepsize_max 69 assert sampler.stepsize.get_value() <= sampler.stepsize_max
69 70