ĐẠI HỌC QUỐC GIA HÀ NỘI
TRƯỜNG ĐẠI HỌC KHOA HỌC TỰ NHIÊN
TRẦN THỊ BÍCH NGỌC
PHƯƠNG PHÁP MCMC
VÀ MỘT SỐ ỨNG DỤNG
TÓM TẮT LUẬN VĂN THẠC SỸ TOÁN HỌC
Chuyên ngành : LÝ THUYẾT XÁC SUẤT VÀ THỐNG KÊ TOÁN HỌC
Mã số : 60 46 01 06
HÀ NỘI, 2014
LỜI MỞ ĐẦU
Luận văn này với mục đích trình bày về phương pháp MCMC và một
số ứng dụng của nó.Luận văn được xây dựng dựa trên lý thuyết về suy
luận Bayes,tích phân Monte Carlo và xích Markov
Luận văn gồm có 4 chương:
Chương 1. Tổng quan.
Suy luận Bayes: giới thiệu về suy luận Bayes, các đặc điểm của mô hình
Bayes, các tiên nghiệm Jeffreys.
Tích phần Monte Carlo: Bài toán tích phân Monte Carlo, xấp xỉ Monte
Carlo, Monte Carlo thông qua lấy mẫu theo trọng số.
Phương pháp sinh biến ngẫu nhiên: Phương pháp biến đổi, phương pháp
chấp nhận - bác bỏ, phương pháp tỷ số đều.
Xích Markov: Các định nghĩa và kí hiệu, Sự hội tụ của các phân phối,
giới hạn của giá trị trung bình.
Chương 2. Mẫu Gibbs.
Giới thiệu về phương pháp lấy mẫu Gibbs và ví dụ cho trường hợp biến
ngẫu nhiên nhiều chiều.
Thuật toán mở rộng dữ liệu:mô tả thuật toán và một số ví dụ tương
ứng.
Chương 3. Thuật toán Metropolis- Hastings.
Thuật toán Metropolis- Hasting: Khái niệm, mẫu độc lập, xích bước
ngẫu nhiên.

là mẫu lấy từ Bernoulli (θ) với không gian mẫu X = {0, 1}
và hàm khối xác suất
Pr (X = 1 |θ) = θ;
Pr (X = 0 |θ) = 1 − θ
(1.1)
trong đó X là biến ngẫu nhiên Bernoulli với X = 1 nếu thành công, và
X = 0 nếu thất bại.
Ta viết N =

n
i=1
x
i
là số quan sát thành công trong n phép thử
Bernoulli.
Khi đó N |θ ∼ B (n, θ) là phân phối nhị thức với cỡ n và xác suất thành
công θ. Xác suất nghịch đảo của θ cho bởi x
1
, x
2
, , x
n
được hiểu như phân
phối hậu nghiệm,được xem như là phân phối Beta, Beta(1+N,1+n-N) với
hàm mật độ xác suất
1
B(1 + N, 1 + n −N)
θ
(1+N)−1
(1 −θ)

ở đó
L (θ |X ) ∝ f (X |θ)
trong đó δ được gọi là thống kê hợp lý của δ với X đã cho.
1.1.2 Các tiên nghiệm Jeffreys
Đối với các trường hợp khi thông tin này là không sẵn có hoặc không
dễ xác định bằng một phân phối xác suất chính xác, đặc biệt là đối với
các bài toán với số chiều cao, khi đó phương pháp thường được sử dụng là
phương pháp Jeffreys, với việc giả thiết tiên nghiệm có dạng:
π
J
(θ) ∝ |I (θ)|
1
2
(θ ∈ Θ)
(1.6)
Trong đó I (θ) là lượng thông tin Fisher. Phân phối hậu nghiệm tương ứng
của θ cho bởi X như sau:
π
J
(µ |X ) = N (X, 1) (1.7)
5
1.2 Tích phân Monte Carlo
1.2.1 Bài toán
Cho ν là độ đo xác suất trên σ - trường Borel X với không gian mẫu
X ⊆ R
d
, trong đó R
d
là không gian Euclide d-chiều. Một khó khăn thường
gặp trong bài toán là ước tính tích phân dạng:

i
) (1.10)
có thể được sử dụng để tính xấp xỉ (1.10) vì h
n
hội tụ tới (1.10) hầu chắc
chắn theo luật số lớn. Khi h (X) có phương sai hữu hạn, sai số của xấp xỉ
này có thể được mô tả bằng định lý giới hạn trung tâm, nghĩa là:
h
n
− E
f
[h (X)]

nV ar (h (X))
∼ N (0, 1)
Tương tự V ar (h (X)) có thể được xấp xỉ bằng phương sai mẫu:
1
n −1
n

i=1

h (X
1
) −h
n

2
Phương pháp xấp xỉ tích phân qua các mẫu mô phỏng được biết đến như
là phương pháp Monte Carlo

1.3 Phương pháp sinh biến ngẫu nhiên
Thuật toán 1.1 (Hàm phân bố ngược liên tục)
1, Sinh ra một biến ngẫu nhiên đều U.
2, Tính toán và đưa ra kết quả X = F
−1
(U) trong đó F
−1
(.) là hàm số
ngược của hàm phân bố liên tục F (.).
Thuật toán 1.2 (Hàm phân bố ngược rời rạc)
1, Sinh ra biến ngẫu nhiên đều U.
2, Tìm X thỏa mãn F (X − 1) < U ≤ F (X).
3, Trả lại giá trị X.
1.3.1 Phương pháp biến đổi
Các phương pháp biến đổi tốt hơn thu được bằng cách dựa vào phân
phối mục tiêu f(x).
Công thức Phép biến đổi Phân phối
Mũ X = −ln(U) X ∼ Expo(1)
Cauchy X = tan (πU − π/2)) X ∼ Cauchy(0, 1)
Beta X
i
ind
∼
Gamma (α
i
) , i = 1, 2
X
1
X
1

tiếp từ C
h
,ta có thể lấy mẫu một cách gián tiếp qua C
h
như sau:
(i) Sinh ra những điểm có tính đều trên một miền mở rộng và dễ dàng để
lấy mẫu D ⊇ C
h
và
(ii) Thu thập những điểm thuộc vào miền C
h
. Miền mở rộng D như vậy
có thể được xây dựng bằng một phân phối có thể lấy mẫu một cách
đơn giản với hàm mật độ g (x) thoả mãn
f(x)
g(x)
bị chặn trên bởi một số
hằng số hữu hạn M. Vì vậy C
h
là đóng trong miền:
C
g
= {(x, u) : 0 ≤ u ≤ g (x)} ⊂ R
d+1
(1.13)
với h (x) ∝ f (x). Phân phối g (x) được gọi là phân phối công cụ và
f (x) là phân phối mục tiêu. Tóm lại, thuật toán AR dùng để sinh các
số ngẫu nhiên từ f (x) bằng cách sử dụng phân phối công cụ g (x),
trong đó :
sup

1.3.3 Phương pháp tỷ số đều
Thuật toán 1.5.
Lặp lại hai bước sau cho đến khi giá trị trả về trong bước 2:
1, Sinh (Y, Z) có độ lệch đều trên miền D ⊇ C
(Y,Z)
h
2, Nếu (Y, Z) ∈ C
(Y,Z)
h
, trả về giá trị X = x(Y, Z) là độ lệch mong
muốn.
Thuật toán này có tỉ số chấp nhận
r =

C
(Y,Z)
h
dydz

D
dydz
=

X
h (x) dx
J
x,u
(z, y)




∂x
∂z




là hệ số Jacobi của các phép biến đổi.
1.4 Xích Markov
Xích Markov là một dãy các biến ngẫu nhiên {X
i
, i = 0, 1, 2 } với tính
Markov được cho bởi trạng thái hiện tại, trạng thái tương lai, trạng thái
quá khứ là độc lập, nghĩa là với mọi tập đo được A ⊆ X:
Pr (X
t+1
∈ A |X
0
= x
0
, , X
t
= x
t
) = Pr (X
t+1
∈ A |X
t
= x
t

∈ Bi.o. |X
0
= x) = 1 với hầu hết π (x)
(b) Xích là hồi quy Harris nếu Pr (X
n
∈ Bi.o. |X
0
= x) = 1 với hầu hết
π(x)
Định nghĩa 1.2. Các dạng ergodic khác nhau được cho như sau:
(a) Một xích Markov được gọi là ergodic nếu nó là Harris dương hồi quy
và không tuần hoàn.
(b) Cho H
B
là thời điểm chạm của tập B. Một xích ergodic với phân phối
dừng π (x) được gọi là ergodic cấp 2 nếu:

B
E
x

H
2
B

π (dx) < ∞
với mọi H ∈ X thỏa mãn π (H) > 0
(c) Một xích ergodic với phân phối dừng π (x) được gọi là ergodic hình học
nếu tồn tại một hàm số thực không âm M thỏa mãn E (|M (X)|) < ∞
và một hằng số dương r < 1 sao cho:

→ E
f
(h (X)) h.c.c
.
Định lý 1.3. Giả sử rằng X
n
là ergodic bậc 2 với phân phối cân bằng
f (x) và giả sử h (x) có giá trị thực và bị chặn. Khi đó tồn tại một số thực
σ
h
sao cho phân phối của
√
n

h
n
− E
f
(h (X))

hội tụ yếu tới phân phối
chuẩn với kỳ vọng bằng 0 và phương sai σ
2
h
với mọi phân phối ban đầu.
Giả thiết về tính bị chặn của h(x) có thể được bỏ nếu xích là ergodic
đều và E
f

h

h
với mọi phân phối ban đầu.
11
Chương 2
MẪU GIBBS
Chương này đề cập đến thuật toán đơn giản nhất của phương pháp
MCMC: thuật toán lấy mẫu Gibbs và một số ứng dụng của nó.
2.1 Mẫu Gibbs
Giả sử rằng ta muốn sinh các số ngẫu nhiên từ hàm mật độ mục tiêu
f (x), x ∈ X ⊆ R
d
. Ta tiến hành phân hoạch vector d-chiều x vào K khối
và viết x = (x
1
, , x
K
)

trong đó K ≤ d và dim (x
1
) + dim (x
K
) = d
với dim (x
k
) là số chiều của x
k
.
Ta kí hiệu
f


x
j
1
, , x
j
k−1
, y
j
k+1
, , y
j
K

f
j
k

y
j
k


x
j
1
, , x
j
k−1
, y

1, Sinh ra x
(t)
1
∼ f
1

x
1



x
(t−1)
2
, , x
(t−1)
K

.
.
.
k, Sinh ra
x
(t)
k
∼ f
k

x
k


x
(t)
1
, , x
(t)
K−1

Ví dụ 2.1. (Phân phối chuẩn của biến ngẫu nhiên nhiều chiều)
Để minh hoạ cho mẫu Gibbs, ta dùng phân phối chuẩn hai chiều
p (x) = N

µ,


; µ = [µ
1
, µ
2
] = [0, 0] ;

=

1 ρ
12
ρ
21
1

=


µ
2
+ ρ
12

x
(t)
1
− µ
1

,

1 −ρ
2
12

13
Hình 2.1: Mẫu Gibbs đối với phân phối chuẩn hai chiều
2.2 Thuật toán mở rộng dữ liệu
Thuật toán DA: mẫu Gibbs hai bước
Lấy θ
(0)
∈ Θ và lặp lại với t = 1, 2,
Bước I Chỉ ra X
(t)
mis
∼ f
mis

ngược và một số ứng dụng của thuật toán này trong việc xác định điểm
thay đổi
3.1 Thuật toán Metropolis – Hastings
3.1.1 Khái niệm
Thuật toán lấy mẫu Metropolis (hay mẫu Metropolis) có thể tóm tắt
như sau:
Định nghĩa 3.1. Mẫu Metropolis
1, Sinh ra y từ q (y|x
t
)
2, Tính toán tỷ số chấp nhận
α (x
t
, y) = min

1,
f (y)
f (x
t
)

Đặt
x
t+1
= y
với xác suất α (x
t
, y) và x
t+1
= x

= y
với xác suất α (x
t
, y) và x
t+1
= x
t
với xác suất 1 −α (x
t
, y)
3.1.2 Mẫu độc lập
Mẫu độc lập có thể xem xét như là dạng tổng quát của thuật toán chấp
nhận – bác bỏ. Dễ dàng nhận thấy rằng xích độc lập là bất khả quy và
không tuần hoàn nếu: {x : x ∈ X, f (x) > 0} ⊆ {x : x ∈ X, g (x) > 0}
Định lý 3.1. Xích độc lập là ergodic đều nếu tồn tại một hằng số M sao
cho
f (x) ≤ Mg (x)
(x ∈ {x : f(x) > 0})
3.1.3 Xích bước ngẫu nhiên
Xích bước ngẫu nhiên được tạo nên bằng cách lấy phân phối có điều
kiện có dạng:
q (x, y) = q (y −x)
Nghĩa là, bước nhảy đề xuất có hướng và khoảng cách xuất phát từ trạng
thái hiện tại x
t
là độc lập với x
t
Những phân phối hình cầu đơn giản, ví
dụ như phân phối chuẩn tắc, phân phối Student, phân phối đều với hình
cầu có tâm tại O, phân phối ellip thường là những lựa chọn phổ biến nhất

x


x
(t−1)

Xác suất chấp nhận α = min

1,
p(x
∗
)q
(
x
(t−1)
|x
∗
)
p
(
x
(t−1)
)
q
(
x
∗
|
x
(t−1)

Ta giới hạn khoảng θ
1
và θ
2
là [0,8] và đặt λ
1
= 0, 5; λ
2
= 0, 1; λ +
0, 01; max(λ
1
, λ
2
) = 8
17
Hình 3.1: Thuật toán MH đối với cập nhật từng khối đối với phân phối mũ hai chiều
Ví dụ 3.2. Phân phối chuẩn hai chiều
p (x) = N

µ,


; µ = [µ
1
, µ
2
] = [0, 0] ;

=


2
21

x
(t)
2
∼ N

µ
2
+ ρ
12

x
(t)
1
− µ
1

,

1 −ρ
2
12

Hình bên dưới minh hoạ phân phối mẫu (bên trái) và phân phối mục
tiêu(bên phải)thông qua thuật toán MH đối với cập nhật từng khối.Chúng
ta có thể thấy rằng cập nhật từng khối mô tả khá tốt việc tạo ra các mẫu
từ phân phối mục tiêu.
18

j

q

x
(t−1)
i
|x
∗
i

p

x
(t−1)
i
, x
(t−1)
j

q

x
∗
i



x
(t−1)

] = [0, 0] ;
19

=

1 ρ
12
ρ
21
1

=

1 0, 8
0, 8 1

x
(t)
1
∼ N

µ
1
+ ρ
21

x
(t−1)
2
− µ

Hình 3.3: Thuật toán MH với cập nhật từng thành phần đối với phân phối chuẩn hai
chiều
3.3 Các dạng khác nhau của thuật toán Metropolis
- Hastings
3.3.1 Thuật toán chạm và chạy
Định nghĩa 3.3. Thuật toán nhấn và chạy
1, Sinh ra
d ∼ g (d)
(d ∈ O)
và
λ ∼ l (λ |d, x)
trên X
x,d
và tính toán một xác suất chấp nhận MH α (x, y) trong đó
x = x
(t)
20
2, Sinh ra U từ Unif (0, 1) và đặt:
X
(t+1)
=

x + λd, nếuU ≤ α (x, y)
x, nếu ngược lại
3.3.2 Thuật toán Langevin
Định nghĩa 3.4. Thuật toán Langevin
1. Chỉ ra một phương trình mới:
x
∗
= x

2
∇log f (x
∗
)



2
/2σ
2

f

x
(t)

exp

−



x
∗
− x
(t)
−
σ
2
2

= ω (y
i
, x)
với i = 1, 2, , k
2. Chọn y = y
j
từ {y
1
, y
2
, , y
k
} theo xác suất tỷ lệ với ω
i
, i = 1, 2, , k.
Sinh ra x
∗
1
, , x
∗
k−1
từ q (. |y ). Đặt x
∗
k
= x và tính toán ω
∗
i
= ω (x
∗
i

1. Chọn mẫu M
k
∗
với xác suất q

k
(t)
, k
∗

2. Đưa ra u
1
, , u
s
∼ ψ
k
(t)
→k
∗
(u)
3. Đặt (θ
k
∗
, u
∗
) = T

θ
(t)
k

(t)
k
|Y

q

k
(t)
, k
∗

ψ
k
(t)
→k
∗
(u)







∂ (θ
∗
k
∗
, u
∗

,u
)
là Jacobi của phép biến đổi (3.10)
5. Đặt X
(t+1)
= (k
∗
, θ
∗
k
∗
) với xác suất min (1, r) và X
(t+1)
= X
t
với xác
suất còn lại
3.4.2 Xác định điểm thay đổi
Xét ứng dụng sau đây của RJMCMC cho bài toán xác định điểm thay
đổi. Đặt Z = (z
1
, , z
n
) là dãy quan sát độc lập. Đặt ϑ = (ϑ
1
, , ϑ
n−1
) là
chỉ số của điểm thay đổi, một vector nhị phân với ϑ
c

r
(.) là phân phối Gauss với các tham số
µ
r
, σ
2
r
chưa biết.
Đặt ϑ
(k,l)
là mẫu thứ l sinh ra tại bước lặp t, trong đó k chỉ số điểm
thay đổi của mẫu. Các mẫu tiếp theo có thể được tạo ra theo các bước sau
đây:
22
a, Đặt j = k −1 j = k hoặc j = k + 1 theo xác suất q
k,j
trong đó q
k,k
=
1
3
với k
min
≤ k ≤ k
max
q
k
min
,k
min

(k,l)
bằng dịch chuyển ’tử’
23
Chương 4
Phương pháp biến phụ trợ MCMC
Trong chương này ta xét sự tồn tại các công thức phụ trợ MCMC.
4.1 Mô phỏng nhiệt luyện
Định nghĩa 4.1. Thuật toán mô phỏng nhiệt luỵện
1. Khởi tạo mô phỏng tại nhiệt độ T
1
và một mẫu bất kỳ x
0
2. Tại mỗi nhiệt độ T
i
, mô phỏng của phân phối f (x, T
i
) với N
i
bước lặp
sử dụng một mẫu MCMC. Thông qua mẫu cuối cùng tới mức nhiệt
độ thấp hơn tiếp theo như là mẫu khởi tạo.
4.2 Mô phỏng điều hoà nhiệt
Định nghĩa 4.2. Mô phỏng điều hoà nhiệt
1. Sinh ra một số ngẫu nhiên U ∼ Uniform [0, 1] và xác định giá trị
của j theo ma trận truyền đề nghị (q
ij
)
2. Nếu j = i
t
đặt i

−H (x)

1
T
j
−
1
T
i
t

q
j,i
t
q
i
t
,j

trong đó

Z
i
chỉ ước lượng của Z
i
. Nếu nó được chấp nhận đặt i
t+1
= j.
Ngược lại, đặt i
t+1

) f

y





θ

q (θ
t
|θ

) f (y |θ
t
)
f (x |θ
t
) f (θ
t
) f

y




θ



∼ q (θ

|θ, x)
2. Sinh ra một biến phụ trợ y ∼ f (y |θ

) với xác suất min {1, r (θ, θ

, y |x)}
trong đó:
r (θ, θ

, y |x) =
π (θ

) f (x |θ

) f (y |θ) q (θ |θ

)
π (θ) f (x |θ) f (y |θ

) q (θ

|θ)
(4.1)
25