PosteriorSampler

---

In this tutorial, we will perform sampling when given an ensemble. Previously, in the Ensemble tutorial we constructed ensembles for two lambda values and saved them to pickle files. Please read the section lambda values in the previous tutorial for more information.

---

import pickle, os
import numpy as np
import biceps

outdir = 'results'
n_lambdas = 2
lambda_values = np.linspace(0.0, 1.0, n_lambdas)
ensembles = []
for lam in lambda_values:
    with open(outdir+"/ensemble_%s.pkl"%lam,'rb') as file:
        ensemble = pickle.load(file)
    ensembles.append(ensemble)

Next, we need to specify number of steps for BICePs sampling. We recommend to run at least 1M steps for converged Monte Carlo samplings. Checking the convergence of the MCMC trajectory can be done simply using the submodule, biceps.Convergence.

Using the first ensemble of the list of ensembles…

nsteps = 100000 # number of steps of MCMC simulation
sampler = biceps.PosteriorSampler(ensembles[0])

Then, we sample for 100 k steps…

sampler.sample(nsteps=nsteps, burn=1000, print_freq=10000, verbose=True)

Step            State   Para Indices            Avg Energy      Acceptance (%)
0               [87]    [196, 127, 195]         7.150           56.94   False
10000           [67]    [235, 150, 184]         8.713           69.42   False
20000           [67]    [267, 141, 202]         8.858           73.08   True
30000           [68]    [243, 137, 196]         9.303           72.81   True
40000           [27]    [229, 150, 199]         7.451           72.53   True
50000           [25]    [230, 134, 197]         9.089           72.14   True
60000           [65]    [276, 129, 182]         9.437           72.21   False
70000           [18]    [224, 147, 198]         11.604          72.58   True
80000           [21]    [302, 117, 182]         11.992          73.06   True
90000           [79]    [241, 135, 199]         7.123           73.16   True

Accepted 72.88811881188118 %


Accepted [ ...Nuisance paramters..., state] %
Accepted [32.68811881 31.30891089 31.30891089  8.89108911] %

Now, process and save the results for the next step, Analysis.

lam = lambda_values[0]
filename = os.path.join(outdir,'traj_lambda%2.2f.npz'%(lam))
sampler.traj.process_results(filename)

Then, we will do the same thing for the second ensemble in the list of ensembles. Note that we separated them for simplicity of the tutorial.

lam = lambda_values[1]
sampler = biceps.PosteriorSampler(ensembles[1])
sampler.sample(nsteps=nsteps, burn=1000, verbose=False)
filename = os.path.join(outdir,'traj_lambda%2.2f.npz'%(lam))
sampler.traj.process_results(filename)
biceps.toolbox.save_object(sampler, filename.replace(".npz", ".pkl"))

 99%|██████████████████████████████████████████████████████████▍| 100000/101000 [00:01<00:00, 61468.68it/s]

Accepted 65.55940594059406 %


Accepted [ ...Nuisance paramters..., state] %
Accepted [32.73762376 31.24950495 31.24950495  1.57227723] %

Conclusion

In this tutorial, we used the biceps.PosteriorSampler class to perform MCMC sampling given the ensemble Python object. In the next tutorial, we will analyze our trajectory data.

# NOTE: The following cell is for pretty notebook rendering

from IPython.core.display import HTML
def css_styling():
    styles = open("../../../theme.css", "r").read()
    return HTML(styles)
css_styling()