Tài liệu Phương pháp MCMC và vài ứng dụng

ĐẠ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 LUẬN VĂN THẠC SỸ KHOA 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 NGƯỜI HƯỚNG DẪN KHOA HỌC: TS. NGUYỄN THỊNH HÀ NỘI, 2014 Mục lục LỜI MỞ ĐẦU 5 BẢNG KÝ HIỆU 7 1 TỔNG QUAN 1.1 Suy luận Bayes . . . . . . . . . . . . . . . . . . . 1.1.1 Đặc điểm mô hình Bayes . . . . . . . . . . 1.1.2 Các tiên nghiệm Jeffreys . . . . . . . . . . 1.2 Tích phân Monte Carlo . . . . . . . . . . . . . . . 1.2.1 Bài toán . . . . . . . . . . . . . . . . . . . 1.2.2 Xấp xỉ Monte Carlo . . . . . . . . . . . . . 1.2.3 Monte Carlo thông qua lấy mẫu theo trọng 1.3 Phương pháp sinh biến ngẫu nhiên . . . . . . . . . 1.3.1 Phương pháp biến đổi . . . . . . . . . . . . 1.3.2 Phương pháp chấp nhận - bác bỏ . . . . . . 1.3.3 Phương pháp tỷ số đều . . . . . . . . . . . 1.4 Xích Markov . . . . . . . . . . . . . . . . . . . . . 1.4.1 Các định nghĩa và kí hiệu . . . . . . . . . . 1.4.2 Sự hội tụ của phân phối . . . . . . . . . . . 1.4.3 Giới hạn của giá trị trung bình . . . . . . . . . . . . . . . . . . . số . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 8 9 9 10 10 11 12 13 13 14 15 16 18 19 19 2 MẪU GIBBS 21 2.1 Mẫu Gibbs . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 2.2 Thuật toán mở rộng dữ liệu . . . . . . . . . . . . . . . . . . 24 3 THUẬT TOÁN METROPOLIS-HASTINGS 27 3.1 Thuật toán Metropolis – Hastings . . . . . . . . . . . . . . 27 3.1.1 Khái niệm . . . . . . . . . . . . . . . . . . . . . . . 27 2 3.2 3.3 3.4 3.1.2 Mẫu độc lập . . . . . . . . . . . . . . . . . . . . . 3.1.3 Xích bước ngẫu nhiên . . . . . . . . . . . . . . . . Thuật toán Metropolis- Hasting cho các phân phối nhiều chiều . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.1 Cập nhật từng khối . . . . . . . . . . . . . . . . . 3.2.2 Cập nhật từng thành phần . . . . . . . . . . . . . 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 . . . . . . . . . . . . . . 3.3.2 Thuật toán Langevin . . . . . . . . . . . . . . . . 3.3.3 Thuật toán đa phép thử MH . . . . . . . . . . . . Thuật toán bước nhảy ngược MCMC cho bài toán lựa chọn mô hình Bayes . . . . . . . . . . . . . . . . . . . . . . . . 3.4.1 Thuật toán bước nhảy ngược MCMC . . . . . . . 3.4.2 Xác định điểm thay đổi . . . . . . . . . . . . . . 4 Phương pháp biến phụ trợ MCMC 4.1 Mô phỏng nhiệt luyện . . . . . . . 4.2 Mô phỏng điều hoà nhiệt . . . . . 4.3 Thuật toán Moller . . . . . . . . . 4.4 Thuật toán trao đổi . . . . . . . . Tài liệu tham khảo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 . 30 . . . . . . . 30 30 34 36 36 37 38 . 39 . 39 . 43 . . . . 46 48 49 51 53 56 3 LỜI CẢM ƠN Luận văn này được hoàn thành với sự hướng dẫn tận tình và cũng hết sức nghiêm khắc của TS. Nguyễn Thịnh. Thầy đã dành nhiều thời gian quý báu của mình để hướng dẫn cũng như giải đáp các thắc mắc của tôi trong suốt cả quá trình làm luận văn. Tôi muốn tỏ lòng biết ơn chân thành và sâu sắc nhất tới người thầy của mình. Tôi cũng muốn gửi tới toàn thể các thầy cô Khoa Toán - Cơ - Tin học trường Đại học Khoa học Tự nhiên, Đại học Quốc gia Hà Nội, các thầy cô đã đảm nhận giảng dạy khóa Cao học 2011 - 2013, đặc biệt là các thầy cô tham gia giảng dạy nhóm Xác suất thống kê 2011 - 2013 lời cám ơn chân thành đối với công lao dạy dỗ trong suốt thời gian của khóa học. Tôi xin cám ơn gia đình, bạn bè, đồng nghiệp và các anh chị em trong nhóm Xác suất thống kê 2011 - 2013 đã quan tâm, giúp đỡ, tạo điều kiện và động viên tinh thần để tôi có thể hoàn thành được khóa học này. 4 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. Thuật toán Metropolis - Hasting đối với phân phối nhiều chiều: giới thiệu ứng dụng của thuật toán Metropolis - Hasting đối với các biến ngẫu nhiên nhiều chiều bằng cập nhật từng khối, cập nhật từng thành phần. Các dạng khác nhau của thuật toán Metropolis - Hasting: Thuật toán chạm và chạy, thuật toán Langevin, thuật toán đa phép thử MH. Chương 4. Phương pháp biến phụ trợ MCMC. 5 Giới thiệu về mặt lý thuyết một vài thuật toán của phương pháp MCMC có sử dụng các biến phụ trợ: Phương pháp mô phỏng nhiệt luyện, mô phỏng điều chỉnh nhiệt,Moller, thuật toán trao đổi, phương pháp lấy mẫu MH kép. Do thời gian gấp rút và kiến thức còn hạn chế nên luận văn không thể tránh khỏi những thiếu sót, vì vậy, rất mong nhận được những ý kiến đóng góp của các thầy cô và bạn bè đồng nghiệp, xin trân trọng cám ơn. Hà Nội, tháng 11 năm 2014 6 BẢNG KÝ HIỆU MCMC: Xích Markov Monte Carlo AD: Thuật toán mở rộng dữ liệu AR: Thuật toán chấp nhận - bác bỏ h.c.c: hầu chắc chắn MTH: thuật toán đa phép thử Metropolis - Hastings MTM: thuật toán đa phép thử Metropolis RJMCMC: Thuật toán bước nhảy ngược MCMC 7 Chương 1 TỔNG QUAN 1.1 Suy luận Bayes Suy luận Bayes là một công thức suy luận xác suất. Với ưu điểm là tính toán đơn giản và cùng với những phát triển gần đây của các phương pháp xích Markov Monte Carlo(MCMC) cho việc tính xấp xỉ tích phân có số chiều cao mà suy luận Bayes ngày càng được sử dụng rộng rãi. Suy luận Bayes được bắt nguồn từ Thomas Bayes (1764), người đã rút ra xác suất nghịch đảo của xác suất thành công θ trong một dãy các phép thử độc lập Bernoulli, trong đó θ được lấy từ phân phối đều trên khoảng (0,1). Ví dụ 1.1. (Mô hình Bernoulli với tiên nghiệm đã biết) Giả sử rằng θ ∼ U nif (0, 1) là phân phối đều trên khoảng (0,1),và x1 , x2 , ..., xn 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. Pn Ta viết N = i=1 xi 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 x1 , x2 , ..., xn đượ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 θ(1+N )−1 (1 − θ)(1+n−N )−1 (0 ≤ θ ≤ 1) (1.2) B(1 + N, 1 + n − N ) 8 trong đó B (◦ ,◦ ) là kí hiệu của hàm Beta 1.1.1 Đặc điểm mô hình Bayes Theo như những nghiên cứu toán học đã biết thì để xác định mô hình Bayes ta cần : (i) Chỉ rõ một mô hình lấy mẫu từ dữ liệu quan sát X, có điều kiện trên một đại lượng chưa biết θ. (X ∈ X , θ ∈ Θ) X ∼ f (X |θ ) (1.3) ở đó f (X |θ ) là hàm mật độ xác suất, và (ii) Chỉ rõ một phân phối biên,được gọi là phân phối tiên nghiệm hay đơn giản là tiên nghiệm π (θ) của θ: θ ∼ π (θ) (θ ∈ Θ) (1.4) Phân tích dữ liệu dựa trên kết quả những suy luận ở trên nhằm mục đích rút gọn tính toán tích phân đối với phân phối hậu nghiệm, hay nói gọn là hậu nghiệm, π (θ |X ) = R π (θ) L (θ |X ) π (θ) L (θ |X ) dθ (θ ∈ Θ) (1.5) ở đó 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 Một cách tự nhiên ta thấy rằng việc chỉ rõ mô hình Bayes chẳng khác gì việc tổng hợp các thông tin có thể trong thực tế theo quan điểm xác suất chính xác. Đồng thời, việc chỉ rõ mô hình xác suất đối với dữ liệu quan sát X là việc làm tất yếu. Thêm vào đó khi xét mô hình lấy mẫu của dữ liệu quan sát X đối với đại lượng chưa biết θ suy luận Bayes yêu cầu tiên nghiệm cho θ phải được xác định rõ ràng. Trong trường hợp thông tin tiên 9 nghiệm của θ là sẵn có và có thể biết một cách chính xác bởi một phân phối xác suất thì điều này là hiển nhiên. Tuy nhiên, đố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: 1 πJ (θ) ∝ |I (θ)| 2 (θ ∈ Θ) (1.6) Trong đó I (θ) là lượng thông tin Fisher. Ví dụ 1.2. Giả sử rằng ta xét một mẫu được lấy từ phân phối N (µ, 1) Thông tin Fisher thu được như sau: Z+∞ I (µ) = φ (x − µ) dx = 1 −∞ Trong đó   1 2 φ (x − µ) = (2π) exp − (x − µ) 2 1 2 là hàm mật độ của N (µ, 1). Điều này dẫn đến tiền nghiệm Jeffreys của θ là πJ (θ) ∝ 1 (−∞ < µ < +∞) (1.7) Ta thu được 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.2 1.2.1 (1.8) Tích phân Monte Carlo 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 ⊆ Rd , trong đó Rd 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: Z Eν [h (X)] = h (x) ν (dx) (1.9) X 10 Trong đó h(x) là hàm đo được. Giả sử rằng ν có hàm mật độ xác suất f (x) thì (1.9) có thể được viết thành: Z Ef [h (X)] = h (x) f (x) dx (1.10) X Ví dụ 1.3. Để ước lượng xác suất P r (X ∈ S) với S ∈ X , h (x) hàm chỉ tiêu là: h (x) = Ix∈S với  1, nếu x ∈ S h (x) = 0, nếu ngược lại , và tính toán phân phối thành phần fY (y) từ phân phối đồng thời fX,Y (x, y).   Khi đó thay vào trong (1.10) ta được là EfX fY |X (y|x) , trong đó fX (x) là hàm mật độ của thành phần X,và fY |X (y |x) là hàm mật độ có điều kiện của Y đối với X đã biết. 1.2.2 Xấp xỉ Monte Carlo Ta kí hiệu X1 , ..., Xn là một mẫu kích thước n lấy từ hàm mật độ xác suất f (x) trong (1.10). Khi đó trung bình mẫu của h (X) là: n 1X hn = h (Xi ) n i=1 (1.11) có thể được sử dụng để tính xấp xỉ (1.10) vì hn 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à: hn − Ef [h (X)] p ∼ N (0, 1) nV ar (h (X)) Tương tự V ar (h (X)) có thể được xấp xỉ bằng phương sai mẫu: n
2 1 X h (X1 ) − hn n − 1 i=1 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 11 1.2.3 Monte Carlo thông qua lấy mẫu theo trọng số Trong trường hợp ta gặp khó khăn khi sinh trực tiếp các mẫu từ f (x), ta có thể sử dụng phương pháp lấy mẫu theo trọng số, phương pháp này dựa trên phép đồng nhất sau đây: R R (x) Ef [h (X)] = h (x) f (x) dx = h (x) fg(x) g (x) dx X X = Eg [h (X) f (X) /g (X)] Trong đó g (x) là hàm mật độ xác suất trên X và g(x) > 0 với mọi x mà tại đó f (x) > 0. Đồng nhất thức này cho chỉ ra rằng các mẫu có các hàm mật độ khác nhau xuất phát từ f (x) cũng có thể được xấp xỉ (1.10). Lý thuyết Monte Carlo áp dụng được trong trường hợp này vì:   h i f (X) e Ef [h (X)] = Eg h (X) = Eg h (X) g (X) trong đó f (x) e h (x) = h (x) g (x) g (x) Ước lượng của Ef [h (X)] bây giờ trở thành: m 1 X f (x1 ) h= h (xi ) m i=1 g (xi ) (1.12) trong đó x1 , ..., xn là các mẫu độc lập cùng phần phối sinh ra từ g (x). So (xi ) sánh với (1.11), với mỗi i = 1, ..., m xi có trọng số wi = fg(x . Chính vì lý i) do đó mà phương pháp này được gọi là phương pháp lấy mẫu theo trọng số. Vấn đề mấu chốt của phương pháp này là chọn g (x) thỏa mãn cả tính đơn giản trong việc sinh ra các mẫu Monte Carlo và độ chính xác trong ước lượng Ef [h (X)] bằng cách kiểm soát các sai số Monte Carlo. Với độ tin cậy Monte Carlo,ta cần chọn g (x) để cực tiểu phương sai của e h (X) với X ∼ g (x). Người ta chứng minh được rằng hàm g(x) thoả mãn điều kiện trên là: |h (x)| f (x) g ∗ (x) = R |h (y)| f (y) dy X 12 1.3 Phương pháp sinh biến ngẫu nhiên Phương pháp MC dựa trên việc lấy mẫu từ các phân phối xác suất. Mặt khác,dựa vào phân phối đều U nif (0, 1) ta có thể sinh được các số ngẫu nhiên của một phân phối xác suất bất kỳ .Do đó phương pháp sinh một mẫu độc lập cùng phân phối từ phân phối đều đơn giản nhất U nif (0, 1) là rất quan trọng bởi vì toàn bộ các phương pháp lấy mẫu đều dựa trên các số ngẫu nhiên đều được sinh ra. 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. Tuy nhiên,thuật toán này nhìn chung có tính toán phức tạp. Các phương pháp hiệu quả hơn sẽ được mô tả trong phần sau của luận văn. Sau đây là một số phương pháp thường được sử dụng để lấy mẫu từ các các phân phối trong trường hợp công thức hàm phân bố ngược không thể áp dụng được. 1.3.1 Phương pháp biến đổi Phương pháp biến đổi dựa trên phép biến đổi của các biến ngẫu nhiên,thuật toán 1.1 và 1.2 là một ví dụ. Tuy nhiên,ngoại trừ một vài trường hợp như phân phối mũ và phân phối Bernoulli thì thuật toán 1.1 và 1.2 thường không hiệu quả. 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). Sau đây là một số ví dụ thường được sử dụng trong thực hành. 13 Công thức Phép biến đổi Mũ X = −ln(U ) Cauchy X = tan (πU − π/2)) Beta 1.3.2 ind Xi ∼ Gamma (αi ) , i = 1, 2 Phân phối X ∼ Expo(1) X ∼ Cauchy(0, 1) X1 X1 +X2 ∼ Beta (α1 , α2 ) Phương pháp chấp nhận - bác bỏ Phương pháp chấp nhận - bác bỏ (AR) rất hữu ích trong việc sinh các số ngẫu nhiên khi các phương pháp biến đổi trực tiếp không tồn tại hoặc tính toán không hiệu quả. Ta mô tả phương pháp AR thông qua một đối số hình học. Xét mẫu có phân phối d - chiều với không gian mẫu X ⊆ Rd . Theo định nghĩa về hàm mật độ, miền phía dưới đường cong/mặt phẳng của hàm mật độ Cf = {(x, u) : 0 ≤ u ≤ f (x)} ⊂ Rd+1 (1.13) bằng một đơn vị thể tích.Do đó nếu (X,U) là đều trong miền Cf thì X ∼ f (x). Chú ý rằng X ∼ f (x) vẫn đúng khi f (x) trong (1.13) được làm bội bởi một hằng số dương tùy ý, nghĩa là: Ch = {(x, y) : 0 ≤ u ≤ h (x)} ⊂ Rd+1 (1.14) trong đó h (x) ∝ f (x),bởi sự thay đổi tỷ lệ trên U sẽ không ảnh hưởng đến phân phối biên của X. Điều này có nghĩa là ta có thể sinh ra X bằng các điểm mô phỏng phân phối đều trên Cf hoặc Ch . Khi ta gặp khó khăn để lấy mẫu một cách trực tiếp từ Ch ,ta có thể lấy mẫu một cách gián tiếp qua Ch 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 ⊇ Ch và (ii) Thu thập những điểm thuộc vào miền Ch . 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 (x) đơn giản với hàm mật độ g (x) thoả mãn fg(x) bị chặn trên bởi một số hằng số hữu hạn M. Vì vậy Ch là đóng trong miền: Cg = {(x, u) : 0 ≤ u ≤ g (x)} ⊂ Rd+1 14 (1.15) 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 x h (x) ≤M <∞ g (x) Thuật toán 1.3(Chấp nhận -bác bỏ) Lặp lại 2 bước sau cho đến khi một giá trị được trả về trong bước 2: 1, Sinh ra X từ g(x) và U từ U nif (0, 1). 2, Nếu U ≤ f (X) M g(X) , trả lại giá trị X (như là độ lệch ngẫu nhiên từ f (x)). Trong trường hợp hàm số h(x) khó ước lượng, ta sử dụng hàm số kẹp s(x) 0 ≤ s (x) ≤ h (x) có tính toán đơn giản hơn để rút gọn việc tính toán h(x). Thuật toán 1.4 (Chấp nhận - bác bỏ với hàm số kẹp). Lặp lại hai bước sau đây cho đến khi một giá trị xuất ra trong bước 2: 1, Sinh ra X từ g (x) và U từ U nif (0, 1). S(X) 2, Nếu U ≤ Ms(X) g(X) hoặc M g(X) < U ≤ một độ lệch ngẫu nhiên từ f (x)). Do đó trong trường hợp này U ≤ 1.3.3 s(X) M g(X) , h(X) M g(X) trả lại giá trị X (như là thuật toán không ước lượng h (x) Phương pháp tỷ số đều Phương pháp tỷ số đều là phương pháp thông dụng để sinh các số ngẫu nhiên của nhiều phân phối thông dụng như phân phối Gamma, chuẩn, và student-t. Ý tưởng tổng quát của phương pháp tỷ số đều là tìm ra một cặp phép biến đổi khả vi U = u(Y ) và X = x(Z, Y ) với U = u(Y ) tăng thực sự để 15 thoả mãn (1.14) và do đó với một hằng số Jacobi thì (Y, Z) cũng đều trên tập ảnh tương ứng của Ch :  (Y,Z) Ch = (y, z) : u−1 (0) ≤ y = u−1 (u) ≤ u−1 (h (x (z, y))) ⊂ Rd+1 (1.16) trong đó u−1 (.) là hàm số ngược của u(.). Điều này dẫn tới thuật toán bác bỏ tổng quát như sau: 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: (Y,Z) 1, Sinh (Y, Z) có độ lệch đều trên miền D ⊇ Ch (Y,Z) 2, Nếu (Y, Z) ∈ Ch muốn. , trả về giá trị X = x(Y, Z) là độ lệch mong Thuật toán này có tỉ số chấp nhận R R (Y,Z) dydz dx C X h (x) R r = Rh = Jx,u (z, y)  D dydz D dydz Trong đó   ∂x  Jx,u (z, y) =  ∂z  0 ∂x ∂y ∂u ∂y      ∂x   0   = u (y)   ∂z   là hệ số Jacobi của các phép biến đổi. 1.4 Xích Markov Trong trường hợp việc sinh các mẫu độc lập cùng phân phối từ phân phối mục tiêu π là không thể thực hiện được, các mẫu phụ thuộc {Xi } có thể sử dụng thay thế, với điều kiện là trung bình mẫu (1.11) hội tụ tới (1.10). Xích Markov là một dãy các biến ngẫu nhiên {Xi , 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 (Xt+1 ∈ A |X0 = x0 , ..., Xt = xt ) = Pr (Xt+1 ∈ A |Xt = xt ) 16 (1.17) với thời gian t = 0, 1, ....Để thuận lợi cho việc trình bày ta sử dụng kí hiệu π (dy) để chỉ độ đo xác suất π trên (X, X ) cho cả trường hợp biến ngẫu nhiên rời rạc và liên tục. Với biến liên tục X, hàm mật độ f (x) của nó là đạo hàm Radon - Nikodym của độ đo xác suất π (dx) đối với độ đo Lebesgue. Đối vớibiến ngẫu nhiên rời rạc X, hàm mật độ f (x), là đạo hàm của π (dx) đối với độ đo đếm. Do vậy, kí hiệu Pt (dx) cho phân phối của Xt đối với trạng thái X tại thời điểm t. Xuất phát với phân phối ban đầu P0 (dx), xích Markov {Xt } khai triển như sau: Pt+1 (dy) = ∫ Pt (dx) Pt (x, dy) (1.18) X Phân phối Pt (x, dy) là độ đo xác suất đối với Xt+1 trong đó Xt = x cho trước và được gọi là phân phối hạch chuyển dịch tại thời điểm t. Trong thực tế, đây là hàm mật độ có điều kiện của Xt+1 với Xt = x cho trước. Một lớp các xích Markov sơ cấp thường được sử dụng trong MCMC là các xích Markov dừng, trong đó Pt (x, dy) = P (x, dy) với t = 1, 2, .... Trong trường hợp này (1.18) trở thành: Z Pt+1 (dy) = Pt (dx) P (x, dy) (1.19) (1.20) X và Pt (dx) là xác định duy nhất bởi phân phối ban đầu P0 (dx) và hạch chuyển dịch P (x, dy). Tổng quát lên ta kí hiệu P n (x, .) là phân phối có điều kiện của Xt0 +n với Xt0 = x cho trước. Ý tưởng cơ bản cho việc tạo ra các xích Markov để xấp xỉ Eπ (h (X)) là xây dựng một hạch chuyển dịch P (x, dy) với π (dx) là phân phối dừng, nghĩa là P (x, dy) và π (dx) thỏa mãn điều kiện cân bằng: Z π (dy) = π (dx) P (x, dy) (1.21) X Trong trường hợp phân phối mục tiêu π có hàm mật độ xác suất f (x) và hạch dịch chuyển P (x, dy) có hàm mật độ điều kiện p (y |x), điều kiện cân bằng có thể được viết dưới dạng: Z f (y) = p (y |x) f (x) dx χ 17 Chú ý: Nếu với hầu hết π(x) và với mọi tập đo được A ta có: lim Pr (Xt ∈ A |X0 = x) = π (A) thì π(dx) được gọi là phân phối cân t→∞ bằng của xích Markov. 1.4.1 Các định nghĩa và kí hiệu Định nghĩa 1.1. Cho Xn là một xích bất khả quy với phân phối dừng π (.) và kí hiệu {An i.o} là một dãy xuất hiện thường xuyên vô hạn, nghĩa P là i IAi = ∞ với xác suất 1 (a) Xích là hồi quy nếu với mọi B thoã mãn π (B) > 0,thì Pr (Xn ∈ Bi.o. |X0 = x) > 0 với mọi x và P r (Xn ∈ Bi.o. |X0 = x) = 1 với hầu hết π (x) (b) Xích là hồi quy Harris nếu P r (Xn ∈ Bi.o. |X0 = x) = 1 với hầu hết π(x) Để xác định các dạng khác của ergodic, ta sử dụng khái niệm tổng biến thiên khoảng cách giữa hai độ đo trên X và khái niệm thời điểm chạm. Tổng biến thiên khoảng cách giữa hai độ đo trên (X, X ) xác định bằng tổng biến thiên chuẩn của độ đo λ trên (X, X ) kλk = sup λ (A) − inf λ (A) A∈X A∈X (1.22) Thời điểm chạm của tập con B ∈ X là biến ngẫu nhiên: HB = inf {t ≥ 0 : Xt ∈ B} Trong đó cận dưới đúng của tập rỗng tiến tới ∞ Đị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 HB 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: Z
Ex HB2 π (dx) < ∞ B 18 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: kP n (x, .) − πk ≤ M (x) rn ∀x (d) Xích trong (c) được gọi là ergodic đều nếu tồn tại một hằng số M và một hằng số dương r < 1 sao cho kP n (x, .) − πk ≤ M rn 1.4.2 Sự hội tụ của phân phối Tổng biến thiên khoảng cách giữa hai độ đo trên (X, X ) đã được sử dụng để mô tả sự hội tụ của một xích Markov trong định lý sau đây (Định lý 1 của Tierney, 1994) Định lý 1.1. Giả sử rằng P (x, dy) có π(x) là bất khả quy và dừng. Khi đó P (x, dy) là hồi quy dương và π (dx) là phân phối dừng duy nhất của P (x, dy). Nếu P (x, dy) cũng không tuần hoàn thì với hầu hết π (x): kP n (x, .) − πk → 0 với k.k là tổng biến thiên khoảng cách. Nếu P (x, dy) là hồi quy Harris thì nó hội tụ với mọi x 1.4.3 Giới hạn của giá trị trung bình Định lý 1.2. Giả sử rằng Xn là ergodic với phân phối cân bằng f (x) và giả sử h (x) có giá trị thực và Ef (|h (X)|) < ∞. Khi đó với bất kỳ phân phối ban đầu, hn → Ef (h (X)) h.c.c . Định lý 1.3. Giả sử rằng Xn 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 hn − Ef (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 σh2 với mọi phân phối ban đầu. 19 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à Ef h2 (X) < ∞ Định lý 1.4. Giả sử rằng Xn là ergodic đều với phân phối cân bằng f (x)
và giả sử h (x) có giá trị thực và Ef h2 (X) < ∞. Khi đó tồn tại một số
√ thực σh sao cho phân phối của n hn − Ef (h (X)) hội tụ yếu tới phân phối chuẩn với kỳ vọng 0 và phương sai σh2 với mọi phân phối ban đầu. 20