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          Negative 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
    185          if self._nmoves > 0:
    186              return float(self._accepted) / float(self._nmoves)
    187          else:
    188              return 0.0
    189
    190      @property
    191 -    def last_move_accepted(self):
    192          """
    193          Information whether the last MC move was accepted or not.
    194          """
    195          return self._last_move_accepted
    196
    197      @property
    198 -    def temperature(self):
    199          return self._temperature
    200
    201
    202 -class HMCSampler(AbstractSingleChainMC):
    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.):
    241
    242          super(HMCSampler, self).__init__(pdf, state, temperature)
    243
    244          self._timestep = None
    245          self.timestep = timestep
    246          self._nsteps = None
    247          self.nsteps = nsteps
    248          self._mass_matrix = None
    249          self._momentum_covariance_matrix = None
    250          self._integrator = integrator
    251          self._gradient = gradient
    252
    253          self._setup_matrices(mass_matrix)
    254
    255          self._propagator = self._propagator_factory()
    256
    257 -    def _setup_matrices(self, mass_matrix):
    258
    259          self._d = len(self.state.position)
    260
    261          self._mass_matrix = mass_matrix
    262          if self.mass_matrix is None:
    263              self.mass_matrix = InvertibleMatrix(numpy.eye(self._d), numpy.eye(self._d))
    264
    265          self._momentum_covariance_matrix = self._temperature * self.mass_matrix
    266
    267 -    def _propagator_factory(self):
    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
    280 -    def _propose(self):
    281
    282          current_state = self.state.clone()
    283          current_state = augment_state(current_state, self.temperature, self.mass_matrix)
    284          proposal_state = self._propagator.generate(current_state, self._nsteps).final
    285
    286          return SimpleProposalCommunicator(current_state, proposal_state)
    287
    288 -    def _hamiltonian(self, state):
    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
    305 -    def _calc_pacc(self, proposal_communicator):
    306
    307          cs = proposal_communicator.current_state
    308          ps = proposal_communicator.proposal_state
    309
    310          pacc = csb.numeric.exp(-(self._hamiltonian(ps) - self._hamiltonian(cs)) /
    311                                 self.temperature)
    312
    313          if self.state.momentum is None:
    314              proposal_communicator.proposal_state.momentum = None
    315          else:
    316              proposal_communicator.proposal_state.momentum = self.state.momentum
    317
    318          return pacc
    319
    320      @property
    321 -    def timestep(self):
    322          return self._timestep
    323
    324      @timestep.setter
    325 -    def timestep(self, value):
    326          self._timestep = float(value)
    327          if "_propagator" in dir(self):
    328              self._propagator.timestep = self._timestep
    329
    330      @property
    331 -    def nsteps(self):
    332          return self._nsteps
    333
    334      @nsteps.setter
    335 -    def nsteps(self, value):
    336          self._nsteps = int(value)
    337
    338      @property
    339 -    def mass_matrix(self):
    340          return self._mass_matrix
    341      @mass_matrix.setter
    342 -    def mass_matrix(self, value):
    343          self._mass_matrix = value
    344          if "_propagator" in dir(self):
    345             self._propagator.mass_matrix = self._mass_matrix
    346
    347
    348 -class RWMCSampler(AbstractSingleChainMC):
    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.):
    378
    379          super(RWMCSampler, self).__init__(pdf, state, temperature)
    380          self._stepsize = None
    381          self.stepsize = stepsize
    382          if proposal_density == None:
    383              self._proposal_density = lambda x, s: x.position + \
    384                                       s * numpy.random.uniform(size=x.position.shape,
    385                                                                low=-1., high=1.)
    386          else:
    387              self._proposal_density = proposal_density
    388
    389 -    def _propose(self):
    390
    391          current_state = self.state.clone()
    392          proposal_state = self.state.clone()
    393          proposal_state.position = self._proposal_density(current_state, self.stepsize)
    394
    395          return SimpleProposalCommunicator(current_state, proposal_state)
    396
    397 -    def _calc_pacc(self, proposal_communicator):
    398
    399          current_state = proposal_communicator.current_state
    400          proposal_state = proposal_communicator.proposal_state
    401          E = lambda x: -self._pdf.log_prob(x)
    402
    403          pacc = csb.numeric.exp((-(E(proposal_state.position) - E(current_state.position))) /
    404                                 self.temperature)
    405
    406          return pacc
    407
    408      @property
    409 -    def stepsize(self):
    410          return self._stepsize
    411
    412      @stepsize.setter
    413 -    def stepsize(self, value):
    414          self._stepsize = float(value)
    415
    416
    417 -class AbstractNCMCSampler(AbstractSingleChainMC):
    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):
    439
    440          self._protocol = None
    441          self.protocol = protocol
    442          self._reverse_protocol = None
    443          self.reverse_protocol = reverse_protocol
    444
    445          pdf = self.protocol.steps[0].perturbation.sys_before.hamiltonian
    446          temperature = self.protocol.steps[0].perturbation.sys_before.hamiltonian.temperature
    447
    448          super(AbstractNCMCSampler, self).__init__(pdf, state, temperature)
    449
    450 -    def _pick_protocol(self):
    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
    463 -    def _propose(self):
    464
    465          protocol = self._pick_protocol()
    466
    467          gen = NonequilibriumStepPropagator(protocol)
    468
    469          traj = gen.generate(self.state)
    470
    471          return NCMCProposalCommunicator(traj)
    472
    473 -    def _accept_proposal(self, proposal_state):
    474
    475          if self.state.momentum is not None:
    476              proposal_state.momentum *= -1.0
    477          else:
    478              proposal_state.momentum = None
    479
    480          super(AbstractNCMCSampler, self)._accept_proposal(proposal_state)
    481
    482      @property
    483 -    def protocol(self):
    484          return self._protocol
    485      @protocol.setter
    486 -    def protocol(self, value):
    487          self._protocol = value
    488
    489      @property
    490 -    def reverse_protocol(self):
    491          return self._reverse_protocol
    492      @reverse_protocol.setter
    493 -    def reverse_protocol(self, value):
    494          self._reverse_protocol = value
    495
    496
    497 -class SimpleProposalCommunicator(object):
    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):
    515
    516          self._current_state = current_state
    517          self._proposal_state = proposal_state
    518
    519      @property
    520 -    def current_state(self):
    521          return self._current_state
    522
    523      @property
    524 -    def proposal_state(self):
    525          return self._proposal_state
    526
    527
    528 -class NCMCProposalCommunicator(SimpleProposalCommunicator):
    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
    537 -    def __init__(self, traj):
    538
    539          self._traj = None
    540          self.traj = traj
    541
    542          super(NCMCProposalCommunicator, self).__init__(traj.initial, traj.final)
    543