Page 366 - Numerical Methods for Chemical Engineering

Markov chains and processes; Monte Carlo methods 355

The acceptance probability is related to the generating probability distribution P(q)by
the condition of detailed balance,
[k] [k] (new) [k] (new)
P q γ q → q α q → q
(new) (new) [k] (new) [k]
= P q γ q → q α q → q (7.233)
To understand the origin and meaning of this relation, consider the following thought
experiment. Let us say that we were generating not one sequence, but many sequences
independently and in parallel. At each step, the number of sequences in our ensemble at or
near q should be proportional to P(q). Therefore, from one step to the next in simulating
our ensemble of sequences, the number of sequences that move away from any point q
to any other point should be balanced by the number moving to q from all other points.
The condition of detailed balance goes yet further, and states that the number of sequences
making the move q [k] → q (new) must be balanced by the number making the move q (new) →
[k]
q . This is represented mathematically as the flux balance (7.233). From this condition,
we obtain the following rule for the ratio of the acceptance probabilities,
[k] (new) (new) (new) [k]
α q → q P q γ q → q
= (7.234)
(new) [k] [k] [k] (new)
α q → q P q γ q → q
We often pick a generating process that is symmetric,
[k] (new) (new) [k]
γ q → q = γ q → q (7.235)
For example, we can displace each component at random,
q m (new) ← q m [k] + m (u m − 0.5) (7.236)

where u m is a random variable (independent for each component) that is uniformly dis-
tributed between zero and one. m is some maximum allowable displacement that may be
tuned dynamically to optimize the fraction of moves that are accepted to the range 0.1–0.5.
With a symmetric process for generating the moves, the ratio of acceptance probabilities
becomes
[k] (new) (new)
α q → q P q
= (7.237)
(new) [k] [k]
α q → q P q
A choice that satisfies this condition is
[k] (new) (new) [k]
α q → q = min 1, P q /P q (7.238)
which yields the Metropolis Monte Carlo method, based upon iterating the following
procedure:

generate a proposed move q [k] →q (new) at random from γ (q [k] →q (new) )
[k]
compute P(q (new) )/P(q )
generate a random variable u, uniformly distributed on [0, 1]
[k]
if u ≤ min{1, P(q (new) )/P(q )}, q [k+1] = q (new) ; else, q [k+1] = q [k]
To simulate a physical system using this method, we start at some initial state, and perform
some large number of Monte Carlo steps to equilibrate the system. Initially, we may start at