7' Metropolis and Other Algorithms

1
7. Metropolis and Other Algorithms
2
Markov Chain and Monte Carlo

• Markov chain theory describes a particularly
  simple type of stochastic processes. Given a
  transition matrix, W, the invariance distribution
  P can be determined.
• Monte Carlo is a computer implementation of a
  Markov chain. In Monte Carlo, P is given, we
  need to find W such that P P W.

3
A Computer Simulation of a Markov Chain

• Let ?0, ?1, , ?n , be a sequence of
  independent, uniformly distributed random
  variable between 0 and 1
• Partition the set 0,1 into subintervals Aj such
  that Aj(p0)j, and Aij such that AijW(i
  -gtj), and define two functions
• G0(u) j if u ?Aj
• G(i,u) j if u ?Aij, then the chain is
  realized
• 4. X0 G0(?0), Xn1 G(Xn,?n1) for n 0

4
Partition of Intervals
u
A1
A2
A3
A4
1
0
p1
p1p2
p1p2p3
A1 is the subset 0,p1), A2 is the subset p1,
p1p2), , such that length Aj pj.
5
Markov Chain Monte Carlo

• Generate a sequence of states X0, X1, , Xn, such
  that the limiting distribution is given P(X)
• Move X by the transition probability
• W(X -gt X)
• Starting from arbitrary P0(X), we have
• Pn1(X) ?X Pn(X) W(X -gt X)
• Pn(X) approaches P(X) as n go to 8

6
Necessary and sufficient conditions for
convergence

• Ergodicity
• Wn(X - gt X) gt 0
• For all n gt nmax, all X and X
• Detailed Balance
• P(X) W(X -gt X) P(X) W(X -gt X)

7
Taking Statistics

• After equilibration, we estimate

It is necessary that we take data for each sample
or at uniform interval. It is an error to omit
samples (condition on things).
8
Choice of Transition Matrix W

• The choice of W determines a algorithm. The
  equation
• P PW or P(X)W(X-gtX)P(X)W(X-gtX)
• has (infinitely) many solutions given P.
• Any one of them can be used for Monte Carlo
  simulation.

9
Metropolis Algorithm (1953)

• Metropolis algorithm takes
• W(X-gtX) T(X-gtX) min(1, P(X)/P(X))
• where X ? X, and T is a symmetric stochastic
  matrix
• T(X -gt X) T(X -gt X)

10
The Paper (7500 citations from 1988 to 2003)
THE JOURNAL OF CHEMICAL PHYSICS VOLUME 21,
NUMBER 6 JUNE, 1953
Equation of State Calculations by Fast Computing
Machines
NICHOLAS METROPOLIS, ARIANNA W. ROSENBLUTH,
MARSHALL N. ROSENBLUTH, AND AUGUSTA H.
TELLER, Los Alamos Scientific Laboratory, Los
Alamos, New Mexico AND EDWARD TELLER,
Department of Physics, University of Chicago,
Chicago, Illinois (Received March 6, 1953) A
general method, suitable for fast computing
machines, for investigating such properties as
equations of state for substances consisting of
interacting individual molecules is described.
The method consists of a modified Monte Carlo
integration over configuration space. Results
for the two-dimensional rigid-sphere system have
been obtained on the Los Alamos MANIAC and are
presented here. These results are compared to
the free volume equation of state and to a
four-term virial coefficient expansion.
1087
11
(No Transcript)
12
Model Gas/Fluid
A collection of molecules interact through some
potential (hard core is treated), compute the
equation of state pressure p as function of
particle density ?N/V.
(Note the ideal gas law) PV N kBT
13
The Statistical Mechanics Problem

• Compute multi-dimensional integral
• where potential energy

14
Importance Sampling

• , instead of choosing configurations randomly,
  , we choose configuration with a probability
  exp(-E/kBT) and weight them evenly.
• - from M(RT)2 paper

15
The M(RT)2

• Move a particle at (x,y) according to
• x -gt x (2?1-1)a, y -gt y (2?2-1)a
• Compute ?E Enew Eold
• If ?E 0 accept the move
• If ?E gt 0, accept the move with probability
  exp(-?E/(kBT)), i.e., accept if
• ?3 lt exp(-?E/(kBT))
• Count the configuration as a sample whether
  accepted or rejected.

16
A Priori Probability T

• What is T(X-gtX) in the Metropolis algorithm?
  And Why it is symmetric?

17
The Calculation

• Number of particles N 224
• Monte Carlo sweep 60
• Each sweep took 3 minutes on MANIAC
• Each data point took 5 hours

18
MANIAC the Computer and the Man
Seated is Nick Metropolis, the background is the
MANIAC vacuum tube computer
19
Summary of Metropolis Algorithm

• Make a local move proposal according to T(Xn -gt
  X), Xn is the current state
• Compute the acceptance rate
• r min1, P(X)/P(Xn)
• Set

20
Metropolis-Hastings Algorithm

• where X?X. In this algorithm we remove the
  condition that T(X-gtX) T(X-gtX)

21
Why Work?

• We check that P(X) is invariant with respect to
  the transition matrix W. This is easy if
  detailed balanced is true. Take
• P(X) W(X -gt Y) P(Y) W(Y-gtX)
• Sum over X, we get
• ?XP(X)W(X-gtY) P(Y) ?XW(Y-X) P(Y)

?X W(Y-gtX) 1, because W is a stochastic matrix.
22
Detailed Balance Satisfied

• For X?Y, we have W(X-gtY) T(X-gtY) min1,
  P(Y)T(Y-gtX)/(P(X)T(X-gtY))
• So if P(X)W(X-gtY) gt P(Y)T(Y-gtX), we get
  P(X)W(X-gtY) P(X) T(X-gtY), and P(Y)W(Y-gtX)
  P(Y) T(Y-gtX)
• Same is true then inequality is reversed.

Detailed balance means P(X)W(X-gtY) P(Y) W(Y-gtX)
23
Ergodicity

• The unspecified part of Metropolis algorithm is
  T(X-gtX), the choice of which determines if the
  Markov chain is ergodic.
• Choice of T(X-gtX) is problem specific. We can
  adjust T(X-gtX) such that acceptance rate r 0.5

24
Gibbs Sampler or Heat-Bath Algorithm

• If X is a collection of components, X(x1,x2, ,
  xi, , xN),
• and if we can compute
• P(xix1,,xi-1,xi1,..,xN),
• we generate the new configuration by sampling xi
  according to the above conditional probability.