Lecture Notes 25: MCMC - Metropolis Hastings Algorithm

Acceptance-Rejection Method (AR) 1 14.384 Time Series Analysis, Fall 2007 Professor Anna Mikusheva Paul Schrimpf, scribe Novemeber 29, 2007 Lecture 25 MCMC: Metropolis Hastings Algorithm A good reference is Chib and Greenberg (The American Statistician 1995). Recall that the key object in Bayesian econometrics is the posterior distribution: f (YT |θ)p(θ) p(θ|YT ) = R ˜ f (YT |θ)dθ̃ It is often difficult to compute this distribution. In particular, the integral in the denominator is difficult. So far, we have gotten around this by using conjugate priors – classes of distributions for which we know the form of the posterior. Generally, it’s easy to compute the numerator, f (YT |θ)p(θ), but it is hard to compute R ˜ the normalizing constant, the integral in the denominator, f (YT |θ)dθ̃. One approach is to try to compute this integral in some clever way. Another, more common approach is Markov Chain Monte-Carlo (MCMC). The goal here is to generate a random sample θ1 , .., θN from p(θ|YT ). We can then use moments from this sample to approximate moments of the posterior distribution. For example, 1 X θn E(θ|YT ) ≈ N There are a number of methods for generating random samples from an arbitrary distribution. Acceptance-Rejection Method (AR) We start with the simplest one. The goal is to simulate ξ ∼ π(x). What we know: 1) a function, f (x), such that π(x) = f (x) The constant k is unknown (that is, f is a k pdf up to an unknown normalization). 2)we can simulate draws from some candidate pdf h(x); 3)there is a known constant c such that f (x) ≤ ch(x) We simulate draws from π(x) as follows: 1. Draw z ∼ h(x), u ∼ U [0, 1] 2. If u ≤ f (z) ch(z) , accept the draw ξ = z. Otherwise discard the draw and repeat (1) The intuition of the procedure is the following: Let v = uch(z) and imagine the joint distribution of (v, z). It has support under the graph of ch(z) with a uniform density (it is uniform on {(v, z) : z ∈ Support(h), 0 ≤ v ≤ ch(z)}). Then, it is fairly easy to see that if we accept ξ = z, the joint distribution of (v, ξ) is uniform over the support {(v, ξ) : ξ ∈ Supportt(π), f (ξ) ≥ v ≥ 0}. Then (for the same reason that h(z) is the marginal density of (v, z)), the marginal density of ξ will be f (ξ) k . More formally, Proof. Let ρ be the probability of rejecting a single draw. Then, ¶ µ f (z1 ) P (ξ ≤ x) =P z1 ≤ x, u1 ≤ (1 + ρ + ρ2 + ...) ch(z1 ) µ ¶ · µ ¶ ¸ 1 f (z) 1 f (z1 ) |z 1{z≤x} = Ez P u ≤ = P z1 ≤ x, u1 ≤ 1−ρ ch(z1 ) 1−ρ ch(z) Z x Z x Z x 1 f (z) f (z) = h(z)dz = dz = π(z)dz 1 − ρ −∞ ch(z) c(1 − ρ) −∞ −∞ Markov Chains 2 The last line is due to the fact that there exists the unique constant that normalizes f to be a pdf. Since 1 is this constant. the left hand side is a cdf, then c(1−ρ) A major drawback of this method is that is may lead us to reject many draws before we finally accept f (z) could be very small one. This can make the procedure inefficient. If we choose c and h(z) poorly, then ch(z) for many z. It will be especially difficult to choose a good c and h() when we do not know much about π(z). Markov Chains A Markov Chain is a stochastic process where the distribution of xt+1 only depends on xt , P (xt+1 ∈ A|xt , xt−1 , ...) = P (xt+1 ∈ A|xt ) ∀A. Definition 1. A transition kernel is a function, P (x, A), such that, for every x it is a probability measure in the second argument: P (x, A) = P (xt+1 ∈ A|xt = x) It gives the probability of moving from x into the set A. The transition kernel may have atoms, in particular, we would be considering cases with non-zero probability of (not moving) staying: P (x, {x}) 6= 0. We want to study the behavior of a sequence of draws x1 → x2 → ... where we move around according to a transition kernel. Suppose the distribution of xt is P (t) , then the distribution of y = xt+1 is Z (t+1) P (y)dy = P (t) (x)P (x, dy)dx. < DefinitionR 2. A distribution π ∗ is called aninvariant measure (with respect to transition kernel P (x, A)) if π ∗ (y)dy = < π ∗ (x)P (x, dy)dx. Under some regularity conditions, a transition kernel P (x, A) has a unique invariant distribution π ∗ ; and a marginal distribution P (t) of xt - an element in Markov chain with the transitional kernel P (x, A) converges to its invariant distribution π ∗ as t → ∞. That is, if one would run a Markov chain long enough then the distribution of the draw is close to π ∗ . Generally, if the transition kernel is irreducible (it can reach any point from any other point) and aperiodic (not periodic, i.e. the greatest common denominator of {n : y can be reached from x in n steps} is 1), then it converges to an invariant distribution. A classical Markov chain problem is to find π ∗ given P (x, A). The MCMC has an inverse problem. Assume we want to simulate a draw from π ∗ (which we know up to a constant multiplier). We need to find a transition kernel P (x, dy) such that π ∗ is its invariant measure. Let’s suppose that π ∗ is continuous. We will consider the class of kernels (∗) P (x, dy) = p(x, y)dy + r(x)∆x (dy), here ∆x (dy) is a unit mass measure concentrated at point x: ∆x (A) = I{x ∈ A}. So, the transition kernel (*) says that we can stay at x with probability r(x), otherwise y is distributed according to some p R df proportional Rto p(x, y). Notice, Rthat p(x, y) isn’t exactly a density because it doesn’t integrate to 1. P (x, dy) = 1 = p(x, y)dy + r(x); p(x, y)dy = 1 − r(x). Definition 3. A transition kernel is reversible if π(x)p(x, y) = π(y)p(y, x) Theorem 4. If a transition kernel is reversible, then π is invariant. Metropolis-Hastings 3 Proof. We need to check that the definition of invariant distribution is satisfied ¶ Z Z µZ Z π(x)P (x, A)dx = p(x, y)dy π(x)dx + r(x)∆x (A)π(x)dx < < Z< Z A Z = p(x, y)π(x)dxdy + r(x)π(x)dx ZA Z< ZA = p(y, x)π(y)dxdy + r(x)π(x)dx A < A µZ ¶ Z Z = π(y) p(y, x)dx dy + r(x)π(x)dx < A ZA Z = π(y)(1 − r(y))dy + r(x)π(x)dx = π(A) A A Metropolis-Hastings The goal: we want to simulate a draw from the distribution π which we know up to a constant. That is, we can compute a function proportional to π, f (x) = kπ(x). We will generate a Markov chain with transition kernel of the form (*), that will be reversible for π. Then if the chain will run long enough the element of the chain will have distribution π. The main question is how to generate such a Markov chain? Supp R ose we have a Markov chain in state x. Assume that we can draw y ∼ q(x, y), a pdf with respect to y (so q(x, y)dy = 1). Consider using this q as a transition kernel. Notice that if π(x)q(x, y) > π(y)q(y, x) then the chain won’t be reversible (we would move from x to y too often). This suggests that rather than always moving to the new y we draw, we should only move with some probability, α(x, y). If we construct α(x, y) such that π(x)q(x, y)α(x, y) = π(y)q(y, x)α(y, x) then we will have a reversible transition kernel with invariant measure π. We can take: π(y)q(y, x) α(x, y) = min{1, } π(x)q(x, y) We can calculate α(x, y) because although we do not know π(x), we do know f (x) = kπ(x), so we can compute the ratio. In summary, the Metropolis-Hastings algorithm is: given xt we move to xt+1 by 1. Generate a draw, y, from q(xt , ·) 2. Calculate α(xt , y) 3. Draw u ∼ U [0, 1] 4. If u < α(xt , y), then xt+1 = y. Otherwise xt+1 = xt This produces a chain with P (x, dy) = q(y, x)α(y, x)dy + r(x)∆x (dy), Z r(x) = 1 − q(y, x)α(y, x)dy. Then the marginal distribution of xt will converge to π. In practice, we begin the chain at an arbitrary x0 , run the algorithm many, say M times, then use the last N < M draws as a sample from π. Note that although the marginal distribution of the xt is π, the xt are autocorrelated. This is not a problem for computing moments from the draws (although the higher the autocorrelation, the more draws we need to get the same accuracy), but if we want to put standard errors on these moments, we need to take the autocorrelation into account. Metropolis-Hastings 4 Choice of q() • Random walk chain: q(x, y) = q1 (y − x), i.e. y = x + ², ² ∼ q1 . This can be a nice choice because q (x,y) π (y) if q1 is symmetric, q1 (z) = q1 (−z), then q(y,x) drops out of α(x, y) = min{1, π(x) }. Popular such q1 are normal and U [−a, a]. Note that there is a tradeoff between step-size in the chain and rejection probability when choosing σ 2 = E²2 . Choosing σ 2 too large will lead to many draws of y from low probability areas (low π), and as a result we will reject lots of draws. Choosing σ 2 too small will lead us to accept most draws, but not move very much, and we will have difficulty covering the whole support of π. In either case, the autocorrelation in our draws will be very high and we’ll need more draws to get a good sample from π. • Independence chain: q(x, y) = q1 (y) • If there is an additional information that π(y) ∝ ψ(y)h(y) where ψ is bounded and we can sample from ψ(y) q(x, y) = h(y). This also simplifies α(x, y) = min{1, ψ(x) } • Autocorrelated y = a + B(x − a) + ² with B < 0, this leads to negative autocorrelation in y. The hope is that this reverses some of the positive autocorrelation inherent in the procedure.