1 """
2 Various Monte Carlo equilibrium sampling algorithms, which simulate only one Markov chain.
3
4 Here is how to sample from a PDF using the L{HMCSampler} class. In the following
5 snippet we draw 5000 samples from a 1D normal distribution and plot them:
6
7
8 >>> import numpy
9 >>> from csb.io.plots import Chart
10 >>> from csb.statistics.pdf import Normal
11 >>> from csb.statistics.samplers import State
12 >>> from csb.statistics.samplers.mc.singlechain import HMCSampler
13
14 >>> initial_state = State(numpy.array([1.]))
15 >>> grad = lambda q, t: q
16 >>> timestep = 1.5
17 >>> nsteps = 30
18 >>> nsamples = 5000
19
20 >>> sampler = HMCSampler(Normal(), initial_state, grad, timestep, nsteps)
21
22 >>> states = []
23 >>> for i in range(nsamples):
24 sampler.sample()
25 states.append(sampler.state)
26
27 >>> print('acceptance rate:', sampler.acceptance_rate)
28 0.8
29
30 >>> states = [state.position[0]for state in states]
31 >>> chart = Chart()
32 >>> chart.plot.hist([numpy.random.normal(size=5000), states], bins=20, normed=True)
33 >>> chart.plot.legend(['numpy.random.normal', 'HMC'])
34 >>> chart.show()
35
36
37 First, several things which are being needed are imported.
38 As every sampler in this module implements a Markov Chain, an initial state has to be
39 chosen. In the following lines, several parameters are set:
40 - the gradient of the negative log-probability of the PDF under consideration
41 - the integration timestep
42 - the number of integration steps to be performed in each iteration, that is, the HMC
43 trajectory length
44 - the number of samples to be drawn
45
46 The empty list states is initialized. It will serve to store the samples drawn.
47 In the loop, C{sampler.sample()} is repeatedly called. After each call of C{sampler.sample()},
48 the current state of the Markov Chain is stored in sampler.state and this state is appended
49 to the sample storage list.
50
51 Then the acceptance rate is printed, the numeric values are being extracted from the
52 L{State} objects in states, a histogram is created and finally plotted.
53 """
54
55 import numpy
56 import csb.numeric
57 import csb.core
58
59 from abc import ABCMeta, abstractmethod
60
61 from csb.statistics.samplers import State
62 from csb.statistics.samplers.mc import AbstractMC, MCCollection, augment_state
63 from csb.statistics.samplers.mc.propagators import MDPropagator
64 from csb.statistics.samplers.mc.neqsteppropagator import NonequilibriumStepPropagator
65 from csb.numeric.integrators import FastLeapFrog
66 from csb.numeric import InvertibleMatrix
70 """
71 Abstract class for Monte Carlo sampling algorithms simulating
72 only one ensemble.
73
74 @param pdf: probability density function to sample from
75 @type pdf: subclass of L{csb.statistics.pdf.AbstractDensity}
76
77 @param state: Initial state
78 @type state: L{State}
79
80 @param temperature: Pseudo-temperature of the Boltzmann ensemble
81 M{p(x) = 1/N * exp(-1/T * E(x))} with the
82 pseudo-energy defined as M{E(x) = -log(p(x))}
83 where M{p(x)} is the PDF under consideration
84 @type temperature: float
85 """
86
87 __metaclass__ = ABCMeta
88
89 - def __init__(self, pdf, state, temperature=1.):
98
102
104 """
105 Draw a sample.
106
107 @rtype: L{State}
108 """
109
110 proposal_communicator = self._propose()
111 pacc = self._calc_pacc(proposal_communicator)
112
113 accepted = None
114 if numpy.random.random() < pacc:
115 accepted = True
116 else:
117 accepted = False
118
119 if accepted == True:
120 self._accept_proposal(proposal_communicator.proposal_state)
121
122 self._update_statistics(accepted)
123 self._last_move_accepted = accepted
124
125 return self.state
126
127 @abstractmethod
129 """
130 Calculate a new proposal state and gather additional information
131 needed to calculate the acceptance probability.
132
133 @rtype: L{SimpleProposalCommunicator}
134 """
135 pass
136
137 @abstractmethod
139 """
140 Calculate probability with which to accept the proposal.
141
142 @param proposal_communicator: Contains information about the proposal
143 and additional information needed to
144 calculate the acceptance probability
145 @type proposal_communicator: L{SimpleProposalCommunicator}
146 """
147 pass
148
150 """
151 Accept the proposal state by setting it as the current state of the sampler
152 object
153
154 @param proposal_state: The proposal state
155 @type proposal_state: L{State}
156 """
157
158 self.state = proposal_state
159
161 """
162 Update the sampling statistics.
163
164 @param accepted: Whether or not the proposal state has been accepted
165 @type accepted: boolean
166 """
167
168 self._nmoves += 1
169 self._accepted += int(accepted)
170
171 @property
173 """
174 Negative log-likelihood of the current state.
175 @rtype: float
176 """
177 return -self._pdf.log_prob(self.state.position)
178
179 @property
181 """
182 Acceptance rate.
183 """
184
185 if self._nmoves > 0:
186 return float(self._accepted) / float(self._nmoves)
187 else:
188 return 0.0
189
190 @property
192 """
193 Information whether the last MC move was accepted or not.
194 """
195 return self._last_move_accepted
196
197 @property
199 return self._temperature
200
203 """
204 Hamilton Monte Carlo (HMC, also called Hybrid Monte Carlo by the inventors,
205 Duane, Kennedy, Pendleton, Duncan 1987).
206
207 @param pdf: Probability density function to be sampled from
208 @type pdf: L{csb.statistics.pdf.AbstractDensity}
209
210 @param state: Inital state
211 @type state: L{State}
212
213 @param gradient: Gradient of the negative log-probability
214 @type gradient: L{AbstractGradient}
215
216 @param timestep: Timestep used for integration
217 @type timestep: float
218
219 @param nsteps: Number of integration steps to be performed in
220 each iteration
221 @type nsteps: int
222
223 @param mass_matrix: Mass matrix
224 @type mass_matrix: n-dimensional L{InvertibleMatrix} with n being the dimension
225 of the configuration space, that is, the dimension of
226 the position / momentum vectors
227
228 @param integrator: Subclass of L{AbstractIntegrator} to be used for
229 integrating Hamiltionian equations of motion
230 @type integrator: L{AbstractIntegrator}
231
232 @param temperature: Pseudo-temperature of the Boltzmann ensemble
233 M{p(x) = 1/N * exp(-1/T * E(x))} with the
234 pseudo-energy defined as M{E(x) = -log(p(x))}
235 where M{p(x)} is the PDF under consideration
236 @type temperature: float
237 """
238
239 - def __init__(self, pdf, state, gradient, timestep, nsteps,
240 mass_matrix=None, integrator=FastLeapFrog, temperature=1.):
256
266
268 """
269 Factory which produces a L{MDPropagator} object initialized with
270 the MD parameters given in __init__().
271
272 @return: L{MDPropagator} instance
273 @rtype: L{MDPropagator}
274 """
275
276 return MDPropagator(self._gradient, self._timestep,
277 mass_matrix=self._mass_matrix,
278 integrator=self._integrator)
279
287
289 """
290 Evaluates the Hamiltonian consisting of the negative log-probability
291 and a quadratic kinetic term.
292
293 @param state: State on which the Hamiltonian should be evaluated
294 @type state: L{State}
295
296 @return: Value of the Hamiltonian (total energy)
297 @rtype: float
298 """
299
300 V = lambda q: -self._pdf.log_prob(q)
301 T = lambda p: 0.5 * numpy.dot(p.T, numpy.dot(self.mass_matrix.inverse, p))
302
303 return T(state.momentum) + V(state.position)
304
319
320 @property
322 return self._timestep
323
324 @timestep.setter
326 self._timestep = float(value)
327 if "_propagator" in dir(self):
328 self._propagator.timestep = self._timestep
329
330 @property
333
334 @nsteps.setter
336 self._nsteps = int(value)
337
338 @property
340 return self._mass_matrix
341 @mass_matrix.setter
343 self._mass_matrix = value
344 if "_propagator" in dir(self):
345 self._propagator.mass_matrix = self._mass_matrix
346
349 """
350 Random Walk Metropolis Monte Carlo implementation
351 (Metropolis, Rosenbluth, Teller, Teller 1953; Hastings, 1970).
352
353 @param pdf: Probability density function to be sampled from
354 @type pdf: L{csb.statistics.pdf.AbstractDensity}
355
356 @param state: Inital state
357 @type state: L{State}
358
359 @param stepsize: Serves to set the step size in
360 proposal_density, e.g. for automatic acceptance
361 rate adaption
362 @type stepsize: float
363
364 @param proposal_density: The proposal density as a function f(x, s)
365 of the current state x and the stepsize s.
366 By default, the proposal density is uniform,
367 centered around x, and has width s.
368 @type proposal_density: callable
369
370 @param temperature: Pseudo-temperature of the Boltzmann ensemble
371 M{p(x) = 1/N * exp(-1/T * E(x))} with the
372 pseudo-energy defined as M{E(x) = -log(p(x))}
373 where M{p(x)} is the PDF under consideration
374 @type temperature: float
375 """
376
377 - def __init__(self, pdf, state, stepsize=1., proposal_density=None, temperature=1.):
388
396
407
408 @property
410 return self._stepsize
411
412 @stepsize.setter
414 self._stepsize = float(value)
415
418 """
419 Implementation of the NCMC sampling algorithm (Nilmeier et al., "Nonequilibrium candidate Monte
420 Carlo is an efficient tool for equilibrium simulation", PNAS 2011) for sampling from one
421 ensemble only.
422 Subclasses have to specify the acceptance probability, which depends on the kind of
423 perturbations and propagations in the protocol.
424
425 @param state: Inital state
426 @type state: L{State}
427
428 @param protocol: Nonequilibrium protocol with alternating perturbations and propagations
429 @type protocol: L{Protocol}
430
431 @param reverse_protocol: The reversed version of the protocol, that is, the order of
432 perturbations and propagations in each step is reversed
433 @type reverse_protocol: L{Protocol}
434 """
435
436 __metaclass__ = ABCMeta
437
438 - def __init__(self, state, protocol, reverse_protocol):
449
451 """
452 Picks either the protocol or the reversed protocol with equal probability.
453
454 @return: Either the protocol or the reversed protocol
455 @rtype: L{Protocol}
456 """
457
458 if numpy.random.random() < 0.5:
459 return self.protocol
460 else:
461 return self.reverse_protocol
462
472
481
482 @property
484 return self._protocol
485 @protocol.setter
487 self._protocol = value
488
489 @property
491 return self._reverse_protocol
492 @reverse_protocol.setter
494 self._reverse_protocol = value
495
498 """
499 This holds all the information needed to calculate the acceptance
500 probability in both the L{RWMCSampler} and L{HMCSampler} classes,
501 that is, only the proposal state.
502 For more advanced algorithms, one may derive classes capable of
503 holding the neccessary additional information from this class.
504
505 @param current_state: Current state
506 @type current_state: L{State}
507
508 @param proposal_state: Proposal state
509 @type proposal_state: L{State}
510 """
511
512 __metaclass__ = ABCMeta
513
514 - def __init__(self, current_state, proposal_state):
518
519 @property
521 return self._current_state
522
523 @property
525 return self._proposal_state
526
529 """
530 Holds all information (that is, the trajectory with heat, work, Hamiltonian difference
531 and jacbian) needed to calculate the acceptance probability in the AbstractNCMCSampler class.
532
533 @param traj: Non-equilibrium trajectory stemming from a stepwise protocol
534 @type traj: NCMCTrajectory
535 """
536
543