Tài liệu Xích markov và thuật toán metropolis

BỘ GIÁO DỤC VÀ ĐÀO TẠO TRƯỜNG ĐẠI HỌC SƯ PHẠM HÀ NỘI 2 KHOA TOÁN Nguyễn Thị Thanh Hà XÍCH MARKOV VÀ THUẬT TOÁN METROPOLIS KHÓA LUẬN TỐT NGHIỆP ĐẠI HỌC Hà Nội – Năm 2017 BỘ GIÁO DỤC VÀ ĐÀO TẠO TRƯỜNG ĐẠI HỌC SƯ PHẠM HÀ NỘI 2 KHOA TOÁN Nguyễn Thị Thanh Hà XÍCH MARKOV VÀ THUẬT TOÁN METROPOLIS Chuyên ngành: Toán Ứng dụng KHÓA LUẬN TỐT NGHIỆP ĐẠI HỌC NGƯỜI HƯỚNG DẪN KHOA HỌC: TS. TRẦN VĨNH ĐỨC Hà Nội – Năm 2017 Mục lục Lời cảm ơn 1 Lời cam đoan 2 Lời mở đầu 3 1 Kiến thức chuẩn bị 5 1.1 1.2 Biến ngẫu nhiên và phân phối xác suất . . . . . . . . . 5 1.1.1 Biến ngẫu nhiên . . . . . . . . . . . . . . . . . . 5 1.1.2 Hàm phân bố xác suất . . . . . . . . . . . . . . 7 1.1.3 Phân loại biến ngẫu nhiên . . . . . . . . . . . . 9 Mô phỏng biến ngẫu nhiên . . . . . . . . . . . . . . . . 10 1.2.1 Các phương pháp chung . . . . . . . . . . . . . 10 1.2.2 Phương pháp Monte Carlo . . . . . . . . . . . . 14 2 Xích Markov 17 2.1 Định nghĩa . . . . . . . . . . . . . . . . . . . . . . . . 17 2.2 Ma trận chuyển trạng thái . . . . . . . . . . . . . . . . 20 2.3 Phân bố dừng . . . . . . . . . . . . . . . . . . . . . . . 24 2.4 Phân bố giới hạn . . . . . . . . . . . . . . . . . . . . . 31 2.5 Một số tính chất quan trọng và định nghĩa liên quan . 32 i Khóa luận tốt nghiệp Đại học Nguyễn Thị Thanh Hà 2.5.1 Tính Ergodic . . . . . . . . . . . . . . . . . . . 32 2.5.2 Tính khả nghịch . . . . . . . . . . . . . . . . . 35 3 Thuật toán Metropolis 37 3.1 Giới thiệu xích Markov Monte Carlo . . . . . . . . . . 37 3.2 Thuật toán Metropolis - Hasting . . . . . . . . . . . . 38 3.2.1 Giới thiệu . . . . . . . . . . . . . . . . . . . . . 38 3.2.2 Thuật toán . . . . . . . . . . . . . . . . . . . . 39 3.2.3 Ví dụ 40 . . . . . . . . . . . . . . . . . . . . . . . Kết luận 46 Tài liệu tham khảo 47 ii Khóa luận tốt nghiệp Đại học Nguyễn Thị Thanh Hà Lời cảm ơn Em hoàn thành khóa luận tốt nghiệp dưới sự hướng dẫn và chỉ bảo tận tình của Ts. Trần Vĩnh Đức. Thầy cũng dành thời gian quý báu của mình để hướng dẫn cũng như giải đáp thắc mắc của em trong quá trình làm bản khóa luận này. Em muốn tỏ lòng biết ơn và gửi lời cảm ơn sâu sắc nhất tới thầy! Em cũng xin chân thành cảm ơn các thầy cô giáo Khoa Toán trường Đại học Sư phạm Hà Nội 2, đặc biệt là tổ Ứng dụng, đã tạo điều kiện thuận lợi cho em trong quá trình học Đại học và thực hiện bản khóa luận này. Hà Nội, ngày 24/04/2017 Tác giả khóa luận Nguyễn Thị Thanh Hà 1 Khóa luận tốt nghiệp Đại học Nguyễn Thị Thanh Hà Lời cam đoan Khóa luận tốt nghiệp "Xích Markov và thuật toán Metropolis" được hoàn thành do sự cố gắng, nỗ lực tìm hiểu nghiên cứu của bản thân cùng với sự hướng dẫn tận tình của TS.Trần Vĩnh Đức. Tôi xin cam đoan khóa luận tốt nghiệp này không trùng lặp với kết quả của các tác giả khác. Hà Nội, ngày 24/04/2017 Tác giả khóa luận Nguyễn Thị Thanh Hà 2 Khóa luận tốt nghiệp Đại học Nguyễn Thị Thanh Hà Lời mở đầu Thuật toán Metropolis được John von Neumann, Stan Ulam và Nick Metropolis xây dựng vào năm 1946 còn được biết đến với tên là phương pháp Monte Carlo. Thuật toán thuật toán này đưa ra một cách thức hiệu quả để đi “đến gần” lời giải cho các bài toán quá phức tạp, khó có thể giải một cách chính xác. Đượ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. Ý 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. Khóa luận gồm ba chương. Chương 1 "Kiến thức chuẩn bị" trình bày một số các hiểu biết biến ngẫu nhiên, hàm phân phối xác suất đồng thời trình bày về mô phỏng biến ngẫu nhiên trong đó có các phương pháp chung và phương pháp Monte Carlo. Chương 2 "Một vài hiểu biết về xích Markov" trình bày một số định nghĩa cơ bản nhất về xích Markov cũng như các tính chất của xích. Chương 3 "Thuật toán Metropolis" thảo luận một số nội dung như xích Markov Monte Carlo (MCMC) là phương pháp Monte Carlo kết 3 Khóa luận tốt nghiệp Đại học Nguyễn Thị Thanh Hà hợp với xích Marko và thuật toán Metropolis là một trường hợp đặc biệt của phương pháp MCMC ở trên. Bên cạnh đó trong chương này còn trình bày một số ví dụ để làm rõ ứng dụng của thuật toán Metropolis. 4 Chương 1 Kiến thức chuẩn bị Thuật toán Metropolis là sự kết hợp của phương pháp mô phỏng biến ngẫu nhiên với xích Markov, chính vì vậy trong chương này trình bày các kiến thức cơ bản về biến ngẫu nhiên, hàm phân phối xác suất và phương pháp mô phỏng biến ngẫu nhiên. Đây là các kiến thức nền tảng để ta đi xây dựng thuật toán Metropolis. 1.1 1.1.1 Biến ngẫu nhiên và phân phối xác suất Biến ngẫu nhiên Định nghĩa 1.1. Một không gian xác suất (Ω, F, P ) là một không gian được trang bị một độ đo với độ đo toàn thể bằng 1 (nghĩa là P (Ω)=1). Trong đó, thành phần đầu, Ω được xem không gian mẫu, là một tập không rỗng, với các phần tử thường được biết như là các "kết quả" hay "trạng thái tự nhiên" (ví dụ trạng thái sấp hay ngửa của đồng tiền,...). Một trạng thái tự nhiên luôn tồn tại với một xác suất nào đó. Một phần tử của Ω thường được ký hiệu bởi ω. 5 Khóa luận tốt nghiệp Đại học Nguyễn Thị Thanh Hà Thành phần thứ hai, F, là một tập hợp mà các phần tử của nó được gọi là các sự kiện (event). Các sự kiện là các tập con của Ω. Tập F phải thỏa mãn một vài điều kiện, đặc biệt nó phải là một σ-đại số. Cùng với nhau, Ω và F tạo thành một không gian đo được. Một sự kiện là một tập hợp các kết quả hay trạng thái tự nhiên mà ta có thể xác định xác suất của nó. Thành phần thứ ba, P , được gọi là "độ đo xác suất", hay "xác suất". Nó là một hàm số từ F vào tập số thực, cho tương ứng mỗi sự kiện với một xác suất có giá trị nằm giữa 0 và 1. Nó cần thỏa mãn các điều kiện, đó là nó phải là một độ đo và P (Ω) = 1. Các độ đo xác suất thường được viết đậm có gạch, ví dụ P hay Q Định nghĩa 1.2. Một biến ngẫu nhiên là một ánh xạ từ Ω đến một tập hợp khác, thông thường là tập các số thực. Đặc biệt, nó phải là một hàm đo được. Điều này có nghĩa là, ví dụ, nếu X là một biến ngẫu nhiên thực, thì tập hợp các kết quả sao cho X là dương, {ω ∈ Ω : X(ω) > 0}, phải là một sự kiện. Người ta thường tóm gọn cách viết {ω ∈ Ω : X(ω) > 0} thành {X > 0} và viết P (X > 0) thay vì viếtP ({X > 0}). Nếu biến ngẫu nhiên X chỉ nhận hữu hạn giá trị thì nó được gọi là biến ngẫu nhiên đơn giản. Biến ngẫu nhiên còn được gọi là đại lượng ngẫu nhiên. Đối với biến ngẫu nhiên người ta chỉ quan tâm xem nó có nhật một giá trị nào đó hoặc nhận trong một khoảng nào đó với xác suất là bao nhiêu? Ví dụ về các đại lượng ngẫu nhiên: 6 Khóa luận tốt nghiệp Đại học Nguyễn Thị Thanh Hà • Số chấm xuất hiện khi ta gieo một con xúc xắc, • Số người mua hàng ở một của hàng ở một khoảng thời gian nào đó, • Số lần bắn trúng tâm của một xạ thủ trong một cuộc thi nào đó. 1.1.2 Hàm phân bố xác suất Định nghĩa 1.3. Giả sử (Ω, F, P ) là một không gian xác suất, X : Ω → R là biến ngẫu nhiên. Khi đó hàm số FX (x) = P (X < x) = P (ω : X(ω) < x) được gọi là hàm phân bố xác suất của X. Nhận xét: FX (x) = P [X −1 (−∞, x)] = PX [(−∞, x)]. Hàm phân bố xác suất có các tính chất sau: • 0 ≤ FX (x) ≤ 1 với mọi x ∈ R, • FX (x) là hàm không giảm liên tục phải, • FX (−∞) = lim FX (x) = 0; FX (+∞) = lim FX (x) = 1, x→−∞ x→+∞ • P {a < X ≤ b} = FX (a) − FX (b), • P {X > a} = 1 − FX (a) ; P {X < a} = FX (a− ) với FX (a− ) = lim FX (x). x 1/2. Tính P {X ≤ 1/4}, P {0 < X ≤ 1/4}, P {X = 0} và P {0 ≤ X ≤ 1/4}. Lời giải. • P {X ≤ 1/4} = F (1/4) = 1/4 + 1/2 = 3/4, • P {0 < X ≤ 1/4} = F (1/4) − F (0) = 3/4 − 1/2 = 1/4, • P {X = 0} =P {X ≤ 0} - P {X < 0} = F (0)−F (0− ) = 1/2−0 = 1/2, • P {0 ≤ X ≤ 1/4} = P {X = 0} + P {0 < X ≤ 1/4} = 1/2+3/4− 1/2 = 3/4. Ví dụ 1.1.2. Chọn ngẫu nhiên 3 viên bi từ một hộp gồm 6 bi đen và 4 bi trắng. Gọi X là số bị trắng trong số bi vừa chọn thì X là biến ngẫu nhiên rời rạc. Tìm hàm phân bố xác suất của X. Lời giải. Vì ta lấy ra 3 viên bi từ hộp bên trong có 6 bi đen và 4 bi trắng ở trên nên số bi trắng có thể có là 0 viên, 1 viên, 2 viên hoặc 3 viên nên X ∈ {0, 1, 2, 3}. Ta có: C63 P {X = 0} = 3 = 5/30 = 1/6, C10 C62 C41 P {X = 1} = = 15/30 = 1/2, 3 C10 8 Khóa luận tốt nghiệp Đại học Nguyễn Thị Thanh Hà C61 C42 = 9/30 = 3/10, 3 C10 C 0C 3 P {X = 1} = 6 3 4 = 1/30. C10 Vậy ta có hàm phân bố xác suất là : P {X = 2} = F (x) = 1.1.3    0         5/30    khi x < 0, khi 0 ≤ x < 1, khi 1 ≤ x < 2, 20/30       29/30       1 nếu 2 ≤ x < 3, nếu x ≥ 3. Phân loại biến ngẫu nhiên Biến ngẫu nhiên được chia làm 2 loại đó là biến ngẫu nhiên rời rạc và biến ngẫu nhiên liên tục. Định nghĩa 1.4. Biến ngẫu nhiên được gọi là biến ngẫu nhiên rời rạc nếu nó chỉ nhận một số hữu hạn hoặc đếm được các giá trị. Đối với biến ngẫu nhiên rời rạc, tập hợp tất cả các giá trị có thể có của nó có thể được liệt kê bằng một dãy hữu hạn hoặc vô hạn x1 , x2 , ..., xn , ... Tập hợp các giá trị có thể có của biến ngẫu nhiên X được ký hiệu là X(Ω). Định nghĩa 1.5. Biến ngẫu nhiên được gọi là biến ngẫu nhiên liên tục nếu hàm phân phối F (x) của nó là hàm liên tục và tồn tại hàm số p(x) sao cho: 1. p(x) ≥ 0, với mọi x ∈ R, 9 Khóa luận tốt nghiệp Đại học 2. F (x) = +∞ R Nguyễn Thị Thanh Hà p(t)dt, với mọi x ∈ R. −∞ Hàm p(x) nêu trên được gọi là hàm mật độ xác suất của X. Từ định nghĩa suy ra: 1. Với mọi a, b thỏa mãn −∞ ≤ a < b ≤ +∞ ta có P (a < X < b) = Rb p(x)dx, a 2. +∞ R p(x)dx = 1, −∞ 3. p(x) = F (x) tại mọi điểm x mà p(x) liên tục. 1.2 Mô phỏng biến ngẫu nhiên Mô phỏng được ứng dụng rộng rãi trong kinh tế, kỹ thuật, vật lý và nhiều lĩnh vực khác. Theo từ điển Oxford năm 1976: "Mô phỏng có nghĩa là giả cách, làm ra như vẻ, hành động như, bắt chước như, làm giả các điều kiện của tình huống nào đó để thông qua một mô hình với mục đích huấn luyện hoặc tiện lợi". Về mặt ý nghĩa kỹ thuật, mô phỏng (hay đúng hơn là phương pháp mô phỏng) hàm chứa việc áp dụng một mô hình nào đó để tạo ra kết quả chứ không phải là thử nghiệm một hệ thống thực tế nào đó để tạo ra kết quả. 1.2.1 Các phương pháp chung a, Phương pháp ngược 10 Khóa luận tốt nghiệp Đại học Nguyễn Thị Thanh Hà Mệnh đề 1.1. Giả sử biến ngẫu nhiên có hàm phân phối xác suất F (x) đều trên (0,1). Ta thấy rằng nếu giả sử cho biến ngẫu nhiên X liên tục và F (x) liên tục, đơn điệu tăng khi 0 ≤ F (x) ≤ 1. Khi đó tồn tại hàm ngược F −1 . Như vậy để tạo X ta có thể dùng thuật toán như sau: 1, Tạo U ∼ U (0, 1), 2, Lấy X=F −1 (U ). Dễ dàng thấy X có phân phối F , thật vậy chúng ta thấy rằng bất kỳ số thực x nào đó thì P (X ≤ x) = F (x) và do đó F có hàm ngược nên ta có: P (X ≤ x) = P [F −1 ≤ x] = P [U ≤ F (x)] = F (x) ở đây U ∼ U (0, 1). Mệnh đề 1.2. Giả sử hàm F được được xác định bởi: F (U ) = min {x : F (x) ≥ u}, Khi đó nếu U ∼ U (0, 1) thì X= F (U ) tuân theo phân phối định nghĩa bởi F . Sử dụng mệnh đề này chúng ta có thể tạo ra biến số rời rạc rồi chỉ lấy một giá trị xác định x1 , x2 , ... với x1 ≤ x2 ≤ .... Trong trường hợp này: F (x) = P (X ≤ x) = P x1 ≤x p(xi ), ở đây p(xi ) = P (X = xi ). Thuật toán được mô tả như sau: (1). Tạo U ∼ U (0, 1), (2). Xác định số nguyên nhỏ nhất I sao cho F (xI ) ≥ U , (3). Lấy X = x1 . 11 Khóa luận tốt nghiệp Đại học Nguyễn Thị Thanh Hà Một số ví dụ: Ví dụ 1.2.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 khoảng (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 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ố λ. Ví dụ 1.2.2. Mô phỏng biến ngẫu nhiên có phân phối Bernoulli (p) và biến ngẫu nhiê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: P =   1 nếu U < p,  0 nếu 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 , X2 , ..., Xn là mội mẫu độc lập cùng phân phối Bernoulli. n P Khi đó Y = Xi có phân phối nhị thức B(n, p). 1 12 Khóa luận tốt nghiệp Đại học Nguyễn Thị Thanh Hà b, Phương phấp chấp nhận - loại bỏ Phương pháp này đòi hỏi tồn tại một hàm g mà không nhỏ hơn hàm mật độ f , tức là g(x) ≥ f (x) với mọi x. Nói chung hàm g không g(x) phải là hàm mật độ nhưng hàm h được định nghĩa bởi h(x) = c là hàm mật độ. Giả sử biến ngẫu nhiên có hàm mật độ h có thể được tạo ra một cách hiệu quả. Khi đó thuật toán được mô tả như sau: (1). Tạo y từ h, (2). Tạo U ∼ U (0, 1), độc lập với y, f (x) (3). Nếu U ≥ lấy X = g. Ngược lại, chuyển về bước 1. g(x) Phương pháp này hiệu quả với c nhỏ. Ví dụ 1.2.3. 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 x2 exp(− ) với mọi x ∈ R+ . f (x) = π 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 r 2 x2 − 2x 2e (x − 1)2 2e exp(− )= exp(− )≤ , π 2 π 2 π r 2e Từ đó, đặt M = dẫn đến π f (x) (x − 1)2 = exp(− ). M g(x) 2 f (x) = g(x) r Thuật toán lấy mẫu bằng phương pháp chấp nhận - 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), 13 Khóa luận tốt nghiệp Đại học Nguyễn Thị Thanh Hà • Bước 2. Nếu u ≤ exp(− (y − 1)2 ) thì X = y. Ngược lại quay trở 2 về bước 1. c, Phương pháp phức hợp Phương pháp này thường được sử dụng khi hàm phân phối F có thể biểu hiện như một tổ hợp lồi cái hàm phân phối khác mà những biến của chúng có thể tạo ra. Giả sử F (x) = ∞ P aj Fj (x), j=1 ở đây ai ≥ 0, ∞ P aj = 1 và mỗi Fj là một hàm phân phối. Trong thực j=1 tế chỉ có một số hữu hạn a dương. Tương tự hàm mật độ f của X có thể biểu diễn như sau: f (x) = ∞ P aj fj (x), j=1 ở đây fj là hàm mật độ. Thuật toán được tổng quát như sau: (1). Tạo ra một số nguyên dương J sao cho P (J = j) = a với j = 1, 2, 3, ... , (2). Giả sử rằng J = j, tạo X với phân phối Fj và quay lại. 1.2.2 Phương pháp Monte Carlo Mô phỏng Monte Carlo là một phương pháp sử dụng các yếu tố ngẫu nhiên hoặc các biến ngẫu nhiên có phân phối đều trên khoảng (0,1) để giải quyết các bài toán mà sự trôi của thời gian không đóng vai trò quyết định. 14 Khóa luận tốt nghiệp Đại học Nguyễn Thị Thanh Hà Phương pháp Monte Carlo là một lớp các thuật toán để giải quyết nhiều bài toán trên máy tính theo kiểu không tất định, thường sử dụng các yếu tố ngẫu nhiên, ngược lại với các thuật toán tất định. Trong toán học, thuật toán Monte Carlo là phương pháp tính bằng số hiệu quả cho nhiều bài toán liên quan đến nhiều biến số mà không dễ dàng giải được bằng các phương pháp khác, chẳng hạn bằng tích phân. Thuật toán còn được sử dụng cho nhiều lớp bài toán tối ưu hóa như trong ngành tài chính. Nhiều khi, phương pháp Monte Carlo được thực hiện hiệu quả hơn với số giả ngẫu nhiên, thay cho số ngẫu nhiên thực thụ vốn rất khó tạo ra bởi may tính. Các số giả ngẫu nhiên có tính tất định, tạo ra từ chuỗi giả ngẫu nhiên có quy luật, có thể sử dụng chạy thử hoặc chạy mô phỏng theo cùng điều kiện như trước. Các số giả ngẫu nhiên chỉ cần "đủ mức ngẫu nhiên", nghĩa là chúng theo phân bố đều hoặc theo phân bố định trước. Phương pháp Monte Carlo thường được thực hiện lặp lại một số lượng lớn các bước đơn giản, song song với nhau, một phương pháp phù hợp cho máy tính. Kết quả càng chính xác khi các bước lặp càng tăng lên. Ví dụ 1.2.4. Giả sử muốn đánh giá tích phân I = Rb a g(x)d(x), ở đây g(x) là hàm thực khả tích. Sử dụng mô phỏng Monte Carlo ta có thể đánh giá gần đúng bài toán này. Lời giải. Giả sử X là một biến ngẫu nhiên có phân phối đều trên [a,b] và xét Y = (b − a)g(X). Khi đó kỳ vọng của biến Y có dạng sau: E(Y ) = E(b − a)g(X) 15 Khóa luận tốt nghiệp Đại học Nguyễn Thị Thanh Hà = (b − a)Eg(X) Rb = (b − a) a g(x)fX dx Rb = a g(x)dx = I, ở đây fX = (b − a)−1 là hàm mật độ xác suất của biến ngẫu nhiên có phân phối đều trên [0,1]. Như vậy bài toán đánh giá tích phân được đưa về đánh giá kỳ vọng E(Y ). Chúng ta sẽ đánh giá I = E(Y ) bằng trung bình mẫu: Y (n) = Y1 + Y2 + ... + Yn , n trong đó X1 , X2 , ...,Xn là các biến ngẫu nhiên độc lập, cùng phân phối đều trên (0,1). Chúng ta chứng minh rằng EY (n) = 1 sao cho Y (n) là một ước lượng không chệch của I và V ar[Y (n)] = n−1 V ar(Y ). Giả sử V ar(Y ) là hữu hạn, suy ra rằng Y (n) sẽ dần tới I khi n đủ lớn với xác suất đơn vị. 16