Mercurial > pylearn

annotate pylearn/sampling/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 """Metropolis-Hastings Monte Carlo
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
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
4 See also Hamiltonian Monte Carlo (aka Hybrid Monte Carlo) in the hmc.py file.
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 """
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
7 import sys
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
8 import logging
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
9 import numpy as np
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
10 from theano import function, shared
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
11 from theano import tensor as TT
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
12 import theano.sparse #installs the sparse shared var handler
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
13 from theano.printing import Print
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
14
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
15 #TODO:
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
16 # Consider how to refactor this code so that the key functions /functionality are provided as
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
17 # theano expressions??
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
18 # It could be limiting that the implementation requires that the position be a list of shared
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
19 # variables, and basically uses procedural logic on top of that.
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 # Theano can only do function optimization on one function at a time, so until these are
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
22 # written as update expressions, they must be compiled as their own functions and cannot be
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
23 # rolled into things the user is doing anyway.
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
24
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
25 class MCMC_sampler(object):
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
26 """Generic Metropolis-Hastings Monte-Carlo sampler
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 The `draw` function advances the markov chain and returns the current sample by calling
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
29 `simulate` and `get_position` in sequence.
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
30
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
31 The stepsize is updated dynamically to
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 """
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
34
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
35 # TODO: anyone has better guesses for these constants??
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
36 stepsize_dec = 0.98
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
37 stepsize_min = 0.001
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
38 stepsize_max = 0.25
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
39 stepsize_inc = 1.02
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
40 target_acceptance_rate=0.7
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
41
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
42 avg_acceptance_slowness = 0.9 # used in geometric avg. 1.0 would be not moving at all
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
43 n_steps=1
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 def __init__(self, positions, energy_fn,
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
46 initial_stepsize=0.01,
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
47 seed=12345):
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
48 """
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
49 :param positions: list of shared ndarray variables.
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
50
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
51 :param energy:
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
52
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
53 callable such that energy_fn(positions)
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
54 returns theano vector of energies.
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
55 The len of this vector is the batchsize.
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
56 """
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
57
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
58 batchsize = positions[0].get_value(borrow=True).shape[0]
955
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
59 self.s_rng = TT.shared_randomstreams.RandomStreams(seed)
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
60 self.positions = positions
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
61 self.prev_energy = shared(np.zeros(batchsize) + float('inf'))
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
62 self.avg_acceptance_rate = 0.5
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
63 self.stepsize = initial_stepsize
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
64
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
65 s_stepsize = TT.scalar('stepsize')
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
66
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
67 new_positions = [p + s_stepsize * self.s_rng.normal(size=p.get_value(borrow=True).shape)
955
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
68 for p in self.positions]
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
69
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
70 # accept-reject according to Metropolis-Hastings
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
71
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
72 energy = energy_fn(new_positions)
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
73 accept = TT.exp(self.prev_energy - energy) - self.s_rng.uniform(size=(batchsize,)) >= 0
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
74
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
75 self.accept_reject_positions = function([s_stepsize], accept.mean(),
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
76 updates=[(self.prev_energy, TT.switch(accept, energy, self.prev_energy))] + [
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
77 (q, TT.switch(accept.dimshuffle(0, *(('x',)*(q.ndim-1))), new_q, q))
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
78 for q, new_q in zip(self.positions, new_positions)])
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
79
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
80 def simulate(self, n_steps=None):
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
81 if n_steps is None:
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
82 n_steps = self.n_steps
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
83 for ss in range(n_steps):
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
84 acceptance_rate = self.accept_reject_positions(self.stepsize)
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
85 self.avg_acceptance_rate = self.avg_acceptance_slowness * self.avg_acceptance_rate \
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
86 + (1.0 - self.avg_acceptance_slowness) * acceptance_rate
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
87 if self.avg_acceptance_rate < self.target_acceptance_rate:
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
88 self.stepsize = max(self.stepsize*self.stepsize_dec,self.stepsize_min)
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
89 else:
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
90 self.stepsize = min(self.stepsize*self.stepsize_inc,self.stepsize_max)
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
91
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
92 def get_position(self):
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
93 return [q.get_value(borrow=True) for q in self.positions]
955
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
94
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
95 def draw(self, n_steps=None):
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
96 """Return the current sample in the Markov chain as a list of numpy arrays
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
97
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
98 The size of the arrays will match the size of the `initial_position` argument to
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
99 __init__.
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
100 """
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
101 self.simulate(n_steps=n_steps)
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
102 return self.get_position()
492473059b37 Adding sampling module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
103