## New Research In

### Physical Sciences

### Social Sciences

#### Featured Portals

#### Articles by Topic

### Biological Sciences

#### Featured Portals

#### Articles by Topic

- Agricultural Sciences
- Anthropology
- Applied Biological Sciences
- Biochemistry
- Biophysics and Computational Biology
- Cell Biology
- Developmental Biology
- Ecology
- Environmental Sciences
- Evolution
- Genetics
- Immunology and Inflammation
- Medical Sciences
- Microbiology
- Neuroscience
- Pharmacology
- Physiology
- Plant Biology
- Population Biology
- Psychological and Cognitive Sciences
- Sustainability Science
- Systems Biology

# Monte Carlo sampling for stochastic weight functions

Contributed by Daan Frenkel, May 18, 2017 (sent for review December 13, 2016; reviewed by David Ceperley and Athanassios Z. Panagiotopoulos)

## Significance

Markov chain Monte Carlo is the method of choice for sampling high-dimensional (parameter) spaces. The method requires knowledge of the weight function (or likelihood function) determining the probability that a state is observed. However, in many numerical applications the weight function itself is fluctuating. Here, we present an approach capable of tackling this class of problems by rigorously sampling states proportionally to the average value of their fluctuating likelihood. We demonstrate that the method is capable of computing the volume of a basin of attraction defined by stochastic dynamics as well as being an efficient method to identify a transition state along a known reaction coordinate. We briefly discuss how the method might be extended to experimental settings.

## Abstract

Conventional Monte Carlo simulations are stochastic in the sense that the acceptance of a trial move is decided by comparing a computed acceptance probability with a random number, uniformly distributed between 0 and 1. Here, we consider the case that the weight determining the acceptance probability itself is fluctuating. This situation is common in many numerical studies. We show that it is possible to construct a rigorous Monte Carlo algorithm that visits points in state space with a probability proportional to their average weight. The same approach may have applications for certain classes of high-throughput experiments and the analysis of noisy datasets.

Monte Carlo simulations aim to sample the states of the system under study such that the frequency with which a given state is visited is proportional to the weight (often “Boltzmann” weight) of that state. The equilibrium distribution of a system, that is, the distribution for which every state occurs with a probability proportional to its (Boltzmann) weight, is invariant under application of a single Monte Carlo step. Algorithms that satisfy this criterion are said to satisfy “balance” (1). Usually, we impose a stronger condition, “detailed balance,” which implies that the average rate at which the system makes a transition from an arbitrary “old” state (

Many simple Monte Carlo (MC) algorithms satisfy in addition microscopic reversibility, which means that

### Monte Carlo Simulations with Noisy Acceptance Rules.

There are many situations where conventional MCMC cannot be used because the quantity that determines the weight of a state

Note that the problem that we are discussing here is different from the cases considered by Bhanot and Kennedy (4) and by Ceperley and Dewing (5). As we discuss below, these earlier papers consider cases where the weights are nonlinear functions of a fluctuating argument (e.g., an action or an energy), in which case the average of the function is not equal to the function of the average argument. In contrast, we consider the case where the probability to sample a point is given rigorously by the average of the stochastic estimator of the weight function.

To give a specific example, we consider the problem of computing the volume of the basin of attraction of a particular energy minimum

However, many minimizers are not deterministic—and hence the oracle function is probabilistic. (In fact, historical evidence suggests that ancient oracles were probabilistic at best). In that case, if we start a number of minimizations at point

We could obtain an estimate for the average weight

First, however, we generalize the concept of a basin volume **5**, we must be able to sample points with a probability

### Naive MC Algorithm.

We first describe a naive (and very inefficient) but rigorous MC algorithm to sample stochastic weight functions. After that, we show how the algorithm can be made more efficient.

Our aim is to construct a MC algorithm that will visit points **5** becomes a free-energy calculation, for which standard techniques exist (2).

Let us consider two points (

Note that in this naive version of the algorithm, the acceptance rule is not the Metropolis rule that considers the ratio of two weights. Here it is the probability itself. Hence, whenever the probability becomes very low, the acceptance of moves decreases proportionally. We address this problem in what follows.

### “Configurational-Bias” Approach.

With the naive algorithm described above, the acceptance of moves becomes small when the system moves into a region of configuration space where

We can then generate a random walk in this extended space, between points that are surrounded by a “cloud” of

The coordinates of the

We now define an extended state space

We can now construct a MCMC to visit (but not sample) backbone points. To this end, we compute the “Rosenbluth weight” of point

We can then construct a MCMC algorithm where the acceptance of a trial move from the old

Note that during a trial move, the state of the old point is not changed, and hence it retains the same trial directions (hence the same set

The approach of Eq. **13** can be easily incorporated into more sophisticated sampling schemes such as parallel tempering (PT) (12, 13), as discussed in *Supporting Information* and shown in Fig. 2.

### Sampling.

We have shown that backbone points will be visited with a probability proportional to its instantaneous Rosenbluth weight **13** for the acceptance of moves between those backbone points. In contrast, Eq. **14** expresses the detailed balance condition for transitions between cloud points. In what follows, we assume that the probability

To achieve the desired sampling of cloud points, we impose that a given cloud point

The approach that we describe here is better than the naive algorithm because it achieves faster “diffusion” through parts of configuration space where

However, even though Rosenbluth-style sampling ensures that all points in space are sampled with the correct frequency, it is not an efficient algorithm. The reason is obvious: To compute the weights

Fortunately, this drawback can be overcome. Rather than sampling one point at a time, we take steps between backbone points sampled according to Eq. **13** and compute the quantity to be sampled for all

For every backbone point

### Combine with “Waste-Recycling” MC.

Efficiency can be further improved by using the approach underlying “waste-recycling” MC (14). This approach allows us to sample all trial cloud points in the sampling, even if the actual trial backbone move is rejected. The approach of ref. 14 allows us to combine the information of the accepted and the rejected states in our sampling. Specifically, we denote the probability to accept a move from an old state

## Numerical Results

### Basin Volume Calculations.

To test the proposed algorithm we compute the basin volume (probability mass) for a stochastic oracle function as defined in Eq. **5**. We choose a few simple oracle functions, for which the integral in Eq. **5** can be solved analytically. The volume calculations were performed using the multistate-Bennett acceptance ratio (MBAR) method (16) described in ref. 9. As described in ref. 9, a high-dimensional volume calculation is in essence a free-energy calculation, where minus the log of the volume plays the role of the free energy.

We compute the dimensionless free-energy difference between a region of known volume

We first tested the method for a deterministic oracle, namely a simple **13** for

Next, we tested the method for a stochastic oracle function defined such that

### Transition-State Finding.

The algorithm that we described above has wider applicability than the specific examples that we discussed. As an illustration of a very different application, we show that our approach can be used to efficiently identify the transition state along a known reaction coordinate.

Note that points in the transition-state ensemble (in the one-dimensional case: just one point) are characterized by the property that the committor has an average value of 0.5. However, any individual trajectory will either be crossing (“1”) or noncrossing (“0”). Hence, the “signal” is stochastic. As an illustration, we consider the (trivial) one-dimensional case of a particle with kinetic energy

## Relation to Earlier Work

The problem of Monte Carlo sampling in the presence of noise has been discussed by Bhanot and Kennedy (4) and Ceperley and Dewing (5).

Bhanot and Kennedy (4) considered how to construct an unbiased estimator of an exponential function (e.g., a ratio of Boltzmann weights) with a fluctuating argument. This method involves constructing an estimator on the basis of a number of independent samples. The method is subject to certain limitations (it is not guaranteed to generate acceptance probabilities between

## Conclusions and Outlook

Thus far, the algorithm described above was presented as a method to perform MC sampling in cases where the weight function itself is fluctuating.

However, we suggest that the method is not limited to numerical sampling: It could be used to steer sampling of experimental control parameters in experiments that study stochastic events (e.g., crystal nucleation, cell death, or even the effect of advertising). Often, the occurrence of the desired event depends on a large number of variables (temperature, pressure, pH, concentration of various components) and we want to select the optimal combination. However, as the desired event itself is stochastic, individual measurements provide little guidance. One might aim to optimize the conditions by accumulating sufficient statistics for individual state points. However, such an approach is expensive. The procedure described in the preceding sections suggests that it may be better to perform experiments in a cloud of state points around a backbone point. We could then accept or reject the trial move to a new backbone state, using the same rule as in Eq. **13**. In this way, the experiment could be made to evolve toward “interesting” regions of parameter space.

## PT

PT (12, 13) is a MC scheme that targets the slow equilibration of systems characterized by large free-energy barriers that prevent the efficient equilibration of a MCMC random walk. In PT,

The configurational bias approach to cloud sampling embodied by Eq. **13** can be easily generalized to PT to find an acceptance rule for the swap of configurations between replicas

We have tested this method both for a deterministic oracle—a simple **20** as shown in Fig. S1 (compare with Fig. 3 of the main text).

## Acknowledgments

D.F. acknowledges support by Engineering and Physical Sciences Research Council (EPSRC) Program Grant EP/I001352/1 and EPSRC Grant EP/I000844/1. K.J.S. acknowledges support by the Swiss National Science Foundation (Grants P2EZP2-152188 and P300P2-161078). S.M. acknowledges financial support from the Gates Cambridge Scholarship.

## Footnotes

- ↵
^{1}To whom correspondence should be addressed. Email: df246{at}cam.ac.uk.

Author contributions: D.F. and S.M. designed research; D.F., K.J.S., and S.M. performed research; K.J.S. and S.M. analyzed data; and D.F. and S.M. wrote the paper.

Reviewers: D.C., University of Illinois at Urbana–Champaign; and A.Z.P., Princeton University.

The authors declare no conflict of interest.

This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.1073/pnas.1620497114/-/DCSupplemental.

## References

- ↵
- ↵.
- Frenkel D,
- Smit B

- ↵
- ↵.
- Bhanot G,
- Kennedy AD

- ↵
- ↵
- ↵.
- Asenjo D,
- Paillusson F,
- Frenkel D

- ↵.
- Martiniani S,
- Schrenk KJ,
- Stevenson JD,
- Wales DJ,
- Frenkel D

- ↵.
- Martiniani S,
- Schrenk KJ,
- Stevenson JD,
- Wales DJ,
- Frenkel D

- ↵.
- Martiniani S,
- Schrenk KJ,
- Ramola K,
- Chakraborty B,
- Frenkel D

*Nat Phys*, in press. - ↵.
- Frenkel D,
- Mooij GCAM,
- Smit B

- ↵
- ↵
- ↵.
- Frenkel D

- ↵
- ↵
- ↵.
- Yan Q,
- de Pablo JJ

- ↵.
- Bunker A,
- Dünweg B

- ↵

## Citation Manager Formats

### More Articles of This Classification

### Physical Sciences

### Related Content

- No related articles found.

### Cited by...

- No citing articles found.