Trong bài viết này, chúng ta sẽ cùng tìm hiểu một cách trực quan và dễ tiếp cận về phương pháp chuỗi Markov Monte Carlo, thường được gọi tắt là MCMC. Đây là công cụ được sử dụng rất phổ biến để khớp các mô hình thống kê Bayes. MCMC có nhiều biến thể khác nhau, nhưng nội dung hôm nay sẽ tập trung sâu vào thuật toán Metropolis-Hastings. Để giữ cho bài viết ngắn gọn và trực quan, một số chi tiết kỹ thuật chuyên sâu sẽ được lược bớt. Bạn nên tham khảo tài liệu hướng dẫn sử dụng phân tích Bayes của Stata để có thêm thông tin chi tiết trước khi áp dụng vào thực tế.
Hãy tiếp tục với ví dụ tung đồng xu từ phần trước về các khái niệm cơ bản của thống kê Bayes. Mục tiêu của chúng ta là tìm ra phân phối hậu nghiệm của tham số theta, đại diện cho xác suất đồng xu xuất hiện mặt ngửa. Phân phối tiền nghiệm được lựa chọn là phân phối beta phẳng, không mang thông tin với các tham số là 1 và 1. Chúng ta sẽ sử dụng hàm hợp lý nhị thức để định lượng dữ liệu thực nghiệm thu được từ quá trình tung đồng xu, cụ thể là có 4 lần ngửa trong tổng số 10 lần tung.
Chúng ta có thể áp dụng phương pháp MCMC kết hợp thuật toán Metropolis-Hastings để tạo ra một mẫu ngẫu nhiên từ phân phối hậu nghiệm của theta. Từ mẫu thu được này, việc ước lượng các đặc trưng như giá trị trung bình của phân phối hậu nghiệm sẽ trở nên dễ dàng.
Kỹ thuật này được cấu thành từ ba phần cơ bản bao gồm phương pháp Monte Carlo, chuỗi Markov và thuật toán Metropolis-Hastings.
Phương pháp Monte Carlo
Thuật ngữ Monte Carlo dùng để chỉ các phương pháp dựa trên việc tạo ra các số giả ngẫu nhiên. Hình 1 dưới đây minh họa một số đặc điểm cơ bản của một thực nghiệm Monte Carlo.
Trong ví dụ này, chúng ta sẽ tạo ra một chuỗi các số ngẫu nhiên từ phân phối chuẩn có trung bình là 0.5 và một phương sai sigma tùy ý. Phân phối chuẩn này được gọi là phân phối đề xuất.
Đồ thị ở góc trên bên phải là đồ thị đường vết, hiển thị các giá trị của theta theo thứ tự được rút ra. Đồ thị này cũng hiển thị phân phối đề xuất được xoay theo chiều kim đồng hồ 90 độ và dịch chuyển sang phải sau mỗi lần rút một giá trị theta mới. Điểm màu xanh lá cây trên đồ thị đường vết biểu thị giá trị hiện tại của theta.
Đồ thị mật độ ở góc trên bên trái là biểu đồ tần suất của mẫu được xoay ngược chiều kim đồng hồ 90 độ. Các trục tọa độ đã được xoay để trục theta của cả hai đồ thị song song với nhau. Điểm màu xanh lá cây trên đồ thị mật độ cũng hiển thị giá trị hiện tại của theta.
Khi tiến hành rút 10,000 giá trị ngẫu nhiên của theta để quan sát cách quá trình này vận hành, chúng ta sẽ nhận ra một số đặc trưng quan trọng. Đầu tiên, phân phối đề xuất không thay đổi từ vòng lặp này sang vòng lặp khác. Thứ hai, đồ thị mật độ ngày càng trông giống với phân phối đề xuất khi kích thước mẫu tăng lên. Thứ ba, đồ thị đường vết có tính dừng, nghĩa là mô hình biến thiên trông giống nhau trên tất cả các vòng lặp.
Mô phỏng Monte Carlo đã tạo ra một mẫu có dạng rất giống với phân phối đề xuất. Trong nhiều trường hợp, đó là tất cả những gì chúng ta cần. Tuy nhiên, để tạo ra một mẫu từ phân phối hậu nghiệm phức tạp, chúng ta sẽ cần thêm các công cụ bổ trợ khác.
Chuỗi Markov Monte Carlo
Một chuỗi Markov là một chuỗi các số mà trong đó mỗi số phụ thuộc trực tiếp vào số đứng ngay trước nó trong chuỗi. Ví dụ, chúng ta có thể rút các giá trị của theta từ một phân phối đề xuất chuẩn có trung bình bằng chính giá trị trước đó của theta.
Khi tiến hành chạy thử nghiệm này, phân phối đề xuất sẽ thay đổi sau mỗi vòng lặp. Điều này tạo ra một đồ thị đường vết có dạng bước đi ngẫu nhiên, nghĩa là độ biến thiên không giống nhau giữa các vòng lặp. Kết quả là đồ thị mật độ thu được không giống với phân phối đề xuất hay bất kỳ phân phối hữu ích nào khác, và chắc chắn nó không giống với phân phối hậu nghiệm mong muốn.
Nguyên nhân khiến phương pháp MCMC đơn thuần thất bại trong việc tạo ra mẫu từ phân phối hậu nghiệm là do chúng ta chưa đưa bất kỳ thông tin nào của phân phối hậu nghiệm vào quá trình tạo mẫu. Chúng ta có thể cải thiện mẫu bằng cách giữ lại các giá trị đề xuất của theta có khả năng xảy ra cao hơn dưới phân phối hậu nghiệm và loại bỏ các giá trị có khả năng xảy ra thấp hơn.
Tuy nhiên, khó khăn hiển nhiên là chúng ta thường không biết dạng hàm cụ thể của phân phối hậu nghiệm để làm căn cứ chấp nhận hay từ chối các giá trị đề xuất. Nếu biết rõ dạng hàm, chúng ta đã có thể tính toán trực tiếp các xác suất mà không cần tạo mẫu ngẫu nhiên. Giải pháp cho bài toán này chính là thuật toán Metropolis-Hastings.
Thuật toán Metropolis-Hastings
Thuật toán Metropolis-Hastings giúp quyết định giá trị đề xuất nào của theta sẽ được chấp nhận hoặc loại bỏ ngay cả khi không biết dạng hàm của phân phối hậu nghiệm. Hãy cùng phân tích thuật toán qua từng bước cụ thể.
Đầu tiên, chúng ta tính toán tỷ lệ r bằng mật độ hậu nghiệm tại theta mới chia cho mật độ hậu nghiệm tại theta cũ. Mặc dù không biết dạng hàm của phân phối hậu nghiệm, trong ví dụ này chúng ta có thể thay thế bằng tích của phân phối tiền nghiệm và hàm hợp lý.
Tiếp theo, chúng ta tính xác suất chấp nhận alpha, đây đơn giản là giá trị nhỏ nhất giữa tỷ lệ r và 1. Bước này là bắt buộc vì giá trị xác suất luôn phải nằm trong khoảng từ 0 đến 1. Nếu xác suất chấp nhận bằng 1, chúng ta chấp nhận giá trị theta mới và trung bình của phân phối đề xuất sẽ trở thành giá trị mới này ở vòng lặp tiếp theo.
Trong trường hợp xác suất chấp nhận nhỏ hơn 1, chúng ta không tự động từ chối giá trị mới. Thay vào đó, chúng ta rút một số ngẫu nhiên u từ phân phối đều trên khoảng từ 0 đến 1. Nếu u nhỏ hơn xác suất chấp nhận, chúng ta vẫn chấp nhận giá trị đề xuất mới. Ngược lại, chúng ta từ chối giá trị mới và giữ lại giá trị cũ từ bước trước.
Khi áp dụng thuật toán Metropolis-Hastings vào chuỗi MCMC với 10,000 vòng lặp, kết quả sẽ có sự thay đổi rõ rệt. Phân phối đề xuất thay đổi qua các vòng lặp, đồ thị đường vết không còn hiện tượng bước đi ngẫu nhiên mà thể hiện tính dừng ổn định, và đồ thị mật độ thu được đã định hình thành một phân phối có ích.
Đồ thị so sánh giữa biểu đồ tần suất của mẫu được tạo ra bằng thuật toán Metropolis-Hastings và phân phối hậu nghiệm lý thuyết Beta với các tham số 5 và 7 cho thấy sự trùng khớp rất cao. Đường biểu diễn mật độ nhân của mẫu bám sát đường phân phối lý thuyết, chứng minh rằng mẫu ngẫu nhiên thu được là một sự xấp xỉ xuất sắc cho phân phối hậu nghiệm thực tế. Từ mẫu này, chúng ta có thể dễ dàng tính toán giá trị trung bình, trung vị, khoảng tin cậy Bayes 95% hoặc xác suất theta rơi vào một khoảng bất kỳ.
Sử dụng MCMC và thuật toán Metropolis-Hastings trong Stata
Hãy quay trở lại với ví dụ tung đồng xu bằng cách sử dụng lệnh bayesmh trong Stata. Chúng ta thiết lập một phân phối tiền nghiệm beta với các tham số 1 và 1 cùng với hàm hợp lý bernoulli.
1bayesmh heads, likelihood(dbernoulli({theta})) prior({theta}, beta(1,1))
2Burn-in ...
3Simulation ...
4Model summary
5------------------------------------------------------------------------------
6Likelihood:
7 heads ~ bernoulli({theta})
8Prior:
9 {theta} ~ beta(1,1)
10------------------------------------------------------------------------------
11Bayesian Bernoulli model MCMC iterations = 12,500
12Random-walk Metropolis-Hastings sampling Burn-in = 2,500
13 MCMC sample size = 10,000
14 Number of obs = 10
15 Acceptance rate = .4454
16Log marginal likelihood = -7.7989401 Efficiency = .2391
17------------------------------------------------------------------------------
18 | Equal-tailed
19 | Mean Std. Dev. MCSE Median [95% Cred. Interval]
20-------------+----------------------------------------------------------------
21 theta | .4132299 .1370017 .002802 .4101121 .159595 .6818718
22------------------------------------------------------------------------------Kết quả đầu ra của Stata sẽ hiển thị chi tiết về quá trình mô phỏng. Dòng chữ burn-in cho biết giai đoạn khởi động (burn-in) gồm 2,500 vòng lặp đầu tiên đã bị loại bỏ để giảm thiểu ảnh hưởng của giá trị khởi đầu ngẫu nhiên của chuỗi. Kích thước mẫu MCMC chính thức dùng để phân tích là 10,000. Tỷ lệ chấp nhận thể hiện tỷ lệ các giá trị đề xuất của theta được giữ lại trong mẫu cuối cùng.
Bảng kết quả cũng cung cấp sai số tiêu chuẩn Monte Carlo, viết tắt là MCSE, giúp đánh giá độ chính xác của việc ước lượng giá trị trung bình hậu nghiệm thực tế.
Kiểm tra sự hội tụ của chuỗi
Khái niệm hội tụ trong ngữ cảnh MCMC khác biệt hoàn toàn so với trong ước lượng cực đại hàm hợp lý. Các thuật toán cực đại hàm hợp lý sẽ lặp cho đến khi tìm được điểm cực đại tối ưu. Trong khi đó, chuỗi MCMC chỉ đơn giản lặp cho đến khi đạt đủ kích thước mẫu yêu cầu rồi dừng lại. Việc chuỗi dừng chạy không đồng nghĩa với việc chúng ta đã thu được một mẫu tối ưu từ phân phối hậu nghiệm. Chúng ta cần kiểm tra chất lượng mẫu bằng các công cụ trực quan hóa thông qua lệnh bayesgraph diagnostics.
1bayesgraph diagnostics {theta}
Đồ thị chẩn đoán cung cấp bốn góc nhìn quan trọng bao gồm đồ thị đường vết, biểu đồ tần suất, đồ thị mật độ nhân và biểu đồ tự tương quan. Đồ thị đường vết thể hiện một mô hình dừng ổn định, đây là tín hiệu rất tốt. Biểu đồ tần suất không xuất hiện các hình thù dị thường hay hiện tượng đa cực trị. Các đường mật độ nhân của toàn bộ mẫu, nửa đầu chuỗi và nửa sau chuỗi gần như trùng khớp nhau, cho thấy tính ổn định cao của chuỗi theo thời gian. Cuối cùng, biểu đồ tự tương quan giảm nhanh chóng về mức không khi độ trễ tăng lên. Tất cả các dấu hiệu này đều xác nhận chuỗi đã hội tụ tốt và mẫu thu được hoàn toàn đáng tin cậy.
✨ Giá trị đắt giá:
Bản chất cốt lõi của thuật toán Metropolis-Hastings là khả năng khám phá phân phối hậu nghiệm phức tạp mà không cần biết hằng số chuẩn hóa của nó. Bằng cách so sánh tỷ lệ giữa hai điểm dữ liệu, chúng ta triệt tiêu được phần mẫu số khó tính toán nhất trong công thức Bayes. Điều này biến một bài toán tích phân đa chiều bất khả thi thành một chuỗi các quyết định so sánh đơn giản, mở đường cho việc ứng dụng thống kê Bayes vào mọi lĩnh vực khoa học phức tạp.
Câu hỏi tư duy hoặc bài tập ứng dụng:
Giả sử bạn đang chạy một chuỗi MCMC trong Stata và nhận thấy tỷ lệ chấp nhận trong kết quả đầu ra rất thấp, chỉ khoảng dưới 5%. Điều này phản ánh vấn đề gì trong việc lựa chọn phân phối đề xuất, và bạn sẽ điều chỉnh phương sai của phân phối đề xuất như thế nào để cải thiện hiệu suất của thuật toán?
