Документ взят из кэша поисковой машины. Адрес оригинального документа : http://hea-www.harvard.edu/AstroStat/Stat310_fMMV/mxl_20060425.pdf Дата изменения: Wed Apr 26 22:12:06 2006 Дата индексирования: Tue Oct 2 04:28:28 2012 Кодировка: Поисковые слова: ngc 6559

MARKOV CHAIN MONTE CARLO:
A Workhorse for Modern Scientific Computation Xiao-Li Meng
Department of Statistics Harvard University
April 26, 2006 1


Introduction

April 26, 2006

2


Applications of Monte Carlo
Physics Chemistry Astronomy Biology Environment Engineering Traffic ...
April 26, 2006 3

Sociology Education Psychology Arts Linguistics History

Economics Finance Management Policy Military Government

Medical Science Business


Monte Carlo
A Chinese version of the previous slide
...
April 26, 2006 4






Monte Carlo Integration
Suppose we want to compute

where f(x) is a probability density. If we have samples x1,...,xn f(x), we can estimate I by

April 26, 2006

5


Monte Carlo Optimization
We want to maximize p(x) Simulate from f(x) p(x). As , the simulated draws will be more and more concentrated around the maximizer of p(x) =2
dens ity

=1

0.0 0.1 0.2 0.3 0.4

dens ity

-5

0 x

5

0.00 0.05 0.10 0.15

-5

0 x

5

dens ity

=20

0.0 e+00

8.0 e-09

-5

0 x

5

April 26, 2006

6


Simulating from a Distribution
What does it mean?
Suppose a random variable () X can only take two values:

Simulating from the distribution of X means that we want a collection of 0's and 1's:

such that about 25% of them are 0's and about 75%of them are 1's, when n, the simulation size is large.

The {xi, i = 1,...,n} don't have to be independent

April 26, 2006

7


Simulating from a Complex Distribution
Continuous variable X, described by a density function f(x)
f( x,y)

Complex:
the form of f(x) the dimension of x
y x

April 26, 2006

8


Markov Chain Monte Carlo
where {U(t), t=1,2,...} are identically and independently distributed.

Under regularity conditions,

So We can treat {x(t), t= N0, ..., N} as an approximate sample from f(x), the stationary/limiting distribution.

April 26, 2006

9


Gibbs Sampler
Target density: We know how to simulate form the conditional distributions For the previous example,
N(,2) Normal Distribution "Bell Curve"

April 26, 2006

10


April 26, 2006

11


Statistical Inference
Point Estimator: Variance Estimator:

Interval Estimator:

April 26, 2006

12


Gibbs Sampler (k steps)
Select an initial value (x For t = 0,1,2, ..., N
Step 1: Draw x Step 2: Draw x ... Step K:Draw xk
(t+1) (t+1) 1 (t+1) 2

1

(0)

,x

2

(0)

,..., x

k

(0)).

from f(x1|x2(t), x3(t),..., xk(t)) from f(x2|x1(t+1), x3(t),..., xk(t)) from f(xk|x1(t+1), x2(t
+1)

,..., xk-1(

t+1)

)

Output {(x1(t), x

2

(t)

,..., xk(t) ): t= 1,2,...,N}

Discard the first N0 draws Use {(x1(t), x2(t),..., xk(t) ): t= N0+1,2,...,N} as (approximate) samples from f(x1, x2,..., xk).

April 26, 2006

13


Data Augmentation
We want to simulate from

But this is just the marginal distribution of

0.0

0.1

0.2

So once we have simulations: {(x(t), y(t): t= 1,2,...,N)}, we also obtain draws: {x(t): t= 1,2,...,N)}

density

0.3

0.4

0.5

0.6

April 26, 2006

0

2 x

4

6

14


A More Complicated Example

f( x,y)

y

x

April 26, 2006

15


Metropolis-Hastings algorithm
Simulate from an approximate distribution q(z1|z2), then Step 0: Select z(0); Now for t = 1,2,...,N, repeat Step 1: draw z1 from q(z1|z2=z(t)) Step 2: Calculate

Step 3: set Discard the first N0 draws

Accept reject

April 26, 2006

16


M-H Algorithm: An Intuitive Explanation
Assume , then

April 26, 2006

17


M-H: A Terrible Implementation
We choose q(z|z2)=q(z)=(x-4)(y-4) Step 1: draw x N(4,1), y N(4,1); Dnote z1=(x,y)

Step 2: Calculate

Step 3: draw u U[0,1] Let

April 26, 2006

18


Why is it so bad?

April 26, 2006

19


M-H: A Better Implementation
Starting from some arbitrary (x(0),y(0)) Step 1: draw x N(x(t),1), y N(y(t),1) "random walk"

Step 2: dnote z1=(x,y), calculate

Step 3: draw u U[0,1] Let

April 26, 2006

20


Much Improved!

April 26, 2006

21


Further Discussion
How large should N0 and N be?

Not an easy problem!
Key difficulty: multiple modes in unknown area We would like to know all (major) modes, as well as their surrounding mass. Not just the global mode We need "automatic, Hill-climbing" algorithms. The Expectation/Maximization (EM) Algorithm, which can be viewed as a deterministic version of Gibbs Sampler.

April 26, 2006

22


Drive/Drink Safely, Don't become a Go to Graduate School, Become a

Statistic;

Statistician

!

April 26, 2006

23