Tutorial: Monte Carlo Sampling

Single Chain

Various Monte Carlo equilibrium sampling algorithms, which simulate only one Markov chain, are defined in csb.statistics.samplers.mc.singlechain.

Here is how to sample from a PDF using the HMCSampler class. In the following snippet we draw 5000 samples from a 1D normal distribution and plot them:

>>> import numpy
>>> from csb.io.plots import Chart
>>> from csb.statistics.pdf import Normal
>>> from csb.statistics.samplers import State
>>> from csb.statistics.samplers.mc.singlechain import HMCSampler
>>> initial_state = State(numpy.array([1.](1.)))
>>> grad = lambda q, t: q
>>> timestep = 1.5
>>> nsteps = 30
>>> nsamples = 5000
>>> sampler = HMCSampler(Normal(), initial_state, grad, timestep, nsteps)
>>> states = []()
>>> for i in range(nsamples):
        sampler.sample()
        states.append(sampler.state)
>>> print('acceptance rate:', sampler.acceptance_rate)
0.8
>>> states = [state.position[0](state.position[0)for state in states]
>>> chart = Chart()
>>> chart.plot.hist([numpy.random.normal(size=5000), states](numpy.random.normal(size=5000),-states), bins=20, normed=True)
>>> chart.plot.legend(['numpy.random.normal', 'HMC']('numpy.random.normal',-'HMC'))
>>> chart.show()

First, several things which are being needed are imported. As every sampler in this module implements a Markov Chain, an initial state has to be chosen. In the following lines, several parameters are set:

the gradient of the negative log-probability of the PDF under consideration
the integration timestep
the number of integration steps to be performed in each iteration, that is, the HMC trajectory length
the number of samples to be drawn

The empty list states is initialized. It will serve to store the samples drawn. In the loop, sampler.sample() is repeatedly called. After each call of sampler.sample(), the current state of the Markov Chain is stored in sampler.state and this state is appended to the sample storage list.

Then the acceptance rate is printed, the numeric values are being extracted from the State objects in states, a histogram is created and finally plotted.

Multi Chain

csb.statistics.samplers.mc.multichain implements several extended-ensemble Monte Carlo sampling algorithms.

Here is a short example which shows how to sample from a PDF using the replica exchange with non-equilibrium switches (RENS) method. It draws 5000 samples from a 1D normal distribution using the RENS algorithm working on three Markov chains being generated by the HMC algorithm:


>>> import numpy
>>> from numpy import sqrt
>>> from csb.io.plots import Chart
>>> from csb.statistics.pdf import Normal
>>> from csb.statistics.samplers import State
>>> from csb.statistics.samplers.mc.multichain import ThermostattedMDRENSSwapParameterInfo
>>> from csb.statistics.samplers.mc.multichain import ThermostattedMDRENS, AlternatingAdjacentSwapScheme
>>> from csb.statistics.samplers.mc.singlechain import HMCSampler

Pick some initial state for the different Markov chains:

>>> initial_state = State(numpy.array([1.](1.)))

Set standard deviations:

>>> std_devs = [1./sqrt(5), 1. / sqrt(3), 1.](1._sqrt(5),-1.-_-sqrt(3),-1.)

Set HMC timesteps, trajectory length and gradients:

>>> hmc_timesteps = [0.6, 0.7, 0.7](0.6,-0.7,-0.7)
>>> hmc_trajectory_length = 20
>>> hmc_gradients = [lambda q, t: 1 / (std_dev ** 2) * q for std_dev in std_devs](lambda-q,-t_-1-_-(std_dev-__-2)-_-q-for-std_dev-in-std_devs)

Set parameters for the thermostatted RENS algorithm:

>>> rens_trajectory_length = 30
>>> rens_timesteps = [0.5, 0.3, 0.1](0.5,-0.3,-0.1)

Set interpolation gradients as a function of the work parameter l:

>>> rens_gradients = [lambda q, l, i=i: (l / (std_devs[i + 1](lambda-q,-l,-i=i_-(l-_-(std_devs[i-+-1) ** 2) + (1 - l) / (std_devs[i](i) ** 2)) * q 
                      for i in range(len(std_devs)-1)]

Initialize HMC samplers:

>>> samplers = [HMCSampler(Normal(sigma=std_devs[i](i)(i)(HMCSampler(Normal(sigma=std_devs[i)), initial_state, hmc_gradients[i](i)(i), hmc_timesteps[i](i)(i),
                hmc_trajectory_length) for i in range(len(std_devs))]

Create swap parameter objects:

params = [ThermostattedMDRENSSwapParameterInfo(samplers[0](0)(ThermostattedMDRENSSwapParameterInfo(samplers[0), samplers[1](1), rens_timesteps[0](0),
          rens_trajectory_length, rens_gradients[0](0)),
          ThermostattedMDRENSSwapParameterInfo(samplers[1](1)(1), samplers[2](2), rens_timesteps[1](1)(1),
          rens_trajectory_length, rens_gradients[1](1))]

Initialize thermostatted RENS algorithm:

>>> algorithm = ThermostattedMDRENS(samplers, params)

Initialize swapping scheme:

>>> swapper = AlternatingAdjacentSwapScheme(algorithm)

Initialize empty list which will store the samples and obtain samples:

>>> states = []()
>>> for i in range(5000):
        if i % 5 == 0:
            swapper.swap_all()
        states.append(algorithm.sample())

Print acceptance rates:

>>> print('HMC acceptance rates:', [s.acceptance_rate for s in samplers](s.acceptance_rate-for-s-in-samplers))
>>> print('swap acceptance rates:', algorithm.acceptance_rates)

Create and plot histogram for first sampler and numpy.random.normal reference:

>>> chart = Chart()
>>> rawstates = [state[0](0)(state[0).position[0](0) for state in states]
>>> chart.plot.hist([numpy.random.normal(size=5000), rawstates](numpy.random.normal(size=5000),-rawstates), bins=30, normed=True)
>>> chart.plot.legend(['numpy.random.normal', 'RENS + HMC']('numpy.random.normal',-'RENS-+-HMC'))
>>> chart.show()

For ReplicaExchangeMC (RE), the procedure is easier because apart from the two sampler instances the corresponding RESwapParameterInfo objects take no arguments.

Every replica exchange algorithm in this module (ReplicaExchangeMC, MDRENS, ThermostattedMDRENS) is used in a similar way. A simulation is always initialized with a list of samplers (instances of classes derived from AbstractSingleChainMC) and a list of AbstractSwapParameterInfo objects suited for the algorithm under consideration. Every AbstractSwapParameterInfo object holds all the information needed to perform a swap between two samplers. The usual scheme is to swap only adjacent replicae in a scheme:

   1 <--> 2, 3 <--> 4, ...
   2 <--> 3, 4 <--> 5, ...
   1 <--> 2, 3 <--> 4, ...

This swapping scheme is implemented in the AlternatingAdjacentSwapScheme class, but different schemes can be easily implemented by deriving from AbstractSwapScheme. Then the simulation is run by looping over the number of samples to be drawn and calling the AbstractExchangeMC.sample method of the algorithm. By calling the AbstractSwapScheme.swap_all method of the specific AbstractSwapScheme implementation, all swaps defined in the list of AbstractSwapParameterInfo objects are performed according to the swapping scheme. The AbstractSwapScheme.swap_all method may be called for example after sampling intervals of a fixed length or randomly.