Mercurial > pylearn

annotate pylearn/sampling/tests/test_mcmc.py @ 1447:fbe470217937

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Use .get_value() and .set_value() of shared instead of the .value property
author Pascal Lamblin <lamblinp@iro.umontreal.ca>
date Wed, 16 Mar 2011 20:20:02 -0400
parents 492473059b37
children
rev   line source
955
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
1 from pylearn.sampling.mcmc import *
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
2
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
3 def _sampler_on_2d_gaussian(sampler_cls, burnin, n_samples):
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
4 batchsize=3
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
5
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
6 rng = np.random.RandomState(234)
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
7
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
8 #
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
9 # Define a covariance and mu for a gaussian
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
10 #
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
11 tmp = rng.randn(2,2)
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
12 tmp[0] += tmp[1] #induce some covariance
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
13 cov = np.dot(tmp, tmp.T)
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
14 cov_inv = np.linalg.inv(cov)
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
15 mu = np.asarray([5, 9.5], dtype=theano.config.floatX)
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
16
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
17 def gaussian_energy(xlist):
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
18 x, = xlist
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
19 return 0.5 * (TT.dot((x-mu),cov_inv)*(x-mu)).sum(axis=1)
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
20
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
21
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
22 position = shared(rng.randn(batchsize, 2).astype(theano.config.floatX))
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
23 sampler = sampler_cls([position], gaussian_energy)
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
24
1447
fbe470217937 Use .get_value() and .set_value() of shared instead of the .value property
Pascal Lamblin <lamblinp@iro.umontreal.ca>
parents: 955
diff changeset
25 print 'initial position', position.get_value(borrow=True)
955
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
26 print 'initial stepsize', sampler.stepsize
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
27
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
28 # DRAW SAMPLES
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
29
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
30 samples = [sampler.draw() for r in xrange(burnin)] #burn-in
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
31 samples = np.asarray([sampler.draw() for r in xrange(n_samples)])
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
32
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
33 assert sampler.avg_acceptance_rate > 0
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
34 assert sampler.avg_acceptance_rate < 1
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
35
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
36 # TEST THAT THEY ARE FROM THE RIGHT DISTRIBUTION
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
37
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
38 # samples.shape == (1000, 1, 3, 2)
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
39
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
40 print 'target mean:', mu
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
41 print 'empirical mean: ', samples.mean(axis=0)[0]
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
42 #assert np.all(abs(mu - samples.mean(axis=0)) < 1)
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
43
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
44
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
45 print 'final stepsize', sampler.stepsize
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
46 print 'final acceptance_rate', sampler.avg_acceptance_rate
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
47
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
48 print 'target cov', cov
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
49 s = samples[:,0,0,:]
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
50 empirical_cov = np.cov(samples[:,0,0,:].T)
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
51 print ''
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
52 print 'cov/empirical_cov', cov/empirical_cov
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
53 empirical_cov = np.cov(samples[:,0,1,:].T)
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
54 print 'cov/empirical_cov', cov/empirical_cov
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
55 empirical_cov = np.cov(samples[:,0,2,:].T)
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
56 print 'cov/empirical_cov', cov/empirical_cov
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
57 return sampler
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
58
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
59 def test_mcmc():
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
60 print ('MCMC')
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
61 sampler = _sampler_on_2d_gaussian(MCMC_sampler, burnin=3000, n_samples=90000)