Tài liệu Thuật toán Metropolis 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 THỊ PHI DOAN THUẬT TOÁN METROPOLIS VÀ ỨNG DỤNG LUẬN VĂN THẠC SỸ Chuyên ngành : LÝ THUYẾT XÁC SUẤT VÀ THỐNG KÊ TOÁN HỌC Mã số : 60 460 106 Giáo viên hướng dẫn: TS. TRẦN MẠNH CƯỜNG HÀ NỘI, 2014 Mục lục 1 2 LỜI MỞ ĐẦU 4 BẢNG KÝ HIỆU 6 KIẾN THỨC CHUẨN BỊ 1.1 Các phương pháp mô phỏng biến ngẫu nhiên . . . . 1.1.1 Phương pháp lấy mẫu ngược . . . . . . . . . 1.1.2 Phương pháp lấy mẫu loại trừ . . . . . . . . 1.2 Ước lượng bằng mô phỏng . . . . . . . . . . . . . . 1.2.1 Lấy mẫu quan trọng (Importance Sampling) . 1.3 Xích Markov . . . . . . . . . . . . . . . . . . . . . 1.3.1 Giới thiệu về xích Markov . . . . . . . . . . . 1.3.2 Phân bố dừng . . . . . . . . . . . . . . . . . 1.3.3 Phân bố giới hạn . . . . . . . . . . . . . . . 1.3.4 Xích tối giản và không có chu kì . . . . . . . 1.3.5 Xích khả nghịch . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 7 7 8 10 11 13 14 18 20 24 27 THUẬT TOÁN METROPOLIS-HASTINGS 28 2.1 Giới thiệu về MCMC . . . . . . . . . . . . . . . . . . . . . 28 2.2 Thuật toán Metropolis-Hastings . . . . . . . . . . . . . . . 29 3 ÁP 3.1 3.2 3.3 DỤNG THUẬT TOÁN METROPOLIS 46 Giới thiệu về R . . . . . . . . . . . . . . . . . . . . . . . . 46 Mô hình lõi cứng (hard-core model) . . . . . . . . . . . . . 47 Thuật toán và chương trình . . . . . . . . . . . . . . . . . 48 Tài liệu tham khảo 53 2 LỜI CẢM ƠN Luận văn trên em hoàn thành dưới sự hướng dẫn tận tình và cũng hết sức nghiêm khắc của TS. Trần Mạnh Cường. 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 em trong suốt cả quá trình làm luận văn. Em muốn tỏ lòng biết ơn chân thành và sâu sắc nhất tới Thầy. Em 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 việc giảng dạy khóa Cao học 2011 - 2013 từ lúc chúng em ôn thi đầu vào và trong cả quá trình học tại Khoa, đặ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 nhất vì công lao dạy dỗ chúng em trong suốt thời gian của khóa học. Vì thời gian cũng như kiến thức của em cũng còn nhiều hạn chế nên không tránh khỏi những thiếu sót. Kính mong các Thầy, Cô nhận xét và cho em những lời nhận xét góp ý để bản luận văn của em được hoàn thiện hơn. Và có thể những góp ý quý giá đó sẽ mở hướng cho em trong quá trình học tập, nghiên cứu sau này. 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. 3 LỜI MỞ ĐẦU Đầu thế kỷ XX, nhà vật lý và bác học nổi tiếng người Nga A.A.Markov đã đưa ra một mô hình toán học để mô tả chuyển động của các phân tử chất lỏng trong một bình kín. Sau này mô hình Markov được phát triển và mang tên: Quá trình Markov. Xích Markov là trường hợp riêng của quá trình Markov khi ta có thể đánh số được các trạng thái. Chúng ta đều biết vai trò quan trọng của thuyết Monte Carlo trong việc ước lượng các số nguyên và mô phỏng các quá trình ngẫu nhiên. Bước có tính quyết định nhất trong việc phát triển lý thuyết Monte Carlo hiệu quả là ước lượng (lấy mẫu) từ một phân bố xác suất thích hợp π(x). Ta không thể trực tiếp tạo thành các mẫu độc lập từ π(x), cũng như chọn cách lấy mẫu quan trọng, trong các mẫu ngẫu nhiên được lấy từ một phân bố thử dạng khác (nhưng gần giống) với phân bố mục tiêu và sau đó đánh giá dựa vào tỉ số quan trọng; hoặc xây dựng các mẫu xác suất độc lập dựa trên ý tưởng lấy mẫu Markov Chain Monte Carlo. Cho π(x) = Z −1 exp(−h(x)) là phân bố mục tiêu dựa vào kết quả nghiên cứu (có thể tất cả các hàm phân phối xác suất được viết dạng này), khi mà hằng số chuẩn hóa hay hàm phân bố Z mà chúng ta thường R không biết. Về mặt lý thuyết, Z = exp (−h(x))dx có thể tính nhưng không dễ dàng (và thường khó) hơn vấn đề ban đầu là mô phỏng từ π . Được thúc đẩy bởi các vấn đề tính toán các xác suất vật lý, Metropolis đã giới thiệu ý tưởng cơ bản của việc phát triển một quá trình Markov để đạt được việc lấy mẫu của π . Ý tưởng này sau đó được phát triển thành thuyết Metropolis dù đơn giản nhưng rất hữu ích và hiện nay được sử dụng rộng rãi bởi các nhà nghiên cứu trong rất nhiều lĩnh vực khoa học khác nhau như sinh học, hóa học, khoa học máy tính, kinh tế học, các ngành kỹ thuật, khoa học vật liệu và nhiều lĩnh vực khác. Luận văn gồm có 3 chương: Chương 1: Kiến thức cơ sở- Xích Markov: Ở phần đầu em trình bày 4 các phương pháp cơ bản để mô phỏng biến (mẫu) ngẫu nhiên như phương pháp ngược, phương pháp lấy mẫu quan trọng, phương pháp lấy mẫu loại trừ. Tiếp theo là ước lượng bằng mô phỏng. Phần tiếp theo của chương I là lý thuyết xích Markov. Mỗi phần đều có các ví dụ minh họa để việc tiếp cận vấn đề trở nên dễ dàng hơn. Chương 2: Thuật toán Metropolis-Hastings: Cũng là phần chính của luận văn. Trong chương này, em đề cập tới các kiến thức cơ bản để xây dựng thuật toán Metropolis. Em giới thiệu phương pháp MCMC và nêu thuật toán Metropolis cùng các ví dụ cụ thể áp dụng thuật toán này. Chương 3: Áp dụng thuật toán Metropolis: Trong chương này em giới thiệu về ngôn ngữ lập trình R và các tính năng của nó. Tiếp đó em ứng dụng vào bài toán mô hình lõi cứng (hard-core model) để viết một đoạn chương trình áp dụng ngôn ngữ R cho kết quả cụ thể. Trong chương này em nêu thuật toán và chương trình áp dụng trên máy tính. Hà Nội, tháng 05 năm 2014 Học viên Nguyễn Thị Phi Doan 5 BẢNG KÝ HIỆU MCMC: Markov Chain Monte Carlo : điều phải chứng minh BNN: biến ngẫu nhiên E(X): Kỳ vọng Var(X): Phương sai ⊗ : Kết thúc một ví dụ F(x): hàm phân phối tích lũy f(x): hàm mật độ 6 Chương 1 KIẾN THỨC CHUẨN BỊ 1.1 1.1.1 Các phương pháp mô phỏng biến ngẫu nhiên Phương pháp lấy mẫu ngược Định lý 1.1. Cho hàm phân phối tích lũy F(x) với F −1 là hàm ngược của F được xác định như sau: F −1 (u) = min{x|F (x) ≥ u}. với u ∈ (0, 1]. Cho U là một BNN có phân phối đều U(0,1) và đặt X = F −1 (U ) khi đó hàm phân phối của X là F(x). Ví dụ 1.1 Mô phỏng một biến ngẫu nhiên có phân phối mũ với tham số λ. Hàm phân phối mũ có dạng: F (x) = 1 − exp(−λx) với x ≥ 0. 1 Cho U ∼ U (0, 1) và đặt Y = − log(1 − U ). λ Khi đó Y có phân phối mũ với tham số λ. Hơn thế, nếu U ∼ U (0, 1) thì 1 − U cũng có phân phối U (0, 1) do đó nếu đặt 1 Y = − log(U ) có phân phối mũ với tham số λ. ⊗ λ 7 Ví dụ 1.2 Phân phối Becnoulli và phân phối nhị thức B(n, p). Cho U ∼ U (0, 1). Nếu:  1 khi U < p X= 0 trong các trường hợp khác . Thì X là một phân phối Becnoulli với xác suất thành công là p. Cho n X X1 , ...Xn là các BNN độc lập của phân phối Becnoulli khi đó Y = Xi i=1 là phân phối nhị thức B(n, p). ⊗ Ví dụ 1.3 Mô phỏng biến ngẫu nhiên có phân phối hình học (Geo(p)). Giả sử X nhận giá trị thuộc tập N và P (X = j) = Pj Khi đó F −1 (u) = min{j ∈ N |u ≤ n X pi }. i=1 Nếu X có phân phối hình học X ∼ Geo(p) thì P (X > j) = (1 − p)j điều đó có nghĩa là: n X pi = 1 − (1 − p)j ≥ u ⇔ j > i=1 Kí hiệu [a]: phần nguyên của a, khi đó X = [ học. 1.1.2 log(1 − U ) . log(1 − p) log(U ) ] có phân phối hình log(1 − p) ⊗ Phương pháp lấy mẫu loại trừ Giả sử ta có BNN X với hàm mật độ f (x). Ta chưa mô phỏng được X nhưng ta có thể mô phỏng Y với hàm mật độ cho biết trước là g(y) mà f (x) ≤ M với ∀x. giá của f là tập con giá của g và g(x) Sau đó ta sử dụng mẫu Y để tìm mẫu X. Lặp lại các bước sau tới khi cho kết quả. Bước 1: Mô phỏng Y = y từ g(x) và U = u từ phân phối đều U (0, 1). Chuyển sang bước 2, 8 Bước 2: Nếu u 6 f (y) đặt X = y . Ngược lại quay về bước 1. M.g(y) Mệnh đề 1.1. Biến ngẫu nhiên X được lấy dựa trên phương pháp loại trừ như trên có hàm mật độ f (x). Câu hỏi đặt ra là ta cần bao nhiêu phép lặp trong thuật toán này? Trong f (Y ) 1 mỗi phép lặp ta xây dựng một mẫu với xác suất P (U 6 )= . M.g(Y ) M Vậy trung bình số phép lặp là M. Chú ý 1, M càng nhỏ càng có lợi trong quá trình tính toán vì số bước lặp ít. Vậy ta nên tìm hàm mật độ g gần với hàm f. 2, Nếu giá của hàm f không bị chặn thì để tìm M ta xác định hàm mật độ g có đuôi nặng hơn f. Ví dụ 1.4 Giả sử ta tìm mẫu |X| mà X là BNN có phân phối chuẩn tắc khi đó hàm mật độ r của |Z| là: 2 x2 f (x) = exp(− ) với x ∈ R+ . π 2 Ta đã biết cách xây dựng một mẫu có phân phối mũ. Vì vậy ta chọn hàm g là hàm mật độ của phân phối mũ với tham số 1 và có: r r r 2 x2 − 2x 2e (x − 1)2 2e f (x) = exp(− )= exp(− )≤ . g(x) π 2 π 2 π r 2e f (x) (x − 1)2 Đặt M = . Khi đó = exp(− ). π M.g(x) 2 Theo phương pháp lấy mẫu loại trừ thì ta có thể thực hiện theo các bước sau: Bước 1: Đặt Y = y từ phân bố mũ với tham số 1 và U = u từ phân phối chuẩn U (0, 1). (y − 1)2 ) Bước 2: Nếu u ≤ exp(− ) trở lại X = y . Ngược lại quay về Bước 2 1. ⊗ 9 Ví dụ 1.5 Cho một biến ngẫu nhiên Y với hàm mật độ g(x) được xây dựng trên không gian S. Giả sử A ⊂ S và ta lấy mẫu của biến ngẫu nhiên có điều kiện X = (Y |Y ∈ A) với không gian xác định A. Trong trường hợp lấy mẫu bằng phương pháp loại trừ ta có thể thực hiện bởi việc lấy mẫu X lặp lại cho đến khi mẫu ta cần thuộc A. Chính g(x) xác hơn, X có hàm mật độ f (x) = với x ∈ A. P (Y ∈ A) f (x) 1 f (x) Vậy ≤ = M và = 1[x∈A] . g(x) P (Y ∈ A) M.g(x) Giả sử rằng U là phân phối đều trong một khoảng (0,1). Khi đó  f (Y ) 1 nếu Y ∈ A )= P (U ≤ 0 nếu Y ∈ / A. M.g(Y ) Với phương pháp lấy mẫu loại trừ ta chấp nhận giá trị này nếu Y ∈ A và ngược lại ta loại bỏ. Nếu việc đánh giá hàm f không khó khăn, ngoài việc đánh giá cận trên Mg(x) cho f(x) ta có thể đánh giá cận dưới h(x) và ta có phương pháp lấy mẫu loại bỏ cải biên như sau: 1. Lấy Y = y từ g(y) và U = u từ U (0, 1). h(y) 2. Chấp nhận nếu u ≤ và nhận X = y như là mẫu. Ngược lại M.g(y) chuyển sang bước 3. f (y) và nhận X = y như là mẫu. Ngược lại 3. Chấp nhận nếu u ≤ M.g(y) quay về bước 1. 1 R Sẽ có lợi hơn rất nhiều vì trung bình ta chỉ cần lần nếu M h(x)dx thay bởi hàm h. Hàm h có thể tìm được nhờ khai triển Taylor. 1.2 Ước lượng bằng mô phỏng Trong phần trước, ta tìm cách lấy mẫu của phân bố mục tiêu (target density) dựa trên việc lấy mẫu của một phân bố đề xuất (proposal density). Phần này, dựa trên mẫu của phân bố đề xuất ta tìm ước lượng không chệch cho các đặc trưng của phân bố mục tiêu. 10 1.2.1 Lấy mẫu quan trọng (Importance Sampling) Giả sử ta quan tâm đến tích phân: Z I = Ef (h(x)) = h(x)f (x)dx. S trong đó f là hàm mật độ. Viết tích phân này dưới dạng: Z f (x) h(x)g(x)dx. I= g(x) S trong đó g(x) là hàm mật độ thỏa mãn g(x) > 0 nếu f (x).h(x) 6= 0. Ta sinh mẫu ngẫu nhiên (độc lập có cùng phân bố) (x1 , x2 , ..., xn ) từ g và ước lượng I bởi: n n 1 X f (xi ) 1X Ib = h(xi ) = w(xi )h(xi ). n i=1 g(xi ) n i=1 Thủ tục này gọi là Lấy mẫu quan trọng. Hàm mật độ g được gọi là hàm mật độ đề xuất và Ib gọi là ước lượng quan trọng (importance estimator). f (xi ) là trọng số (importance weight). w(xi ) = g(xi ) Chú ý rằng Ib là một ước lượng không chệch của I. Có hai lý do mà ta quan tâm đến phương pháp này là: 1, Việc lấy mẫu từ f (x)là không thể hoặc rất khó. 2, Biến ngẫu nhiên H(X), trong đó X có mật độ f(x), có phương sai n 1X h(xi ) có sai số lớn. Vì vậy ước lượng không chệch truyền thống, tức n i=1 Monte-Carlo lớn. Phương sai của ước lượng Ib theo phương pháp này là hữu hạn nếu: f 2 (x) f (X) Eg (h (X). 2 ) = Ef (h2 (X). ) < ∞. g (x) g(X) 2 11 f (x) là không bị chặn. Vì g(x) vậy, nếu có thể ta chọn hàm mật độ đề xuất g có đuôi dày hơn f. Tổng f (x) quát, nếu không bị chặn thì cho dù phương sai của ước lượng là hữu g(x) hạn, thì thủ tục này cũng không hiệu quả vì phương sai của các trọng số lớn. Các trọng số biến đổi nhiều sẽ ảnh hưởng đến giá trị ước lượng. Ví dụ 1.6 Các xác suất đuôi của phân bố chuẩn. Cho p = P [Z > 4] trong đó Z là biến ngẫu nhiên có phân phối chuẩn tắc. Dùng cách tiếp cận đơn giản nhất ta có thể xây dựng mẫu ngẫu nhiên độc lập có phân phối chuẩn Z (1) , ..., Z (M ) và đặt: Do đó phương sai sẽ bằng vô cùng nếu tỉ số M 1 X pb = 1 (i) . M 1 [Z >4] Tuy vậy M có thể rất lớn để ước lượng khác 0. Một cách thay thế là xây dựng một mẫu độc lập từ biến ngẫu nhiên có phân phối mũ với tham số 1 và dịch sang phải 4 đơn vị, khi đó: 1 1 f (x) = √ exp(− x2 + (x + 4)). w(x) = g(x) 2 2π n 1X b Thay cho ước lượng quan trọng I = w(xi h(xi )) người ta dùng ước n i=1 lượng: n X Ib = h(xj ).w(xj ) j=1 n X . w(xj ) j=1 Ước lượng này có 2 thuận lợi: 1, Nó là ước lượng chệch có phương sai nhỏ hơn ước lượng quan trọng được xây dựng từ trước. Chú ý rằng ước lượng này vẫn là ước lượng vững khi x1 , ...xn là độc lập cùng phân bố với hàm mật độ g ta có: n 1 X f (xj ) −→ 1 khi n −→ ∞ n j=1 g(xj ) 12 2, Ta có thể có mẫu quan trọng nếu chỉ biết f (x) và w(x) sai khác một hằng số nhân. Nếu lý do là không tìm được hàm mật độ của mẫu quan trọng dẫn đến phương sai nhỏ thì có một vài phương pháp làm giảm phương sai như: Phương pháp 1 (Sequential Importance Resampling). f (Y (i) ) (1) (n) +) Sinh mẫu dạng Y , ...Y với trọng số wi = , i = 1, ...n. g(Y (i) ) wj +) Sinh mẫu mới X (1) , ..., X (n) bởi Y (1) , ..., Y (n) xác suất n . X wi j=1 Phương pháp 2 (Rejection Control). Xem xét loại bỏ các điểm mẫu có trọng số (importance weight) dưới ngưỡng c. Làm theo cách này sẽ sinh ra ước lượng chệch, tuy nhiên nếu thay đổi các trọng số một cách phù hợp ta sẽ tránh được điều này. Giả sử ta có mẫu quan trọng Y (1) , ..., Y (n) với trọng số w1 , ..., wn thì phương pháp này được thực hiện như sau: wj +) Với j = 1, ..., n chấp nhận Y (i) với xác suất Pj = min{1; }. c (j) +) Nếu Y được chấp nhận ta tính lại trọng số: e= w q.wj . pj w(x) }g(x)dx. c Chú ý rằng q giống nhau với mọi j nên ta không cần tính cụ thể nếu ta dùng ước lượng tỷ số. Hơn nữa, phương pháp này sinh ra mẫu quan trọng f (x) min{g(x), } c ∗ . từ hàm mật độ đề xuất g (x) = q trong đó q = 1.3 R min{1, Xích Markov Xét một hệ nào đó được quan sát tại các thời điểm rời rạc 0, 1, 2....Giả sử các quan sát đó là X0 , X1 ..., Xn , .... Khi đó ta có một dãy các đại lượng ngẫu nhiên (ĐLNN) (Xn ) là trạng thái của hệ tại thời điểm n. Giả thiết rằng mỗi Xn với n = 0, 1, 2.... là một ĐLNN rời rạc. Kí hiệu S là tập các 13 giá trị của các (Xn ). Khi đó S là một tập hữu hạn hay đếm được, các phần tử của nó được kí hiệu là i, j, k....Ta gọi S là không gian trạng thái của dãy. 1.3.1 Giới thiệu về xích Markov Định nghĩa 1.1. Xích Markov Dãy các ĐLNN (Xn ) là một xích Markov nếu với mọi n1 < n2 < ... < nk < nk+1 và với mọi i1 , i2 , ...ik , ik+1 ∈ S ta có: P (Xnk+1 = ik+1 |Xn1 = i1 , Xn2 = i2 , ..., Xnk = ik ) = P (Xnk+1 = ik+1 |Xnk = ik ). Ta coi thời điểm nk+1 là tương lai, nk là hiện tại và n1 , n2 , ..., nk−1 là quá khứ. Như vậy xác suất có điều kiện của một sự kiện B nào đó trong tương lai nếu biết hiện tại và quá khứ của hệ cũng giống như xác suất có điều kiện của B nếu chỉ biết trạng thái hiện tại của hệ. Đó chính là tính Markov của hệ. Hay nói cách khác, tính Markov cũng được phát biểu dưới dạng: Nếu biết trạng thái hiện tại của hệ thì quá khứ và tương lai độc lập với nhau. Giả sử P (Xm+n = j|Xm = i) là xác suất để xích tại thời điểm m ở trạng thái i, tại thời điểm m + n thì chuyển sang trạng thái j sau n bước. Xác suất này nói chung phụ thuộc vào i, j, m, n. Nếu đại lượng này không phụ thuộc vào m ta nói xích thuần nhất. Từ giờ trở đi, ta chỉ xét các xích thuần nhất. Kí hiệu: Pij = P (Xn+1 = j|Xn = i). Pij (n) = P (Xm+n = j|Xm = i). Ta gọi (Pij , i, j ∈ S) là xác suất chuyển sau một bước (xác suất chuyển)còn (Pij (n), i, j ∈ S) là xác suất chuyển sau n bước. Chú ý rằng: X Pij = 1. j∈S X Pij (n) = 1. j∈S 14 Phân bố của X0 được gọi là phân bố ban đầu. Ta kí hiệu ui = P (X0 = i). Định lý 1.2. Phân bố đồng thời của (X0 , X1 , ..., Xn ) được hoàn toàn xác định từ phân bố ban đầu và xác suất chuyển: P (X0 = i0 , X1 = i1 , ..., Xn = in ) = ui0 Pi0 i1 ... Pin−1 in . Chứng minh: Thật vậy, theo công thức nhân xác suất ta có: P (X0 = i0 , X1 = i1 , ..., Xn = in ) = P (X0 = i0 )P (X1 = i1 |X0 = i0 )× P (Xk = ik |X0 = i0 , ..., Xk−1 = ik−1 )× P (Xn = in |X0 = i0 , ..., Xn−1 = in−1 ). Sử dụng tính Markov ta có: P (Xk = ik |X0 = i0 , ..., Xk−1 = ik−1 ) = P (Xk = ik |Xk−1 = ik−1 ) = Pik−1 ik . Do đó: P (X0 = i0 , X1 = i1 , ..., Xn = in ) = ui0 Pi0 i1 ... Pin−1 in . Định lý 1.3. (Phương trình C-K(Chapman-Kolmogorov)) Pi+j (n + m) = X Pik (n)Pkj (m). k∈S Trong trường hợp S có d phần tử, ta kí hiệu P = (Pij ), P (n) = (Pij (n)) là các ma trận vuông cấp d × d. P được gọi là ma trận xác suất chuyển, P(n) được gọi là ma trận xác suất chuyển sau n bước. Khi đó phương trình Chapman-Kolmogorov tương đương với: P (n + m) = P (n) · P (m). vì P = P (1) nên bằng quy nạp ta dễ thấy: 15 P (n) = P n . Gọi ui (n) = P (Xn = i) Kí hiệu véc tơ U (n) = (u1 (n), u2 (n), ..., ud (n)) là véc tơ hàng d-chiều mô tả phân bố của Xn và U = U (0) = (u1 , u2 , ..., ud ) là véc tơ hàng d-chiều mô tả phân bố ban đầu (phân bố của X0 ). Định lý 1.4. Ta có: P (n + m) = U (m) · P n Nói riêng U (n) = U · P n Ví dụ 1.7 Cho (ξn ), n = 0, 12, ... là dãy các ĐLNN độc lập có cùng phân bố. Giả sử P (ξn = i) = ai , i ∈ Z. Đặt: Xn = n X ξi . i=0 Khi đó Xn là một xích Markov với không gian trạng thái Z. Thật vậy: P (Xn+1 = = = = in+1 |X0 = i0 , X1 = i1 , ..., Xn = in ) P (Xn + ξn+1 = in+1 |ξ0 = i0 , ξ1 = i1 − i0 , ..., ξn = in − in−1 ) P (ξn+1 = in+1 − in ) P (Xn+1 = in+1 |Xn = in ). Do đó (Xn ) là một xích Markov với không gian trạng thái Z, xác suất chuyển là Pij = (aj−i ). ⊗ Ví dụ 1.8 Giả sử trung bình trượt của một dãy Bernoulli Yn = đó (Xn ) là dãy các BNN Bernoulli độc lập với: 16 Xn + Xn−1 trong 2 1 P (Xn = 0) = P (Xn = 1) = (∀n = 0, 1, 2, ...). 2 Ta sẽ chỉ ra rằng (Yn ) là một quá trình không có tính Markov Thật vậy, hàm xác suất của Yn được cho bởi: 1 P (Yn = 0) = P (Xn = 0, Xn−1 = 0) = . 4 1 1 P (Yn = ) = P (Xn = 0, Xn−1 = 1) + P (Xn = 1, Xn−1 = 0) = . 2 2 1 P (Yn = 1) = P (Xn = 1, Xn−1 = 1) = . 4 Xác suất có điều kiện của 2 giá trị liên tiếp của (Yn ) là: P (Yn = 1|Yn−1 1 P (Yn = 1, Yn−1 = ) 1 2 = ) = 1 2 P (Yn−1 = ) 2 P (Xn = 1, Xn−1 = 1, Xn−2 = 0) = 1 P (Yn−1 = ) 2 1 3 ( ) 1 = 2 = . 1 4 2 Mặt khác 1 P (Yn = 1|Yn−1 = , Yn−2 2 1 P (Yn = 1, Yn−1 = , Yn−2 = 1) 2 = 1) = = 0. 1 P (Yn−1 = , Yn−2 = 1) 2 (do không có dãy nào của Xn , Xn−1 Xn−2 , Xn−3 dẫn tới kết quả Yn = 1 1, Yn−1 = , Yn−2 = 1), nên ta có: 2 1 1 P (Yn = 1|Yn−1 = , Yn−2 = 1) 6= P (Yn = 1|Yn−1 = ). 2 2 Vậy quá trình này không có tính Markov. 17 1.3.2 Phân bố dừng Định nghĩa 1.2. Phân bố dừng Phân bố ban đầu U = (ui ), i ∈ S là phân bố dừng nếu ta có U (n) = U với ∀n tức là ui (n) = ui với ∀n ,∀i ∈ S . Khi đó dãy (Xn ) có cùng phân bố. Ta cũng có kết quả tương đương sau: U = (ui ) là phân X bố dừng nếu và chỉ nếu: 1. ui > 0 và ui = 1, X i∈S 2. uj = ui Pij với ∀j ∈ S. i∈S Ví dụ 1.9 Cho (Xn ) là xích Markov với không gian trận xác suất chuyển P:  3 1 5 5   1 2 P = 6 3   3 3 8 8 trạng thái S = 1, 2, 3 và ma  1 5   1  6   1 4 Hãy tìm phân bố dừng (nếu có) của xích đó. Lời giải: Đặt U = (x, y, z). Khi đó U là phân bố dừng khi và chỉ khi x, y, z là nghiệm không âm của hệ sau :  1 3 3   x+ y+ z =x   5 6 8      1 2 3 x+ y+ z =y  5 3 8        1x + 1y + 1z = z 5 6 4 18 Giải ra ta được x = 15 18 8 ,y = ,z = . 41 41 41 Vậy phân bố dừng của xích là U = (x = 15 18 8 , y = , z = ). 41 41 41 Ví dụ 1.10 Cho (Xn ) là một xích Markov có hai trạng thái S = {0, 1} với các ma trận xác suất chuyển là:   1−a a P = b 1−b trong đó a, b > 0, 0 < a + b < 1. Ta đi tìm phân bố dừng Đặt U = (x, y). Khi đó U là phân bố dừng khi và chỉ khi x, y là nghiệm không âm của hệ sau:  (1 − a)x + by = x      ax + (1 − b)y = y      x+y =1 Phương trình thứ nhất và thứ hai của hệ tương đương với ax = by hay by x= . a b a Thế vào phương trình thứ ba của hệ ta có x = ,y = . a+b a+b Vậy phân bố dừng của xích là: U =( b a , ). a+b a+b Khi nghiên cứu sâu hơn ta có thể thấy phân bố dừng không phải bao giờ cũng tồn tại. Vậy ta đặt ra câu hỏi: Với điều kiện nào thì tồn tại phân bố dừng và phân bố dừng nếu tồn tại thì có duy nhất không ? Định lý 1.5. Giả sử (Xn ) là một xích Markov với không gian trạng thái S = {1, 2, ...} với xác suất chuyển sau n bước là P (n) = (Pij (n)). Giả sử 19 rằng mọi i, j ∈ S tồn tại giới hạn lim Pij (n) = πj n→∞ và X giới hạn này không X phụ thuộc i. Khi đó: πj ≤ 1 và πj = πi Pij . 1. j∈S j∈S 2. Hoặc πj = 0 với mọi j ∈ S ,hoặc X πj = 1. j∈S 3. Nếu X πj = 1 thì U = (π1 , π2 , ...) là phân bố dừng và là phân bố dừng j∈S duy nhất. Nếu πj = 0 với mọi j ∈ S thì phân bố dừng không tồn tại. 1.3.3 Phân bố giới hạn Định nghĩa 1.3. Phân bố giới hạn Giả sử (Xn ) là xích Markov với không gian trạng thái S = {1, 2, ...} với ma trận xác suất chuyển P = (Pij ) và ma trận xác suất chuyển sau n bước là P (n) = Pij (n). Ta nói rằng xích có phân bố giới hạn nếu với ∀i, j ∈ S tồn tại giới hạn: lim Pij (n) = πj . n→∞ Giới hạn này không phụ thuộc i ∈ S và X πj = 1. Nói cách khác, véc tơ j∈S giới hạn π = (π1 , π2 , ...) lập thành một phân bố xác suất trên S. Ý nghĩa của phân bố giới hạn như sau: Gọi ui (n) = P (Xn = i). Kí hiệu véc tơ U (n) = (u1 (n), u2 (n), ...) là véc tơ hàng d-chiều mô tả phân bố của Xn . Ta có: X P (Xn = j) = P (X0 = i)Pij (n). i∈S Do đó: lim P (Xn = j) = n→∞ X i∈S = X i∈S 20 P (X0 = i) lim Pij (n) n→∞ P (X0 = i)Πj = Πj .