Biometrika (1970), 57, 1, p. 97
Printed in Great Britain

97

Monte Carlo sampling methods using Markov
chains and their applications
B Y W. K. HASTINGS
University of Toronto
SUMMARY

A generalization of the sampling method introduced by Metropolis et al. (1953) is presented along with an exposition of the relevant theory, techniques of application andmethods and difficulties of assessing the error in Monte Carlo estimates. Examples of themethods, including the generation of random orthogonal matrices and potential applications of the methods to numerical problems arising in statistics, are discussed.
1. INTRODUCTION

For numerical problems in a large number of dimensions, Monte Carlo methods are often
more efficient than conventional numerical methods. However, implementation of theMonte Carlo methods requires sampling from high dimensional probability distributions
and this may be very difficult and expensive in analysis and computer time. General methods
for sampling from, or estimating expectations with respect to, such distributions are as
follows.
(i) If possible, factorize the distribution into the product of one-dimensional conditional
distributions from which samples may be obtained.
(ii) Use importance sampling, which may also be used for variance reduction. That is, inorder to evaluate the integralp

J = jf(x)p(x)dz = Ep(f),where p(x) is a probability density function, instead of obtaining independent samples
xv ...,xNfromp(x) and using the estimate Jx = ^/(x^/N, we instead obtain the sample froma distribution with density q(x) and use the estimate <72 = 'L{f(xi)p(xi)}l{q(xi)N}. This may
be advantageous if it is easier to sample from ^(a;) thanp(x), but it is a difficult method to use
in a large number of dimensions, since the values of the weights w{xi) = p(xi)lq{xi) for
reasonable values of N may all be extremely small, or a few may be extremely large. In
estimating the probability of an event A, however, these difficulties may not be as serious
since the only values of w(z) which are important are those for which x e A. Since themethods proposed by Trotter & Tukey (1956) for the estimation of conditional expectations
require the use of importance sampling, the same difficulties may be encountered in their use.
(iii) Use a simulation technique; that is, if it is difficult to sample directly from p(x) orif p(x) is unknown, sample from some distribution q(y) and obtain the sample x values assome function of the corresponding y values. If we want samples from the conditional distribution of x = g{y), given h(y) = h0, then the simulation technique will not be satisfactory
if pr {h(y) = h0) is small, since the condition h(y) = h0 will be rarely if ever satisfied evenwhen the sample size from q(y) is large.
7-2

98

W. K. HASTINGSIn this paper, we shall consider Markov chain methods of sampling that are generalizations of a method proposed by Metropolis et al. (1953), which has been used extensively fornumerical problems in statistical mechanics. An introduction to Metropolis's method and its
applications in statistical mechanics is given by Hammersley & Handscomb (1964, p. 117).
The main features of these methods for sampling from a distribution with density p(x) are:(a) The computations depend on p(x) only through ratios of the form p(x')jp(x), wherez' and x are sample points. Thus, the normalizing constant need not be known, no factorization of p(x) is necessary, and the methods are very easily implemented on a computer. Also,
conditional distributions do not require special treatment and therefore the methods
provide a convenient means for obtaining correlated samples from, for example, the conditional distributions given by Fraser (1968) or the distributions of the elements in a multiway
table given the marginal totals.
(b) A sequence of samples is obtained by simulating a Markov chain. The resulting
samples are therefore correlated and estimation of the standard deviation of an estimate
and assessment of the error of an estimate may require more care than with independent
samples.
2. DEPENDENT SAMPLES USING MARKOV CHAINS

2-1. Basic formulation of the method
Let P = {py} be the transition matrix of an irreducible Markov chain with states 0,1,..., S.
Then, if X(t) denotes the state occupied by the process at time t, we have

If n = (n0, TTlt..., ns) is a probability distribution with ni > 0 for all i, and if/(•) is a function
denned on the states, and we wish to estimate

we may do this in the following way. Choose P so that n is its unique stationary distribution,i.e. 7t = 7iP. Simulate this Markov chain for times t = 1, ...,N and use the estimate

For finite irreducible Markov chains we know that / i s asymptotically normally distributed
and that / - » / in mean square as N->ao (Chung, 1960, p. 99).
In order to estimate the variance of / , we observe that the process X(t) is asymptoticallystationary and hence so is the process Y(t) = f{X(t)}. The asymptotic variance of the meanof such a process is independent of the initial distribution of X(0), which may, for example,
attach probability 1 to a single state, or may be 7t itself, in which case the process isstationary. Thus, if N is large enough, we may estimate var (/), using results appropriate
for estimating the variance of the mean of a stationary process.
Let pj be the correlation of Y(t) and Y(t+j) and let cr2 = v&v{Y(t)}. It is well known(Bartlett, 1966, p. 284) that for a stationary process
r2

2V-1

Monte Carlo methods using Markov chains
and that as iV -+ oo

99

var (F) ~ 2ng{0)IN,

where g(o)) is the spectral density function at frequency w. If the pi are negligible for j ^ j 0 ,then we may use Hannan's (1957) modification of an estimate of var(F) proposed by
Jowett (1955), namely

where

N-J
c,= 2T(t)Y(t+j)l(N-j)

for j > 0 a n d c_j = cy

A satisfactory alternative which is less expensive to compute is obtained by making useof the pilot estimate, corrected for the mean, for the spectral density function at zerofrequency suggested by Blackman & Tukey (1958, p. 136) and Blackman (1965). We divide
our observations into L groups of K consecutive observations each. Denoting the mean of
the ith block byK
we use the estimate

L _

_

4

i)}.

(2)

This estimate has approximately the stability of a chi-squared distribution on (L— 1)
degrees of freedom. Similarly, the covariance of the means of two jointly stationary processes
Y(t) and Z(t) may be estimated bys?z = 2 &- Y) (Zt-Z)I{L(L-1)}.

(3)

il

2-2. Construction of the transition matrix
In order to use this method for a given distribution n, we must construct a Markov chain P
with 7i as its stationary distribution. We now describe a general procedure for doing thiswhich contains as special cases the methods which have been used for problems in statisticalmechanics, in those cases where the matrix P was made to satisfy the reversibility condition
that for all i and j
The property ensures that Yi/niipii = np for all,?', and hence that n is a stationary distributionof P. The irreducibility of P must be checked in each specific application. It is only necessary
to check that there is a positive probability of going from state i to state j in some finite
number of transitions, for all pairs of states i and j .
We assume that p^ has the form
Pii = <lija-ij ( * * J ) ,

(5)

with
Pa = 1 - 2

Pa,

where 0 = fey} is the transition matrix of an arbitrary Markov chain on the states
0,1,..., S and ai:j is given by
n.

100

W. K. HASTINGS

where s^ is a symmetric function of i smdj chosen so that 0 < a,^ ^ 1 for all i andj. With thisform for p^ it is readily verified that Tr^py = Ti^p^, as required. In order to simulate thisprocess we carry out the following steps for each time t:(i) assume that X(t) = i and select a state j using the distribution given by the ith
row of 0 ;(ii) take X(t + 1) = j with probability a{j and X(t + 1) = i with probability 1 — ati.
For the choices of s^ we will consider, only the quantity (7fy<fy*)/(n<?y) enters into thesimulation and we will henceforth refer to it as the test ratio.
Two simple choices for sti are given for all i and j byila

4f = i.
I}

With qtj = qH and slt = s$ we have the method devised by Metropolis et al. (1953) andwith <?y = qit and sfj- = s^' we have Barker's (1965) method.
Little is known about the relative merits of these two choices for s^, but when q^ = q^,
we haveThus we see that if n^ = nt, we will take X(t +1) =j with probability 1 with Metropolis's
method and with probability £ with Barker's method. This suggests that Metropolis's
method may be preferable since it seems to encourage a better sampling of the states.
More generally, we may choose

where the function g(z) is chosen so that 0 < g(x) < 1 + x for 0 < x ^ 1, and g(x) may itself
be symmetric in i and j . For example, we may choose g(x) = 1 + 2(^x)v with the constanty ^ l , obtaining effi* with y = 1 and sff with y = oo.
We may define a rejection rate as the proportion of times t for which X(t+ 1) = X(t).
Clearly, in choosing 0 , high rejection rates are to be avoided. For example, if X(t) = i and
i is near the mode of an unimodal distribution, then 0 should be chosen so that j is not too
far from i, otherwise, ^ / ^ will be small and it is likely that X(t + 1) = i. For each simulation
it is useful to record the rejection rate since a high rejection rate may be indicative of a poor
choice of initial state or transition matrix, or of a 'bug' in the computer program.
We shall apply these methods to distributions n defined mathematically on an infinite
sample space although, when we actually simulate the Markov chains on a digital computer,
we will have a large but finite number of states. When n is continuous, we will have a
discrete approximation to n. Let n(x) dfi(x), p(x, x') d/i(x') and q(x, x')d{t(x') be the probability elements for the distribution n and the Markov processes analogous to P and 0 ,respectively. Let the possible values for x and x' on the computer be x0,..., xs. These values
depend on the word length and on the representation in the computer, float-point, fixed-Monte Carlo methods using Markov chains

101

point, etc. The probability elements may be approximated by nfa) 8/iix^, p{xitand q{xit Xj) S/i{Xj), respectively. Substituting 77p(xi) d/ifa) for •ni etc., we haveand (7^gfy)/(7fy<fyi) is replaced by {n{xi)q{xi,xi))l{n(xi)q(xpXi)']. Therefore, so long as ccti ischosen so that it depends on n and 0 only through the quantity (nri<2fij)/(77V9Vi)> we may use
the densities in place of the corresponding probability mass functions, and we may thinkin terms of continuous distributions ignoring the underlying approximation which will be
more than adequate in most applications.
For simplicity we have considered reversible Markov chains only. An example of an
irreversible P is given in Handscomb (1962), where the states may be subdivided into
finite subsets of states, where the states within a given subset are equally probable. Transitions amongst states within such a subset are made in a cyclic fashion; all other transitions
are reversible.
2-3. Elementary examples
Example 1. Let n be the Poisson distribution with 7ri = A* e~x\i! (i = 0,1,...). For A small
we may use the following choice of 0 to generate samples from n:

Note that ni+1lni = Xj{i + 1) and ~ni_x\Tii = i/A so that the computations may be performed
rapidly. For A large 0 must be chosen so that step sizes greater than unity are permitted
or else very long realizations will be required in order to ensure adequate coverage of thesample space.
Example 2. To sample from the standard normal distribution we define Qk(k = 1,2) inthe following ways. Let X(t) be the state at time t and choose X'(t) to be uniformly distributed on the interval [ekX(t) — A, ekX(t) + &], where A > 0 is a constant, e1 = + 1 ande2 = — 1. To use Metropolis's method with either of these choices of 0 , both of which aresymmetric, at each time t compute the test ratio /?(<) = exp (%[X2(t) — {X'(t)}2]). If fi(t) ^ 1,
set X(t+1) = X'(t); if /?(«) < 1, X(t+1) = X'(t) with probability 0{t) and X(t+1) = X(t)with probability 1—/?(<). Computing time can be saved if we compare X2(t) with{X'(tj}2instead of fi(t) with 1 and compute the exponential only when {X'(t)}2 > X2(t).
Simulations carried out on an IBM 7094II using the above transition matrices yielded
the following estimates X of the mean of the distribution where, in all cases, we chose
N = 1000 and L = 25: with X(0) = 0 and A = 1, we obtained X = - 0-12 and s^ = 0-11
with ek = + 1, and X = — 0-013 and s^ = 0-02 with ek = — 1. The estimated standard deviation of the estimate is smaller with ek = — 1, and is comparable to the theoretical standard
deviation of the mean with 1000 independent observations, 0-031. Here we have an ideal
situation with a symmetric distribution and known mean, but a similar device for variancereduction may sometimes be applicable in other problems of more practical interest. Other
results indicated that there is little to choose between moderate values of A in the range
0-2-1-8, but extreme values of A led to poor estimates. For example, with X(0) = 1-0 and
A = 0-001 we obtained X = 1-001 and s 5 = 0-002.
2-4. Multidimensional distributions
If the distribution n is d-dimensional and the simulated process is X(i) = {X^t),...,
there are many additional techniques which may be used to construct P:

Xd(t)},

102

W. K. HASTINGS(1) In the transition from time t to (t+ 1) all co-ordinates of X(t) may be changed.
(2) In the transition from time t to t + 1 only one co-ordinate of X(<) may be changed,that selection being made at random from amongst the d co-ordinates.
(3) Only one co-ordinate may change in each transition, the co-ordinates being selected
in a fixed rather than a random sequence.
Method (3) was used by Ehrman, Fosdick & Handscomb (1960), who justified the method
for their particular application. A general justification for the method may be obtained in
the following way. Let the transition matrix when co-ordinate k is to be moved by P ; .
(k = 1, ...,d). Assume that the co-ordinates are moved in the order 1, 2,... and that theprocess is observed only at times 0,d,.... The resulting process is a Markov process withtransition matrix P = P x . . . Pd. If, for each k, Pk is constructed so that nPk = n, then nwill be a stationary distribution of P since 7iP = nP1 ...Pd = nP2...Pd=... = n. In
practice, we must check the irreducibility of P to ensure uniqueness of the stationary
distribution. Note that for the validity of this proof it is not necessary that each Pk satisfy
the reversibility conditions (4). Also, in our estimates we may average observed values ofthe function/(.) at all times t of the original process although it would not be desirable to
do so if the function values only change every d steps.
Example 3. To illustrate method (3) we consider sampling from a distribution withprobability density functionp(x) = p{xx, •••,xd) defined over the domain
0 < xr < ... ^ xd < oo.
Let the co-ordinates at time t be x(£) = {x1(<), ...,xd(t)} and assume co-ordinate k (k 4= d) is
to be changed. Let x'k{t) be chosen from the uniform distribution on the interval{xft-i(0> xfc+i(O}> where we assume that xo(t) = 0, and define x'^t) = x^t) (i 4= k). Using
Metropolis's method, we find that the test ratio is p{x'(t)}/p{x(t)}. When k = d, we may
choose x'd(t) from the uniform distribution on the interval
and use the test ratio[p{x'(t)} {xd(t) - z ^A few computer runs using this form of transition matrix have yielded satisfactory resultsfor (i) sampling from the distribution of eigenvalues of a Wishart matrix, and (ii) estimating
the conditional expected values of order statistics needed for the probability plotting
method proposed by Wilk & Gnanadesikan (1968). The latter problem requires estimationof Ep(xk) (k = l,...,d),where

p(x) = c II (A)*z^- 1 exp ( - Aasc)/r(A)(0 < xt < ... < xa),

the pi's are constants and c is a normalizing constant.
Generalizations of these methods are readily obtained by noting that, for example, theco-ordinates need not be moved equally often and that the co-ordinates may be moved ingroups instead of one at a time. For example, in statistical mechanics applications, one
particle and hence three co-ordinates are moved at a time.
When only a few of the co-ordinates are moved at one time, computing effort can usually
be saved by employing a recurrence formula for obtaining/{X(f + 1)} from/{X(f)}. However
it is then necessary to guard against excessive build-up of error. One way to accomplish
this is periodically to calculate/{X(t)} without the aid of the recursion, although an error

Monte Carlo methods using Markov chains

103

analysis in some cases may indicate that this is unnecssary. If, however, it is expensive tocompute either of nt or/(i), it may be preferable to attempt to move all of the co-ordinatesa small distance each at time t.
2-5. Importance sampling
Usually, P is constructed so that it depends only on the ratios Tti\ini as in the methods ofMetropolis and Barker. In this case we need only know the 77/s up to a constant of proportionality, and if n0 +... + ns 4= 1, we are estimatingf ( ) i l i
i=0

/i=0

This expression may be rewritten in the form
i

=

i=0

0

i=0

/

/

i=0

which we recognize as E^{f(i)7riln'^IEJl.(7ril7r'i). Hence, if we set up the Markov chain using
the distribution n' instead of 71 we can estimate / = En(f) by the ratio

This device permits us to use importance sampling with the Markov chain methods assuggested by Fosdick (1963).
If 7T0+... +ns = 1 and n'0 + ...+7T'S = 1, then we may replace the denominator in /by 1, simplifying the estimate. Otherwise we have a ratio estimate and we can estimate its
variance using the usual approximation. Thus, if we denote / by Y\Z, the variance isgiven approximately byvar (Y\Z) = {var (F) - 2/ cov (F, Z) + P var (Z)}l{E(Z)f.
This may be estimated by

4/

where «y, Spg and s~ are obtained as in (2) and (3). For a discussion of the validity of thisapproximation, see Hansen, Hurwitz & Madow (1953, p. 164).
Importance sampling for variance reduction is more easily implemented with the Markov
chain methods than with methods using independent samples, since with the Markov chain
methods it is not necessary to construct the distribution n' so that independent samples
can be obtained from it. If, however, we can obtain independent samples from the distribution n', we may obtain estimates ofEn(f) using the Markov chain methods, which have theadvantage that they do not involve the weights nx(t)ln'x(t)'> s e e the discussion of importancesampling in § 1. To accomplish this we set up the Markov chain P so that it has TT, not n',
as its stationary distribution and choose qi} = v\ for all i and j .

104

W. K. HASTINGS
2-6. Assessment of the error

In using Monte Carlo methods to estimate some quantity we usually use the standard
deviation of the estimate to obtain some indication of the magnitude of the error of theestimate. There are, of course, many sources of error common to all Monte Carlo methods
whose magnitude cannot be assessed using the standard deviation alone. These include:(i) the source of uniform random numbers used to generate the sample, which should beof as high a quality as possible; (ii) the nonnormality of the distribution of the estimate;(iii) computational errors which arise in computing the estimate; (iv) computation errors
which arise in generating the samples (including discretization and truncation of the distribution) ; and (v) errors induced because the sample size is too small, which are often best
overcome by methods other than increasing the sample size; for example, by the use ofimportance sampling. In what follows we shall concentrate upon categories (iv) and (v)above.
In generating successive samples using the Markov chain methods, errors will arise incomputing the new state X(t+ 1) and in computing the test ratio. An error analysis may
sometimes be useful for the computation of X(t +1) (see, for example, § 3), but it is difficult
to assess the effects of using inaccurate test ratios. The situation is also difficult to analyze,
in general, when successive samples are independent and are generated using a factorizationof the probability density function. To see this, letdp(x)=p(x1,x2,...,xd) = YiPMi)»=ibe the joint density function, where Pi(xt) is the conditional density function for x{ given
xx, ...,xi_1. When we attempt to sample from each of the one-dimensional distributions
Pi{Xi) in turn, errors will be introduced and we will instead be sampling from some distribution Pi{xt) = pt{xt) (l + Vi), where yt is a function of xv...,xi
and will not, in general, be
small for all sample points. Consequently, we will be generating samples from a distribution

and, especially when d is large, it will be difficult to assess the error in our estimate inducedby the yt'a. A similar, but more involved, analysis might also be applied to the Markov
chain methods if we consider a single realization of length N as a single sample from a
distribution of dimension Nd; sampling at each step of the Markov chain would correspondto sampling from a single factor above.
If the sample size is not large enough, important regions of the sample space may beinadequately represented. For example, if we are estimating the integral jf(x)p(x)dx bysampling from p(x) and if the major contribution to the value of the integral comes froma small region of low probability in which f(x) has very large values, then we may obtaina very poor estimate and a deceptively small standard deviation even with seemingly large
sample sizes. This difficulty may be encountered with any method of numerical quadratureand there is no substitute for a thorough study of the integrand and a consequent adjustmentof the method if gross errors are to be avoided. With the Markov chain methods the influenceon the result of the choice of initial state and the correlation of the samples, which may beconsiderable if the sample size is small, may be minimized by adopting the following
procedures:

Monte Carlo methods using Markov chains

105(i) Choose a transition matrix 0 so that the sample point in one step may move as largea distance as possible in the sample space, consistent with a low rejection rate. This offerssome protection against the possibility that the sample points for the whole realization
remain near one mode of the distribution which is separated from other modes by a deep
trough, or that the initial state is in a region of low probability.
(ii) If possible, choose the initial state so that it is in a region of high probability, by
sampling from n if this is possible, or set the co-ordinates of the initial state equal to their
expected values, or asymptotic approximations of these.
In view of the many sources of error one or more of the following techniques, depending
upon the application, should be used in practice to aid in assessing the suitability of theMarkov chain methods, the magnitude of the error and the adequacy of the length of therealization and of the standard deviation as a measure of error:(a) Test the method on problems which bear a close resemblance to the problem understudy and for which the results are known analytically.
(6) If the expected value of some function with respect to n is known, this may be estimated for checking purposes, and possibly for variance reduction, while other aspects of n
are under study.
(c) Compare estimates obtained from different segments of the same realization to see if
there is evidence of nonconvergence.
(d) Compare results obtained with and without the use of importance sampling, or using
different choices of P, or using different random numbers, or using the Markov chain method
and some other numerical method, for which adequate assessment of the error may also be
difficult as is often the case with asymptotic results, for example.
The illustrations given by Fox & Mayers (1968), for example, show how even the simplestof numerical methods may yield spurious results if insufficient care is taken in their use, andhow difficult it often is to assess the magnitude of the errors. The discussion above indicatesthat the situation is certainly no better for the Markov chain methods and that they should
be used with appropriate caution.
3. RANDOM ORTHOGONAL MATRICES

We now consider methods for generating random orthogonal and unitary matrices andtheir application to the evaluation of averages over the orthogonal group with respect toinvariant measure, a problem considered analytically by James (1955), and to other related
problems.
Let H be an orthogonal mxm matrix with |H| = 1 and let E^(0) be an elementary orthogonal matrix with elements given by
eu = cos 8, e^ = sin 6, e^ = — sin 6, e^ = cos 6,
eaa = 1(a #= i,j),eafi = 0 otherwise,for some angle 6. We may generate a sequence of orthogonal matrices H(J) (t = 1, 2,...)
using Metropolis's method, which is suitable for estimating averages over the orthogonalgroup with respect to invariant measure {dH}, in the following way:(i) Let H(0) = Ho, where H o is an orthogonal matrix.
(ii) For each t, select i and j at random from the set {1,..., m} with i 4= j . Select 0 from theuniform distribution on [0,2n].
(iii) Let W{t) = E

106

W. K. HASTINGS

Since the measure {dH} is left invariant, {dH'(t)} and {dH(tj} are equal and hence, inMetropolis's method, the new state is always the newly computed value, i.e. H(t +1) = H'(£).
To estimate the integral.
J=\/(H){dH},J0(m)where O(m) denotes the group of orthogonal matrices, we use the estimate

To remcve the restriction that |H| = 1 we may use the following procedure. After i and j
are selected as above, select one of i and j at random, each with probability \, and call the
result i'. Multiply the i'th row of E^(0) by + 1, the sign being selected at random, to forma new matrix which we denote by E^(0); E^(0) is then used in step (iii) above, in place ofE,,(0).
To show that the Markov chain is irreducible we need only show that every orthogonal
matrix H may be represented as the product of matrices of the form E^(0). This, in turn,
may be done by a simple annihilation argument. Choose 6 so that E12(0) H = H a ) has theelement h$ = 0. Choose 0 so that E13(0) H<» = H<2> has elements h$ = h$ = 0. Continue inthis way until we have the matrix H^m-i> with the first column annihilated. Since H (m ~ u is
orthogonal, we must have T^" 11 = ... = h^l~1) = 0 and h$-1) = + 1. Continuing this procedure we may annihilate the remaining off diagonal elements and obtain a diagonal
matrix with elements + 1 on the diagonal. The desired factorization is easily deduced
from this.
We now show that the process, by which the sequence of orthogonal matrices is generated,
is numerically stable. Denote the computed values of H(<) and E y (0) by Hc(£) and Ec(£). Let
H c (0 = H(«) + Ue(t), Ec(«) = E^(0) + Ee(t) and Hc(t + 1) = Ec{t) Hc(t) + F(t). We are interested
in how far He(*) is from zero, which we measure by ||He(<)(|, where the Euclidean norm || A||of a matrix A is defined by,,,

I|A|| = ( S S a% .
Using the fact that || A|| is preserved under orthogonal transformations and the inequalities||AB|| < ||A|| ||B|| and||A + B|| < ||A||+ ||B||, we can easily show that

where k is chosen such that ||EC(<)|| + ||F(f)|| ^ k for all t, and by induction that [|HC(<)|| < Ut,
where Ut is the solution of the difference equation Ut+1 = (l + k)Ut + k with Uo = ||He(0)||.
Solving the difference equation, we arrive at the boundTherefore, unless the word length of the computer is short or t is very large, the accumulation of error may be neglected since, for all t, \\ F(<)|| will be small in view of the few arithmetic
operations involved in forming the elements of the product Ec(<) Hc(t), and ||Ee(£)|| will also
be small if the method below is used for the generation of the four nontrivial elements ofEc(<). For short computer word length or very large t, orthonormalization of the columnsof Hc(<) at regular intervals will prevent excessive accumulation of error.

Monte Carlo methods using Markov chains107

One method for computing cos 8 and sin# required in Ei3-(0) is that of von Neumann.
Let TJX and U2 be independent and uniformly distributed on [0,1]. If XJ\ + V\ 4, 1, computecos0 = (tf|- I7|)/(J7|+ U\), sinS = ± 2Z71E/2/(£/f + U\),
where the sign of sin 6 is chosen at random. If U\ + TJ\ > 1, obtain new values of E/j and U2
until U\+U\^ 1.
Example 4. When/(H) = h\l + ...+ h2mm, then J = 1. Using N = 1000, L = 25, m = 50 and
H(0) = I, we obtained the estimate J = 3-5 with standard deviation 1-5. This estimate ispoor because of the very poor choice for H(0) but the estimated standard deviation gives us
adequate warning. It would be natural to increase the sample size N in order to obtaingreater precision but this may also be achieved by a better choice of H(0) with a saving incomputer time. For problems of the kind being considered here the following choice of H(0)is generally useful and better than that considered above. For j = 1,..., m and for m even, set
cos {(j-l)n),

hrj = V(2/m)cos{(j- 1) (r-2)27r/m)(r = 3,4, ...,\m+ 1),) s m {(j — l) (s — \m - 1) 2 W m } (s =For m odd a similar choice for H(0) may be made. Using this matrix for H(0), we obtained
the estimate J = 0-96 with standard deviation 0-03; the computer time required was
0-08 minutes.
The techniques discussed above may be applied to the following problems:(i) To generate random permutation matrices we set 0 = \n. To generate random
permutations, of m elements, at each time t, we interchange a randomly selected pair ofelements in the latest permutation in order to generate a new permutation.
(ii) To generate k orthogonal vectors in m dimensions we replace H(<) by an m x k matrix
where columns will be the desired vectors.
(iii) To generate random unitary matrices the above procedure for orthogonal matrices
is modified as follows:
In place of E^(0) we use an elementary unitary matrix U^(0, <f>) with elements given by
Uu = cos 6, Vti = e-^ sin 6, Ujt = - e^ sin 6, U^ = cos 6,
Uaa = 1(a + i, j),

Uafi = 0 otherwise.

The matrix U^(0,0) is obtained by multiplying the ith row of Uy(0, <f>) by eiy. Here, 6, <f>and y are chosen to be independent and uniformly distributed on [0,2TT]. Irreducibility may
again be established by an annihilation argument and the proof of numerical stability given
above requires only minor modifications.
(iv) To sample from a distribution defined on a group O = {g}, the transition matrix 0
must generate g'(t) from g(t), and this may be done by arranging that g'(t) = hg{t), where
heG and h is chosen from some distribution on G. If the probability element for g isp(g) {d/i(g)}, where n is the left invariant measure on G, then the ratio corresponding to
77^/77i is p{g'(t)}lp{g(t)}. As above we must ensure that g(t) is close to being a group element.
(v) To generate samples from a distribution with probability element
(Xm)

108

W. K. HASTINGS

we may remove the integral and sample from q(x, H) {dH} dx with, unfortunately, a
consequent increase in the dimension of the space being sampled. The sampling may be
carried out by combining the methods of § 2-4, alternating the changes of x and H, with thetechniques given above for orthogonal matrices. Note that the value of the test ratio forchanging H will no longer be unity. More generally, sampling from any distribution whose
density may be expressed as an average may be accomplished in a similar way. Many of thenormal theory distributions of multivariate analysis given by James (1964) have this form.
In the method described above for orthogonal matrices the successive matrices are
statistically dependent and we now consider a method in which successive matrices are
independent. Let x{j (i = l,...,k;j = l,...,m) be independent standard normal variates
and form the k vectors
xi = (xil,...,xim)(i=l,2,...,k).

If we now apply a Gram-Schmidt orthonormalization to these vectors, we will obtaink vectors y^ with the desired properties, since it is easy to see that the joint distribution of
the y i 's will be unchanged if they are subjected to an arbitrary orthogonal transformation.
For the case k = 1, see Tocher (1963). For each set of k orthogonal vectors generated by thisprocedure we require mk standard normal deviates, k square roots and, approximately,2mk2 arithmetic operations. The Markov chain method, on the other hand, requires only
6k arithmetic operations and the cosine and sine of a uniformly distributed angle and, if
the function/(•) being averaged is evaluated only at times T, 2T,..., we must have T as
large as \mk before the computing times of the two procedures are comparable. The
degree of correlation amongst the sample values for the Markov chain method will depend
upon the particular function/(•) and will determine, ignoring error analysis considerations,which of the two methods is better.
It is a pleasure to thank Messrs Damon Card, Stuart Whittington and Professor John P.
Valleau, who generously shared their considerable knowledge of and experience with theMarkov chain methods in statistical mechanics applications, Professor John C. Ogilvie, who
made several helpful suggestions and Mr Ross D.MacBride, who capably prepared thecomputer programs. I am grateful to the referees whose comments were very helpful in therevision of this paper. The work was supported by the National Research Council of Canada.

REFEKENCESBARKER, A. A. (1965). Monte Carlo calculations of the radial distribution functions for a proton-

electron plasma. Aust. J. Phys. 18, 119-33.
BARTtETT, M. S. (1966). An Introduction to Stochastic Processes, 2nd edition. Cambridge University Press.
BLACKMAN, R. B. (1965). Data Smoothing and Prediction. Reading, Mass.: Addison-Wesley.
BLACKMAN, R. B. & TCTKEY, J. W. (1958). The Measurement of Power Spectra. New York: Dover.
CHUNG, K. L. (1960). Markov Processes with Stationary Transition Probabilities. Heidelberg: Springer Verlag.
ERHMAN, J. R., FOSDICK, L. D. & HANDSCOMB, D. C. (1960). Computation of order parameters in an
Ising lattice by the Monte Carlo method. J. Math. Phys. 1, 547-58.
FOSDICK, L. D. (1963). Monte Carlo calculations on the Ising lattice. Methods Computational Phys.
1, 245-80.
Fox, L. & MAYERS, D. F. (1968). Computing Methods for Scientists and Engineers. Oxford: Clarendon
Press.
FRASER, D. A. S. (1968). The Structure of Inference. New York: Wiley.

Monte Carlo methods using Markov chains

109

BY, J. M. & HANDSCOMB, D. C. (1964). Monte Carlo Methods. London: Methuen.
HANDSCOMB, D. C. (1962). The Monte Carlo method in quantum statistical mechanics. Proc. Gamb.

Phil. Soc. 58, 594-8.
HANNAN, E. J. (1957). The variance of the mean of a stationary process. J. R. Statist. Soc. B 19, 282-5.
HASTSEN, M. H., HUBWITZ, W. N. & MADOW, W. G. (1953). Sample Survey Methods and Theory. New

York: Wiley.
JAMES, A. T. (1955). A generating function for averages over the orthogonal group. Proc. R. Soc. A 229,367-75.
JAMBS, A. T. (1964). Distributions of matrix variates and latent roots derived from normal samples.

Ann. Math. Statist. 35, 475-501.
JOWETT, G. H. (1955). The comparison of means of sets of observations from sections of independent

stochastic series. J. R. Statist. Soc. B 17, 208-27.
METROPOLIS, N., ROSENBLUTH, A. W., ROSENBLUTH, M. N., TELLER, A. H. & TELLEB, E. (1953).

Equations of state calculations by fast computing machines. J. Chem. Phys. 21, 1087-92.
TOCHER, K. D. (1963). The Art of Simulation. London: English Universities Press.
TBOTTEB, H. F. & TUKEY, J. W. (1956). Conditional Monte Carlo for normal samples. In Symposium

on Monte Carlo Methods, ed. H. A. Meyer, pp. 64-79. New York: Wiley.
WILK, M. B. & GNASTADESIKAN, R. (1968). Probability plotting methods for the analysis of

Biometrika 55, 1-17.

[Received February 1969. Revised July 1969]data.