Luận văn Thạc sĩ Khoa học: Thuật toán mô phỏng MCMC thích nghi và ứng dụng

Chia sẻ: Na Na | Ngày: | Loại File: PDF | Số trang:70

Thêm vào BST

Báo xấu

71
lượt xem 8
download

Download Vui lòng tải xuống để xem tài liệu đầy đủ

Mục đích chính của luận văn này là trình bày các phương pháp MCMC cơ bản và hai thuật toán MCMC thích nghi từ bài báo. Đồng thời đưa ra các so sánh giữa các thuật toán MCMC và chứng minh chi tiết các định lý trong bài báo cũng như đưa ra một số ứng dụng của thuật toán.

Chủ đề:

Bình luận(0) Đăng nhập để gửi bình luận!

Lưu

Nội dung Text: Luận văn Thạc sĩ Khoa học: Thuật toán mô phỏng MCMC thích nghi và ứng dụng

ĐẠI HỌC QUỐC GIA HÀ NỘI TRƯỜNG ĐẠI HỌC KHOA HỌC TỰ NHIÊN ------------------- NGUYỄN VĂN TÂN THUẬT TOÁN MÔ PHỎNG MCMC THÍCH NGHI VÀ ỨNG DỤNG Chuyên ngành: Lý thuyết xác suất và thống kê toán học Mã số: 60460106 LUẬN VĂN THẠC SĨ KHOA HỌC NGƯỜI HƯỚNG DẪN KHOA HỌC: TS. TRẦN MẠNH CƯỜNG Hà Nội - 2015
Mục lục Lời nói đầu 3 1 Kiến thức chuẩn bị 5 1.1 Sự hội tụ của dãy đại lượng ngẫu nhiên . . . . . . . . . . . 5 1.2 Dãy mixingale . . . . . . . . . . . . . . . . . . . . . . . . . 6 1.3 Các thuật toán mô phỏng cơ bản . . . . . . . . . . . . . . . 7 1.3.1 Phương pháp biến đổi nghịch đảo . . . . . . . . . . 8 1.3.2 Phương pháp loại bỏ . . . . . . . . . . . . . . . . . 9 1.3.3 Phương pháp lấy mẫu quan trọng . . . . . . . . . . 13 1.4 Xích Markov . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2 Phương pháp MCMC 22 2.1 Giới thiệu . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 2.2 Mẫu Metropolis - Hastings . . . . . . . . . . . . . . . . . . 23 2.3 Một số thuật toán MCMC . . . . . . . . . . . . . . . . . . 29 2.3.1 Mẫu Gibbs . . . . . . . . . . . . . . . . . . . . . . . 29 2.3.2 Mẫu độc lập . . . . . . . . . . . . . . . . . . . . . . 30 2.3.3 Mẫu Metropolis - Hastings du động ngẫu nhiên . . 32 2.3.4 Mẫu Metropolis (thành phần đơn) . . . . . . . . . . 33 3 MCMC thích nghi 34 3.1 Thuật toán Metropolis du động ngẫu nhiên thích nghi . . . 35 3.1.1 Mô tả thuật toán . . . . . . . . . . . . . . . . . . . 35 3.1.2 Tính chất ergodic . . . . . . . . . . . . . . . . . . . 37 3.1.3 So sánh các thuật toán Metropolis với thuật toán AP 38 1
3.2 Thuật toán Metropolis thích nghi . . . . . . . . . . . . . . 42 3.2.1 Mô tả thuật toán . . . . . . . . . . . . . . . . . . . 45 3.2.2 Tính Ergodic . . . . . . . . . . . . . . . . . . . . . . 47 3.2.3 So sánh các thuật toán Metropolis với thuật toán AM 59 3.3 Một số ứng dụng của MCMC thích nghi . . . . . . . . . . . 59 3.3.1 Mô hình mô phỏng GOMOS . . . . . . . . . . . . . 60 3.3.2 Mô hình suy giảm oxy . . . . . . . . . . . . . . . . . 65 Kết quả chính 67 Tài liệu tham khảo 68 2
Lời nói đầu Để tìm hiểu về MC, ta xét bài toán sau: Giả sử ta cần tính tích phân R1 0 h(x)dx. Theo định lý Newton - Leibnitz, nếu F (x) là một nguyên hàm của h(x) thì
1 I = F (x)
= F (1) − F (0).
0 Tuy nhiên, trong nhiều trường hợp, ta không thể tìm được F(x). Giả sử f (x) là hàm mật độ trên [0, 1] sao cho nếu h(x) 6= 0 thì f (x) > 0. Ta viết R1 lại I = 0 fh(x) (x) f (x)dx. Khi đó, chúng ta lấy mẫu độc lập cùng phân phối (1) (n) (x , ..., x ) từ phân phối xác định bởi mật độ f và xét: n ˆ 1X In = h(x(i) )/f (x(i) ). n i=1 Luật số lớn cho ta thấy rằng Iˆn hội tụ với xác suất 1 tới tích phân I khi n tiến tới ∞ nghĩa là Iˆn → I(h.c.c). Như vậy để tính xấp xỉ I, ta phải thực hiện n mô phỏng cho biến ngẫu nhiên X. Các mô phỏng MC cơ bản này có ưu điểm là dễ thực hiện. Tuy nhiên, nó chỉ mô phỏng được đối với các trường hợp đơn giản. Trong nhiều trường hợp phức tạp như số chiều tăng lên (phân phối nhiều chiều) ... thì các MC cơ bản không thể thực hiện được. Đề giải quyết vấn đề này, chúng ta đưa ra một phương pháp gọi là phương pháp MCMC. Ý tưởng chính của phương pháp MCMC là đi xây dựng một xích Markov có tính ergodic mà phân phối dừng là π . Khi đó, chúng ta chạy X lên đến thời gian dài N và ước lượng E(h(Y )) bởi N1 N P n=1 h(Xn ). Định lý ergodic cho ta biết với N đủ lớn, ước lượng trên sẽ gần đến E(h(Y )). Chúng ta thấy rằng việc chọn lựa phân phối đề xuất là quan trọng cho 3
sự hội tụ của thuật toán MCMC. Việc chọn lựa được phân phối đề xuất tốt thường khó thực hiện vì thông tin về mật độ mục tiêu là không có hoặc rất ít. Hơn nữa, trong thuật toán MCMC, phân phối đề xuất được chọn cho mọi bước mô phỏng. Để sử dụng các thông tin đã thu được trong các bước mô phỏng trước để mô phỏng cho bước tiếp theo, chúng ta đưa ra thuật toán MCMC thích nghi. Ở đó, phân phối đề xuất được cập nhật cùng quá trình sử dụng thông tin đầy đủ tích lũy cho đến thời điểm hiện tại. Mỗi lựa chọn phân phối đề xuất thích nghi sẽ cho chúng ta một dạng MCMC thích nghi. Mục đích chính của luận văn này là trình bày các phương pháp MCMC cơ bản và hai thuật toán MCMC thích nghi từ bài báo [6], [7]. Đồng thời đưa ra các so sánh giữa các thuật toán MCMC và chứng minh chi tiết các định lý trong bài báo cũng như đưa ra một số ứng dụng của thuật toán. Luận văn gồm 3 chương. • Chương 1 nhắc lại một số kiến thức bổ trợ về sự hội tụ của dãy đại lượng ngẫu nhiên, dãy mixingale, các thuật toán mô phỏng MC cơ bản và xích Markov. • Chương 2 trình bày về các phương pháp MCMC cơ bản. • Chương 3 trình bày chi tiết về hai phương pháp MCMC thích nghi từ hai bài báo [6] và [7]. Đó là thuật toán Metropolis du động ngẫu nhiên thích nghi ([6]) và thuật toán Metropolis thích nghi ([7]). Chỉ ra tính hội tụ của hai thuật toán và chứng minh tính ergodic của thuật toán Metropolis thích nghi. Sau mỗi thuật toán đều đưa ra sự so sánh giữa các thuật toán MCMC. Đồng thời đưa ra một số ứng dụng thực tế của mô hình MCMC thích nghi. Lời đầu tiên, xin chân thành cảm ơn thầy TS. Trần Mạnh Cường đã nhận hướng dẫn và tận tình giúp đỡ tôi hoàn thành luận văn này. Lòng biết ơn sâu sắc tôi cũng xin được gửi đến các thầy cô trong Trường ĐHKHTN - ĐHQGHN, Khoa Toán - Cơ - Tin đã giúp đỡ tôi hoàn thành khóa học. Hà Nội tháng 12 năm 2015 4
Chương 1 Kiến thức chuẩn bị 1.1 Sự hội tụ của dãy đại lượng ngẫu nhiên Giả sử (Ω, F, P ) là không gian xác suất. Định nghĩa 1.1. Một dãy các đại lượng ngẫu nhiên hay biến ngẫu nhiên (Xn ) được gọi là hội tụ hầu chắc chắn đến biến ngẫu nhiên X nếu: P {ω ∈ Ω : lim Xn (ω) 6= X(ω)} = 0. n→∞ Ký hiệu là limn→∞ Xn = X(h.c.c). Định nghĩa 1.2. Cho dãy (Xn ) các biến ngẫu nhiên. Fn (x), F (x) tương ứng là hàm phân phối của Xn , X . Gọi C(F ) là tập các điểm liên tục của hàm F . Ta nói dãy (Xn ) hội tụ theo phân phối đến X nếu ∀x ∈ C(F ), ta có: lim Fn (x) = F (x). n→∞ d Ký hiệu là Xn → − X. Định nghĩa 1.3. Một dãy các biến ngẫu nhiên (Xn ) được gọi là hội tụ theo xác suất đến biến ngẫu nhiên X nếu ∀ε > 0 ta có : P {ω ∈ Ω : |Xn (ω) − X(ω)| > ε} = 0. P Ký hiệu là Xn − → X. 5
Định nghĩa 1.4. Một dãy các biến ngẫu nhiên (Xn ) được gọi là hội tụ theo trung bình bậc r đến biến ngẫu nhiên X nếu r ≥ 1, E|Xn |r < ∞ ∀n, E|X|r < ∞ và : lim E{|Xn − X|r } = 0. n→∞ r L Ký hiệu là Xn −→ X . Định nghĩa 1.5. (luật số lớn) Cho dãy (Xn ) các biến ngẫu nhiên độc lập cùng phân phối, có cùng kỳ vọng EXi = µ (i = 1, 2, ...). Đặt Sn = X1 +...+Xn n . Ta nói dãy (Xn ) tuân theo luật số lớn nếu Sn sẽ hội tụ theo xác suất đến µ. Định lí 1.6. (định lý giới hạn trung tâm) Cho dãy (Xn ) các biến ngẫu nhiên độc lập cùng phân phối, có cùng kỳ vọng EXi = µ và phương sai √ n −nµ . Khi đó Zn sẽ hội tụ DXi = σ 2 (i = 1, 2, ...). Đặt Zn = X1 +...+X σ n theo phân phối đến biến ngẫu nhiên Z có phân phối chuẩn tắc. 1.2 Dãy mixingale Định nghĩa 1.7. Cho dãy (Xn )n≥1 các biến ngẫu nhiên bình phương khả tích trong không gian xác suất (Ω, F, P ) và dãy (Fn )+∞ n=−∞ là dãy tăng các σ - đại số con của F . Khi đó, (Xn , Fn ) được gọi là dãy mixingale nếu với mọi dãy hằng không âm cn và ψm , trong đó ψm → 0 khi m → ∞, ta có: ||E(Xn |Fn−m )||2 ≤ ψm cn và ||Xn − E(Xn |Fn+m )||2 ≤ ψm+1 cn , với mọi n ≥ 1 và m ≥ 0. Định lí 1.8. [4, tr. 41] Nếu {Xn , Fn } là một mixingale và {bn } là một dãy hằng dương tăng đến ∞ sao cho ∞ X b−2 2 n cn < ∞ và ψn = O(n −1/2 (logn)−2 ) khi n → ∞ n=1 Pn thì b−1 n i=1 Xi → 0(h.c.c). 6
1.3 Các thuật toán mô phỏng cơ bản Các kết quả thống kê thường liên quan đến tích phân. Nhắc lại rằng cả kỳ vọng và xác suất đều nhận được từ tích phân (hoặc tổng). Vì vậy, xét tích phân sau: Z 1 I= h(x)dx 0 Thông thường, người ta tiếp cận dạng tổng Riemann. Chúng ta đánh giá hàm h(x) tại n điểm (x(1) , ..., x(n) ) trong một lưới chính quy và sau đó tính: n 1X I≈ h(x(i) ). n i=1 Tuy nhiên, trong nhiều trường hợp, việc xác định lấy các điểm (x(1) , ..., x(n) ) là không thể hoặc chi phí quá tốn kém, người ta đã đưa ra một cách tiếp cận khác. Đó là quá trình Monte Carlo. Chúng ta bắt đầu bằng việc viết lại tích phân như sau: Z 1 h(x) I= f (x)dx 0 f (x) trong đó f (x) là một mật độ trên [0, 1] sao cho nếu h(x) 6= 0 thì f (x) > 0. Nhưng điều này nghĩa là: I = Ef (h(X)/f (X)), trong đó Ef là ký hiệu của kỳ vọng đối với phân phối xác định bởi f . Bây giờ, chúng ta lấy mẫu độc lập cùng phân phối (x(1) , ..., x(n) ) từ phân phối xác định bởi mật độ f và xét: n 1X Iˆn = h(x(i) )/f (x(i) ). n i=1 Luật số lớn cho ta thấy rằng Iˆn hội tụ với xác suất 1 tới tích phân I khi n tiến tới ∞ nghĩa là Iˆn → I(h.c.c). Hơn nữa, định lý giới hạn trung tâm chỉ ra rằng q (Iˆn − I)/ V ar(Iˆn ) 7
xấp xỉ phân phối chuẩn. Vì vậy phương sai V ar(Iˆn ) cho ta biết về độ chính xác ước lượng của chúng ta và nó có thể được ước lượng như sau: n 1 X vn = (h(xj )/f (xj ) − Iˆn )2 . n(n − 1) j=1 1.3.1 Phương pháp biến đổi nghịch đảo Định lí 1.9. Xét hàm phân phối lũy tích (cdf) F (x). Gọi F −1 là nghịch đảo mở rộng của F , tức là: F −1 (u) = min{x ∈ S : F (x) ≥ u} u ∈ (0, 1] Gọi U là một biến ngẫu nhiên phân phối đều (0, 1) và đặt X = F −1 (U ), khi đó phân phối của X có cdf F (x). (Chú ý rằng đối với hàm phân phối liên tục thì nghịch đảo mở rộng là nghịch đảo thông thường). Bằng định nghĩa của nghịch đảo mở rộng và tính đơn điệu của F , ta có: P (X ≤ x) = P(F −1 (U ) ≤ x) = P (U ≤ F (x)) = F (x). Ví dụ 1.1. Mô phỏng một biến ngẫu nhiên phân phối mũ với tham số λ . Một biến ngẫu nhiên có phân phối mũ với tham số λ có hàm phân phối là: F (x) = 1 − exp(−λx) với x ≥ 0. Gọi U ∼ U (0, 1) (phân phối đều trên (0, 1)) và đặt 1 Y = − log(1 − U ). λ Khi đó Y có phân phối mũ với tham số λ. Điều này có thể đơn giản hóa hơn bằng cách thừa nhận rằng 1 − U cũng là phân phối đều trên (0, 1) và vì thế 1 Y = − log(U ) λ có phân phối mũ với tham số λ. 8
Ví dụ 1.2. Mô phỏng biến ngẫu nhiên có phân phối Bernoulli (p) và biến ngẫu nhiên có phân phối nhị thức B(n, p) Cho U là một biến ngẫu nhiên phân phối đều (0, 1). Nếu ta xét  1 nếu U < p X= 0 ngược lại thì X là biến ngẫu nhiên có phân phối Bernoulli với xác suất thành công p. Cho X1 , ..., Xn là một mẫu độc lập cùng phân phối Bernoulli(p). Khi đó Y = ni=1 Xi có phân phối nhị thức B(n, p). P Ví dụ 1.3. Mô phỏng biến ngẫu nhiên tuân theo phân phối hình học (p) Giả sử X nhận giá trị trong N và P(X = j) = pj . Khi đó: j X −1 F (u) = min{j ∈ N : u ≤ pi }. i=1 Bây giờ, nếu X ∼ G(p) thì P(X > j) = (1 − p)j . Do đó j X pi = 1 − (1 − p)j ≥ u i=1 nếu và chỉ nếu log(1 − u) j≥ . log(1 − p) h i log(U ) Ký hiệu [a] là phần nguyên của a thì X = log(1−p) tuân theo phân phối hình học G(p). 1.3.2 Phương pháp loại bỏ Giả sử chúng ta muốn lấy mẫu X là một biến ngẫu nhiên liên tục với hàm mật độ f (x). Chúng ta không biết cách lấy mẫu từ X nhưng chúng ta biết cách lấy mẫu từ một biến ngẫu nhiên Y tương tự với hàm mật độ g(y). Gọi giá của f là supp(f ) = {x : f (x) > 0}. Nếu ta có supp(f ) ⊆ supp(g) 9
và f (x)/g(x) ≤ M ∀x thì ta có thể lấy mẫu từ Y để tạo ra mẫu cho X . Chúng ta lặp lại các bước sau cho đến khi một mẫu được trả về. • Bước 1: Lấy mẫu Y = y từ g(y) và U = u từ phân phối đều U(0, 1). Sang bước 2. f (y) • Bước 2: Nếu u ≤ M g(y) thì đặt X = y . Ngược lại, quay lại bước 1. Mệnh đề 1.10. Phân phối của biến ngẫu nhiên X được lấy mẫu trong phương pháp loại bỏ như trên có mật độ f (x). Thật vây, ta có f (Y ) P(X ≤ x) = P Y ≤ x|U ≤ M g(Y ) f (Y ) P Y ≤ x, U ≤ M g(Y ) = . f (Y ) P U ≤ M g(Y ) Để tính được xác suất trên, ta cần biết mật độ chung của Y và U . Bởi tính độc lập nên: h(y, u) = g(y)1[0≤u≤1] . Vì vậy: Z x Z f (y)/M g(y) f (Y ) P Y ≤ x, U ≤ = g(y) 1dudy M g(Y ) −∞ 0 Z x Z x f (y) 1 = g(y) dy = f (y)dy −∞ M g(y) M −∞ và ∞ f (Y ) 1 1 Z P U≤ = f (y)dy = . M g(y) M −∞ M Dẫn đến: f (Y ) P Y ≤ x, U ≤ M g(Y ) Z x P(X ≤ x) = = f (y)dy. f (Y ) P U ≤ M g(Y ) −∞ Có bao nhiêu lần lặp trong thuật toán chúng ta dùng đến? Trong mỗi lần 10
lặp, chúng ta tạo ra một mẫu với xác suất P(U ≤ Mf g(Y (Y ) 1 ) ) = M nên tổng số lần lặp tuân theo phân phối hình học với tham số 1/M . Do vậy trung bình cần số lần lặp là M . Chú ý sau đây: 1. Cận M nhỏ hơn thì thuật toán hiệu quả hơn trong tổng số lần lặp. Vì vậy chúng ta nên tìm kiếm một mật độ g gần f . 2. Nếu giá của f không bị chặn thì để có thể tìm thấy cận M , mật độ g cần có đuôi lớn hơn f . Ví dụ 1.4. Giả sử chúng ta muốn lấy mẫu |X| trong đó X là biến ngẫu nhiên phân phối chuẩn tắc. Mật độ của |X| được cho bởi r 2 2 x f (x) = exp − với x ∈ R+ . π 2 Ta đã biết cách lấy mẫu một biến ngẫu nhiên phân phối mũ vì thế chúng ta chọn mật độ g là mật độ của một phân phối mũ với tham số 1. Khi đó: r 2 r x − 2x (x − 1)2 f (x) 2 2e = exp − = exp − g(x) π 2 π 2 r 2e ≤ . π q Từ đó, đặt M = 2e π dẫn đến (x − 1)2 f (x) = exp − . M g(x) 2 Thuật toán lấy mẫu loại bỏ tiến hành như sau: • Bước 1: Lấy mẫu Y = y từ phân phối mũ E(1) và U = u từ phân phối đều U (0, 1). Đến bước 2. 2 • Bước 2: Nếu u ≤ exp − (y−1) 2 thì đặt X = y . Ngược lại, trở lại bước 1. Ví dụ 1.5. Xét một biến ngẫu nhiên Y với mật độ g(x) được xác định trên không gian trạng thái S . Bây giờ, giả sử A ⊂ S và chúng ta muốn lấy 11
mẫu biến ngẫu nhiên có điều kiện X = (Y |Y ∈ A) với không gian trạng thái A. Trong trường hợp này mẫu loại bỏ có thể hoàn thành bởi lấy mẫu lặp đi lặp lại X cho đến khi mẫu của chúng ta nằm trong A. Cụ thể hơn, X có mật độ f (x) = P(Yg(x) ∈A) với x ∈ A. Do đó f (x) 1 f (x) ≤ =M và = 1[x∈A] với x ∈ S. g(x) P(Y ∈ A) M g(x) Giả sử U có phân phối đều trên khoảng đơn vị. Khi đó  1 nếu Y ∈ A P(U ≤ f (Y )/M g(y)) = 0 nếu Y ∈ /A Vì vậy, trong thuật toán lấy mẫu loại bỏ tiêu chuẩn, chúng ta chấp nhận nếu Y ∈ A và ngược lại, chúng ta loại bỏ. Chúng ta không cần lấy mẫu U để đưa ra quyết định này. Nếu đánh giá mật độ mục tiêu f là tốn kém thì phương pháp loại bỏ có thể dùng máy điện toán ít tốn kém hơn. Nếu thêm cận trên M g(x) trên mật độ mục tiêu f (x) thì chúng ta cũng có thể dễ dàng ước lượng cận dưới h(x). Vì thế gọi là thuật toán lấy mẫu loại bỏ hình bao, tiến hành như sau: 1. Giả sử Y = y từ g(y) và U = u từ phần phối đều U (0, 1). 2. Chấp nhận nếu u ≤ h(y)/M g(y) và đặt X = y là một mẫu. Ngược lại, đi đến bước 3. 3. Chấp nhận nếu u ≤ f (y)/M g(y) và trả lại X = y là một mẫu. Ngược lại đi đến bước 1. R Điều này hiệu quả hơn vì trung bình ta cần 1/M h(x)dx lần lặp đánh giá của f được thay thế bởi đánh giá của h. Hàm h có thể được tìm thấy trong ví dụ bởi khai triển Taylor. 12
1.3.3 Phương pháp lấy mẫu quan trọng Trong đoạn trước ta đã đưa ra lấy mẫu loại bỏ, sử dụng mật độ đề xuất để tạo ra mẫu từ mật độ mục tiêu. Trong đoạn này, chúng ta vấn tiếp tục lấy mẫu của mật độ mục tiêu nhưng thay đổi cách đánh giá tạo ra ước lượng không chệch của các đặc tính của mật độ mục tiêu. Nhắc lại cái mà ta đang quan tâm khi thảo luận về phương pháp Monte Carlo là tích phân Z I = Ef (h(X)) = h(x)f (x)dx S với f là một mật độ. Khi đó, ta viết lại tích phân dưới dạng f (x) Z I= h(x)g(x)dx S g(x) trong đó, g là một mật độ sao cho g(x) > 0 với f (x)h(x) 6= 0. Bây giờ, chúng ta tạo ra một mẫu độc lập cùng phân phối (x1 , ..., xn ) từ g và ước lượng I bởi: n n ˆ 1 X f (xi ) 1X I= h(xi ) = w(xi )h(xi ) n i=1 g(xi ) n i=1 Ta gọi cách lấy mẫu này là lấy mẫu quan trọng. Mật độ g được gọi là mật độ đề xuất hoặc mật độ công cụ và trọng số w(xi ) = fg(x (xi ) i) được gọi là trọng số quan trọng. Chú ý rằng Iˆ là một ước lượng không chệch của I . Có hai lý do tại sao chúng ta quan tâm đến biểu diễn mẫu quan trọng: 1. Lấy mẫu từ f (x) là không thể hoặc quá đắt đỏ. 2. h(x), trong đó X ∼ f , có phương sai lớn, vì thế ước lượng không chệch theo quy ước có sai số Monte Carlo (MC) lớn. Phương sai của một ước lượng quan trọng sẽ chỉ hữu hạn nếu ước lượng là bình phương khả tích, tức là f 2 (X) 2 2 f (X) Eg h (X) 2 = Ef h (X) < ∞. g (X) g(X) 13
Do đó, phương sai sẽ thường vô hạn nếu tỷ số f (x)/g(x) không bị chặn. Dẫn đến, nếu có thể, chúng ta nên chọn mật độ đề xuất g có đuôi dày hơn f . Tóm lại, nếu f (x)/g(x) không bị chặn thì thậm chí nếu phương sai của ước lượng thống kê là hữu hạn, thủ tục lấy mẫu là không hiệu quả cũng như phương sai của trọng số quan trọng là lớn. Thay vì ước lượng quan trọng Iˆ = n1 ni=1 w(xi )h(xi ), ước lượng tỷ lệ P sau đây thường được sử dụng Pn j=1 h(xj )w(xj ) I˜ = Pn . j=1 w(x j ) Ước lượng này có hai lợi thế: 1. Nó là ước lượng không chệch, thường có phương sai nhỏ hơn ước lượng quan trọng, đưa vào dễ dàng hơn. Nhưng chú ý rằng ước lượng này vẫn phù hợp đối với x1 , ..., xn độc lập cùng phân phối với mật độ g , ta có n 1X n→∞ f (xj )/g(xj ) −−−−→ 1. n j=1 2. Chúng ta có thể áp dụng lấy mẫu quan trọng ngay cả khi chúng ta biết f (x) và vì thế w(x) chỉ đến một hằng số tỷ lệ. Nếu ta không thể tìm thấy một mật độ quan trọng dẫn đến phương sai nhỏ hợp lý của trọng số quan trọng thì có vài phương pháp lấy mẫu có thể áp dụng để làm giảm phương sai: 1. Phép tính gần đúng đầu tiên được gọi là lấy lại mẫu quan trọng liên tiếp và quá trình này như sau: (a) Lấy một mẫu quan trọng Y (1) , ..., Y (n) với các trọng số quan trọng wi = f (Y (i) )/g(Y (i) ), i = 1, ..., n. (b) Tạo một mẫu mới X (1) , ..., X (n) bằng cách lấy mẫu từ Y (1) , ..., Y (n) trong đó Y j được lấy mẫu với xác suất wj / ni=1 wi . P 2. Phương pháp lấy mẫu thứ hai được gọi là kiểm soát loại bỏ và xem xét loại bỏ bất kỳ điểm mẫu mà có trọng số quan trọng dưới một ngưỡng 14
c cho trước. Loại bỏ những điểm mẫu sẽ đưa ra một độ lệch, nhưng bằng sự thay đổi các trọng số quan trọng thích hợp, độ lệch này có thể tránh được. Cho mẫu quan trọng Y (1) , ..., Y (n) với các trọng số quan trọng w1 , ..., wn , quá trình kiểm soát loại bỏ như sau: (a) Với j = 1, ..., n chấp nhận Y (j) với xác suất pj = min{1, wj /c}. Ngược lại, loại bỏ Y (j) . (b) Nếu Y (j) được chấp nhận tính toán lại thì trọng số quan trọng là R w˜j = qwj /pj , trong đó q = min{1, w(x)/c}g(x)dx. Chú ý vì q như nhau đối với tất cả các điểm mẫu nên ta không cần tính nó rõ ràng nếu ta sử dụng ước lượng tỷ lệ. Hơn nữa, kiểm soát loại bỏ tạo ra một mẫu quan trọng theo mật độ đề xuất min{g(x), f (x)/c} g∗ = . q 1.4 Xích Markov Trong đoạn này, chúng ta đưa ra một số định lý về xích Markov quan trọng cho phương pháp MCMC. Định nghĩa 1.11. Xích Markov. Một dãy đại lượng ngẫu nhiên X = {Xn , n = 0, 1, 2, 3, ...} nhận các giá trị trên tập S được gọi là xích Markov nếu: P(Xn+1 ∈ A|Xn = xn ,Xn−1 = xn−1 , ..., X0 = x0 ) = P(Xn+1 ∈ A|Xn = xn ) với mọi n > 0, A ⊆ S , x0 , x1 , ..., xn ∈ S . Đôi khi tính Markov của xích còn được phát biểu dưới dạng: Nếu biết trạng thái hiện tại (tại thời điểm n) của xích thì quá khứ và tương lai (tại thời điểm n+1) độc lập với nhau. 15
Ví dụ 1.6. Giả sử Xn là thời tiết ngày thứ n. Ta đặt:    0 nếu trời nắng vào ngày thứ n  Xn = 1 nếu trời có mây vào ngày thứ n   2  nếu trời mưa vào ngày thứ n Hình sau chỉ ra các xác suất chuyển cho sự thay đổi thời tiết. Bằng việc lấy mô hình thời tiết như xích Markov, chúng ta giả sử rằng Hình 1.1: Xác suất chuyển của xích thời tiết thời tiết ngày mai được tính theo thời tiết hôm nay, không phụ thuộc vào ngày hôm qua hay bất kỳ ngày trước nào. Định nghĩa 1.12. Xác suất chuyển, Xích thời gian thuần nhất. Một xích Markov X được gọi là xích thuần nhất nếu xác suất chuyển của nó: Z P(Xn+1 ∈ A|Xn = x) = P (x, A) = p(x, y)dy A không phụ thuộc vào n. Ta gọi P(x, A) là nhân chuyển. Trong phạm vi ở đây, chúng ta giả sử rằng nhân chuyển là liên tục tuyệt đối với mọi x ∈ S , tức là nó có một mật độ liên quan hoặc hàm khối xác suất. Vì vấy, cố định x ∈ S , hàm p(x, y) là một mật độ hoặc hàm khối xác suất (pmf). Xác suất chuyển sau n bước của X được định nghĩa bởi Z P(Xn ∈ A|X0 = x) = P (n) (x, A) = p(n) (x, y)dy. A 16