Tuần 1.5. Các thuật toán Markov Chain Monte Carlo

Seminar ngày 13/4 tạm nghỉ vì có lễ ra mắt phòng Data Science bên VIASM. Hôm nay (20/4) mình lại có cơ hội được trình bày tại phòng bộ môn xác suất thống kê. Hai chủ đề chính của ngày hôm nay là về "các thuật toán sinh mẫu Markov Chain Monte Carlo" (MCMC) và mở đầu cho "Thống kê Bayes phi tham số"-phần chính của chuỗi seminar này. Trong bài viết này mình sẽ tóm tắt lại về các thuật toán sinh mẫu MCMC.

Monte Carlo là một thành phố lớn ở Monaco-đất nước giàu nhất thế giới, nổi tiếng với các sòng bạc và các vụ cá cược. Phương pháp Monte Carlo được đặt tên dựa vào sự kiện rằng về lâu về dài, tần suất (thực nghiệm) xảy ra các biến cố sẽ xấp xỉ xác suất lí thuyết (giả dụ như trong roulette, xúc xắc, cò quay Nga...). Cụ thể hơn ta có

Định lí (về tính ergodic). $x^{(1)},x^{(2)},\dots$ là mẫu được sinh ngẫu nhiên từ $X$ và $A$ là một tập con của $\mathbb{R}$ thì

$$ \lim_{n\to \infty} \dfrac{ \# \{x(i) \in A \} }{n} = P(X\in A).$$

Vì thế khi ta đã biết hàm mật độ $f_X$ của $X$ nhưng rất khó để tính tích phân từ $f_X$, ta sẽ dùng phần mềm thống kê (R) để sinh mẫu ngẫu nhiên từ $X$, sau đó sử dụng xấp xỉ trên để tính số xác suất của $X$. Phương pháp đó được gọi là Monte Carlo.

Quay trở lại với Markov Chain Monte Carlo: Là một thuật toán quan trọng trong thống kê Bayes vì nó cho phép ta tính toán xác suất của phân bố hậu nghiệm (cái mà thường rất khó để tính hiển ra trong trường hợp nhiều tham số). 3 phương pháp được nhắc tới tại đây là Gibbs, Metropolis và Metropolis-Hasting.

1. Thuật toán Gibbs.

Thuật toán này thường được dùng với phân bố nhiều chiều, khi mà phân phối đồng thời thì khó tính nhưng phân phối có điều kiện thì có thể tính được dễ dàng. Gọi $\phi = (\phi_1,\dots, \phi_n)$, giả sử đã biết phân phối có điều kiện $p(\phi_i |\phi_1,\dots, \phi_{i-1},\phi_{i+1},\dots, \phi_n)$  

Bước 1. Lấy $\phi^{(0)} = (\phi_1^{(0)},\dots, \phi_s^{(0)})$.

Bước k.($\forall \, k$)  ta lấy mẫu dựa trên các phân phối có điều kiện (và update thông tin liên tục qua từng bước lấy mẫu)

$\phi_1^{(k)} ~ p(\phi_1|\phi_2^{(k-1)},\phi_2^{(k-1)},\dots, \phi_n^{(k-1)}$;

$\phi_2^{(k)} ~ p(\phi_2|\phi_1^{(k)},\phi_3^{(k-1)},\dots, \phi_n^{(k-1)}$;

....

$\phi_n^{(k)} ~ p(\phi_n|\phi_1^{(k)},\phi_2^{(k)},\dots, \phi_{n-1}^{(k)}$;

Cuối cùng đặt $\phi^{(k)} = (\phi_1^{(k)},\dots, \phi_s^{(k)})$, ta có mẫu thứ $k$.

Chú ý rằng việc lấy mẫu thứ $k$ chỉ phụ thuộc vào mẫu thứ $k-1$, vì vậy thuật toán này sinh ra một xích Markov (đó là lí do nó được gọi là MCMC).

Ví dụ: (Phân bố hậu nghiệm của mô hình có phân bố chuẩn $N(\theta,\sigma^2)$)

Từ bài viết tuần 1, ta đã biết rằng để có tính (nửa) liên hợp thì phân bố của $\theta ~ N(\mu_0,\tau_0^2), \tilde{\sigma}^2~  Gamma(\nu_0/2, \nu_0 \sigma_0^2/2)$, phân bố hậu nghiệm được tính trong bài viết tuần 1. Rất khó để tìm phân bố đồng thời của $\theta, \sigma$ nhưng ta có thể tính phân bố điều kiện hậu nghiệm $\tilde{\sigma}|\theta, y_1,\dots, y_n$ có dạng Gamma với tham số khá đơn giản, từ đó triển khai thuật toán lấy mẫu Gibbs.

2. Thuật toán Metropolis.

Khi có thể tính được phân bố đồng thời, ta sử dụng Monte Carlo bình thường, khi tính được các phân bố điều kiên, ta sử dụng Gibbs, còn khi chẳng tính được cái gì, ta đã có thuật toán Metropolis.

Giả sử đã sinh được $\phi^{(0)},\phi^{(1)},\dots, \phi^{(k)}$. Sinh ngẫu nhiên $\phi^{*}$ gần $\phi^{(k)}$ (từ $J(\theta^{*}|\theta^{(s)} = U(\phi^{(k)}-\delta, \phi^{(k)}+\delta)$ hoặc $N(\phi^{(k)},\delta^2)$ với $\delta$ đủ nhỏ chẳng hạn, gọi là phân bố đề nghị). Xét tỉ số "chấp thuận" (sử dụng công thức Bayes)

$$r = \dfrac{p(\theta^{*}|y)}{p(\theta^{(k)}|y)} = \dfrac{p(y|p(\theta^{*}) p(\theta^{*})}{p(y|p(\theta^{(s)}) p(\theta^{(s)})}.$$

Nếu $r>1$, ta có $\theta= \theta^{*}$ dễ xảy ra hơn $\theta= \theta^{(k)}$, chọn luôn $\theta^{(k+1)}=\theta^* $. 

Nếu $r\leq 1$, sinh ngẫu nhiên $u ~ U[0,1]$, chọn $\theta^{(k+1)}=\theta^* $ nếu $r>u$ và $=\theta^{(k)}$ nếu ngược lại.

3. Thuật toán Metropolis-Hasting

Thuật toán này là tổng quát hóa của 2 thuật toán trên. 

Để cho tiện, ta xét biến ngẫu nhiên 2 chiều $(u,v)$, gọi $p_0(u,v) = p(u,v|y)$. Chọn một phân bố đề nghị $J(u^*,v^*|u,v)$ ở bước thứ k, ta chia thành 2 bước nhỏ:

   (1) Cập nhật $u$

       (a) Lấy mẫu $u^* ~ J_u (u|u^{(k)}, v^{(k)})$ (Phân bố biên duyên của $U$);

       (b) Tính $r = \dfrac{p(u^*|v^{(k)}}{p(u^{(k)}|v^{(k)})} \dfrac{J_u (u^{(k)}|u^{*}, v^{(k)})}{J_u (u^{*}|u^{(k)}, v^{(k)})}$;

       (c) Chọn $u^{k+1} = u^*$ với xác suất $\min (1,r)$ và $= u^{(k)}$ với xác suất $\max (0,1-r)$.

  (2) Cập nhật $v$ 

       (a) Lấy mẫu $v^* ~ J_v (v|u^{(k+1)}, v^{(k)})$ (Phân bố biên duyên của $V$, $u$ đã được cập nhật);

       (b) Tính $r = \dfrac{p(v^*|u^{(k)}}{p(v^{(k)}|u^{(k)})} \dfrac{J_v (v^{(k)}|u^{(k+1)}, v^{*})}{J_v (v^{*}|u^{(k+1)}, v^{(k)})}$;

       (c) Chọn $u^{k+1} = u^*$ với xác suất $\min (1,r)$ và $= u^{(k)}$ với xác suất $\max (0,1-r)$.

Dễ thấy thuật toán Metropolis-Hasting là tổng quát hóa của 2 thuật toán trên vì khi cho $J$ là đối xứng $J_u (u^{(k)}|u^{*}, v^{(k)})= J_u (u^{*}|u^{(k)}, v^{(k)}), J_v (v^{(k)}|u^{(k+1)}, v^{*})= J_v (v^{*}|u^{(k+1)}, v^{(k)})$ thì ta nhận được Gibbs. Còn khi  phân bố đề nghị bằng phân bố điều kiên: $J_u(u^*|u^{(k)},v^{(k)}) = p(u^*|v^{(k)})$ ta nhận được thuật toán Metropolis.

Với các mô hình quen thuộc, xích Markov trong bài luôn có tính bất khả quy, aperiodic và recurrent, vậy nó đều tồn tại duy nhất phân bố dừng, và phân bố đó chính là phân bố hậu nghiệm ta cần tìm.