Package csb :: Package statistics :: Package samplers :: Package mc :: Module singlechain
[frames] | no frames]

Source Code for Module csb.statistics.samplers.mc.singlechain

  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 
67 68 69 -class AbstractSingleChainMC(AbstractMC):
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.):
90 91 super(AbstractSingleChainMC, self).__init__(state) 92 93 self._pdf = pdf 94 self._temperature = temperature 95 self._nmoves = 0 96 self._accepted = 0 97 self._last_move_accepted = None
98
99 - def _checkstate(self, state):
100 if not isinstance(state, State): 101 raise TypeError(state)
102
103 - def sample(self):
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
128 - def _propose(self):
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
138 - def _calc_pacc(self, proposal_communicator):
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
149 - def _accept_proposal(self, proposal_state):
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
160 - def _update_statistics(self, accepted):
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
172 - def energy(self):
173 """ 174 Log-likelihood of the current state. 175 @rtype: float 176 """ 177 return self._pdf.log_prob(self.state.position)
178 179 @property
180 - def acceptance_rate(self):
181 """ 182 Acceptance rate. 183 """ 184 return float(self._accepted) / float(self._nmoves)
185 186 @property
187 - def last_move_accepted(self):
188 """ 189 Information whether the last MC move was accepted or not. 190 """ 191 return self._last_move_accepted
192 193 @property
194 - def temperature(self):
195 return self._temperature
196
197 198 -class HMCSampler(AbstractSingleChainMC):
199 """ 200 Hamilton Monte Carlo (HMC, also called Hybrid Monte Carlo by the inventors, 201 Duane, Kennedy, Pendleton, Duncan 1987). 202 203 @param pdf: Probability density function to be sampled from 204 @type pdf: L{csb.statistics.pdf.AbstractDensity} 205 206 @param state: Inital state 207 @type state: L{State} 208 209 @param gradient: Gradient of the negative log-probability 210 @type gradient: L{AbstractGradient} 211 212 @param timestep: Timestep used for integration 213 @type timestep: float 214 215 @param nsteps: Number of integration steps to be performed in 216 each iteration 217 @type nsteps: int 218 219 @param mass_matrix: Mass matrix 220 @type mass_matrix: n-dimensional L{InvertibleMatrix} with n being the dimension 221 of the configuration space, that is, the dimension of 222 the position / momentum vectors 223 224 @param integrator: Subclass of L{AbstractIntegrator} to be used for 225 integrating Hamiltionian equations of motion 226 @type integrator: L{AbstractIntegrator} 227 228 @param temperature: Pseudo-temperature of the Boltzmann ensemble 229 M{p(x) = 1/N * exp(-1/T * E(x))} with the 230 pseudo-energy defined as M{E(x) = -log(p(x))} 231 where M{p(x)} is the PDF under consideration 232 @type temperature: float 233 """ 234
235 - def __init__(self, pdf, state, gradient, timestep, nsteps, 236 mass_matrix=None, integrator=FastLeapFrog, temperature=1.):
237 238 super(HMCSampler, self).__init__(pdf, state, temperature) 239 240 self._timestep = None 241 self.timestep = timestep 242 self._nsteps = None 243 self.nsteps = nsteps 244 245 self._d = len(self.state.position) 246 247 self._mass_matrix = mass_matrix 248 if self.mass_matrix is None: 249 self.mass_matrix = InvertibleMatrix(numpy.eye(self._d), numpy.eye(self._d)) 250 251 self._momentum_covariance_matrix = self._temperature * self.mass_matrix 252 253 self._integrator = integrator 254 self._gradient = gradient 255 256 self._propagator = MDPropagator(self._gradient, self._timestep, 257 mass_matrix=self._mass_matrix, 258 integrator=self._integrator)
259
260 - def _propose(self):
267
268 - def _calc_pacc(self, proposal_communicator):
269 270 current_state = proposal_communicator.current_state 271 proposal_state = proposal_communicator.proposal_state 272 273 E = lambda x: -self._pdf.log_prob(x) 274 K = lambda x: 0.5 * numpy.dot(x.T, numpy.dot(self.mass_matrix.inverse, x)) 275 276 pacc = csb.numeric.exp((-K(proposal_state.momentum) - E(proposal_state.position) 277 + K(current_state.momentum) + E(current_state.position)) / 278 self.temperature) 279 280 if self.state.momentum is None: 281 proposal_communicator.proposal_state.momentum = None 282 else: 283 proposal_communicator.proposal_state.momentum = self.state.momentum 284 285 return pacc
286 287 @property
288 - def timestep(self):
289 return self._timestep
290 291 @timestep.setter
292 - def timestep(self, value):
293 self._timestep = float(value) 294 if "_propagator" in dir(self): 295 self._propagator.timestep = self._timestep
296 297 @property
298 - def nsteps(self):
299 return self._nsteps
300 301 @nsteps.setter
302 - def nsteps(self, value):
303 self._nsteps = int(value)
304 305 @property
306 - def mass_matrix(self):
307 return self._mass_matrix
308 @mass_matrix.setter
309 - def mass_matrix(self, value):
310 self._mass_matrix = value 311 if "_propagator" in dir(self): 312 self._propagator.mass_matrix = self._mass_matrix
313
314 315 -class RWMCSampler(AbstractSingleChainMC):
316 """ 317 Random Walk Metropolis Monte Carlo implementation 318 (Metropolis, Rosenbluth, Teller, Teller 1953; Hastings, 1970). 319 320 @param pdf: Probability density function to be sampled from 321 @type pdf: L{csb.statistics.pdf.AbstractDensity} 322 323 @param state: Inital state 324 @type state: L{State} 325 326 @param stepsize: Serves to set the step size in 327 proposal_density, e.g. for automatic acceptance 328 rate adaption 329 @type stepsize: float 330 331 @param proposal_density: The proposal density as a function f(x, s) 332 of the current state x and the stepsize s. 333 By default, the proposal density is uniform, 334 centered around x, and has width s. 335 @type proposal_density: callable 336 337 @param temperature: Pseudo-temperature of the Boltzmann ensemble 338 M{p(x) = 1/N * exp(-1/T * E(x))} with the 339 pseudo-energy defined as M{E(x) = -log(p(x))} 340 where M{p(x)} is the PDF under consideration 341 @type temperature: float 342 """ 343
344 - def __init__(self, pdf, state, stepsize=1., proposal_density=None, temperature=1.):
345 346 super(RWMCSampler, self).__init__(pdf, state, temperature) 347 self._stepsize = None 348 self.stepsize = stepsize 349 if proposal_density == None: 350 self._proposal_density = lambda x, s: x.position + \ 351 s * numpy.random.uniform(size=x.position.shape, low=-1., high=1.) 352 else: 353 self._proposal_density = proposal_density
354
355 - def _propose(self):
356 357 current_state = self.state.clone() 358 proposal_state = self.state.clone() 359 proposal_state.position = self._proposal_density(current_state, self.stepsize) 360 361 return SimpleProposalCommunicator(current_state, proposal_state)
362
363 - def _calc_pacc(self, proposal_communicator):
364 365 current_state = proposal_communicator.current_state 366 proposal_state = proposal_communicator.proposal_state 367 E = lambda x:-self._pdf.log_prob(x) 368 369 pacc = csb.numeric.exp((-(E(proposal_state.position) - E(current_state.position))) / 370 self.temperature) 371 return pacc
372 373 @property
374 - def stepsize(self):
375 return self._stepsize
376 377 @stepsize.setter
378 - def stepsize(self, value):
379 self._stepsize = float(value)
380
381 382 -class AbstractNCMCSampler(AbstractSingleChainMC):
383 """ 384 Implementation of the NCMC sampling algorithm (Nilmeier et al., "Nonequilibrium candidate Monte 385 Carlo is an efficient tool for equilibrium simulation", PNAS 2011) for sampling from one 386 ensemble only. 387 Subclasses have to specify the acceptance probability, which depends on the kind of 388 perturbations and propagations in the protocol. 389 390 @param state: Inital state 391 @type state: L{State} 392 393 @param protocol: Nonequilibrium protocol with alternating perturbations and propagations 394 @type protocol: L{Protocol} 395 396 @param reverse_protocol: The reversed version of the protocol, that is, the order of 397 perturbations and propagations in each step is reversed 398 @type reverse_protocol: L{Protocol} 399 """ 400 401 __metaclass__ = ABCMeta 402
403 - def __init__(self, state, protocol, reverse_protocol):
404 405 self._protocol = None 406 self.protocol = protocol 407 self._reverse_protocol = None 408 self.reverse_protocol = reverse_protocol 409 410 pdf = self.protocol.steps[0].perturbation.sys_before.hamiltonian 411 temperature = self.protocol.steps[0].perturbation.sys_before.hamiltonian.temperature 412 413 super(AbstractNCMCSampler, self).__init__(pdf, state, temperature)
414
415 - def _pick_protocol(self):
416 """ 417 Picks either the protocol or the reversed protocol with equal probability. 418 419 @return: Either the protocol or the reversed protocol 420 @rtype: L{Protocol} 421 """ 422 423 if numpy.random.random() < 0.5: 424 return self.protocol 425 else: 426 return self.reverse_protocol
427
428 - def _propose(self):
429 430 protocol = self._pick_protocol() 431 432 gen = NonequilibriumStepPropagator(protocol) 433 434 traj = gen.generate(self.state) 435 436 return NCMCProposalCommunicator(traj)
437
438 - def _accept_proposal(self, proposal_state):
439 440 if self.state.momentum is not None: 441 proposal_state.momentum *= -1.0 442 else: 443 proposal_state.momentum = None 444 445 super(AbstractNCMCSampler, self)._accept_proposal(proposal_state)
446 447 @property
448 - def protocol(self):
449 return self._protocol
450 @protocol.setter
451 - def protocol(self, value):
452 self._protocol = value
453 454 @property
455 - def reverse_protocol(self):
456 return self._reverse_protocol
457 @reverse_protocol.setter
458 - def reverse_protocol(self, value):
459 self._reverse_protocol = value
460
461 462 -class SimpleProposalCommunicator(object):
463 """ 464 This holds all the information needed to calculate the acceptance 465 probability in both the L{RWMCSampler} and L{HMCSampler} classes, 466 that is, only the proposal state. 467 For more advanced algorithms, one may derive classes capable of 468 holding the neccessary additional information from this class. 469 470 @param current_state: Current state 471 @type current_state: L{State} 472 473 @param proposal_state: Proposal state 474 @type proposal_state: L{State} 475 """ 476 477 __metaclass__ = ABCMeta 478
479 - def __init__(self, current_state, proposal_state):
480 481 self._current_state = current_state 482 self._proposal_state = proposal_state
483 484 @property
485 - def current_state(self):
486 return self._current_state
487 488 @property
489 - def proposal_state(self):
490 return self._proposal_state
491
492 493 -class NCMCProposalCommunicator(SimpleProposalCommunicator):
494 """ 495 Holds all information (that is, the trajectory with heat, work, Hamiltonian difference 496 and jacbian) needed to calculate the acceptance probability in the AbstractNCMCSampler class. 497 498 @param traj: Non-equilibrium trajectory stemming from a stepwise protocol 499 @type traj: NCMCTrajectory 500 """ 501
502 - def __init__(self, traj):
503 504 self._traj = None 505 self.traj = traj 506 507 super(NCMCProposalCommunicator, self).__init__(traj.initial, traj.final)
508