Giới thiệu về Phương pháp Monte Carlo

Giới thiệu

Hai nhóm bài toán số thường gặp trong quy trình phân tích dữ liệu là bài toán tối ưu hóa và bài toán tích phân. Không phải lúc nào cũng có thể tính giải tích các ước lượng gắn với một mô hình nhất định, nên chúng ta thường phải dùng đến lời giải số. Một cách để tránh vấn đề đó là sử dụng mô phỏng. Ước lượng Monte Carlo đề cập đến việc mô phỏng các mẫu giả định từ một phân phối xác suất nhằm tính các đại lượng quan trọng của phân phối đó.

Ý tưởng cơ bản của Monte Carlo là viết tích phân dưới dạng kỳ vọng theo một phân phối xác suất nào đó, rồi xấp xỉ bằng ước lượng phương pháp moment ($E[g(X)] \approx \overline{g(X)} = \dfrac{1}{n}\sum g(X_{i})$).

Nếu ta có một hàm liên tục $g(\theta)$ và muốn lấy tích phân trên khoảng (a,b), ta có thể viết lại tích phân như một kỳ vọng theo phân phối đều $U \sim U[a,b]$, cụ thể:

Sử dụng ước lượng phương pháp moment, phép xấp xỉ tích phân là:

Trong đó

Là các giá trị được mô phỏng từ một phân phối đều.

Ví dụ 1: Xấp xỉ tích phân hàm mũ

1) Với hàm $f(x)= e^{x}$, tích phân trên khoảng [3,5] là:

Xấp xỉ Monte Carlo của tích phân là:

#       Declaring the desired function
f = function(x){return(exp(x))}
#       Declaring the absolute error function
error = function(x,y){return(abs(x-y))}
#       The actual integral answer
ans = exp(5)-exp(3)

set.seed(6971)
#       number of iterations
n = 10^2
#       simulated uniform data
x= runif(n,3,5)
#       MonteCarlo approximation
MCa= (5-3)*mean(f(x))
#       Approximation error
e = error(ans,MCa)

rest =  data.frame(n = n,MCapprox = MCa,error = e)
set.seed(6971)
for(k in 3:6){
  n = 10^k
  x = runif(n,3,5)
  mca = (5-3)*mean(f(x))
  rest = rbind(rest,c(n,mca,error(ans,mca) ) )
}

kable(rest,digits = 5,align = 'c',caption = "Integral Monte Carlo approximation results",
      col.names =c("Number of simulations","Monte Carlo approximation","Error approximation"))Được hỗ trợ bởi

Xấp xỉ Monte Carlo tổng quát

Trong trường hợp tổng quát, xấp xỉ tích phân cho một phân phối f cho trước là:

Một thuật toán để xây dựng $\widehat{I}$ có thể mô tả theo các bước sau:

1) Sinh từ một phân phối f

2) Tính:

3) Lấy trung bình mẫu: $$\overline{I} = \dfrac{1}{n}\sum_{k=1}^{n}\dfrac{g(\theta_k)}{f(\theta_k)}$$

Trong đoạn tiếp theo, hàm xấp xỉ Monte Carlo đơn giản được trình bày để minh họa cách thuật toán hoạt động, trong đó a và b là các tham số của mật độ đều, n là số lần mô phỏng mong muốn, và f là hàm cần tích phân.

# The simple Monte Carlo function
MCaf = function(n,a,b,f){
x = runif(n,a,b)
MCa = (b-a)*mean(f(x))
return(MCa)
}Được hỗ trợ bởi

Phương pháp Monte Carlo trong phân tích dữ liệu Bayes

Ý tưởng chính của phân tích dữ liệu Bayes là khớp một mô hình ( như mô hình hồi quy hoặc chuỗi thời gian ) theo cách tiếp cận suy luận Bayes. Ta giả định các tham số quan tâm có một phân phối lý thuyết; phân phối này (hậu nghiệm) được cập nhật bằng phân phối của dữ liệu quan sát (hàm hợp lý) và thông tin trước đó hoặc bên ngoài về các tham số (phân phối tiên nghiệm) thông qua định lý Bayes.

$$p(\theta / X) \text{ } \alpha\text{ } p(X/\theta)p(\theta)$$ Trong đó:

$p(\theta/X)$ là phân phối hậu nghiệm của tham số

$p(X/\theta)$ là phân phối lấy mẫu của dữ liệu quan sát ( hàm hợp lý ).

$p(\theta)$ là phân phối tiên nghiệm của tham số.

Vấn đề chính trong cách tiếp cận Bayes là ước lượng phân phối hậu nghiệm. Các phương pháp Markov Chain Monte Carlo ( mcmc ) sinh mẫu từ phân phối hậu nghiệm và xấp xỉ các kỳ vọng, xác suất hoặc lượng phân vị bằng phương pháp Monte Carlo.

Trong hai phần tiếp theo, chúng tôi đưa ra hai ví dụ để xấp xỉ xác suất và lượng phân vị của một phân phối lý thuyết. Quy trình trình bày ở trên là những phương pháp thường dùng trong cách tiếp cận Bayes.

Ví dụ 2: Xấp xỉ xác suất của phân phối gamma

Giả sử ta muốn tính xác suất biến ngẫu nhiên $\theta$ nằm giữa 0 và 5, $P(0 < \theta < 5)$, với $\theta$ có phân phối gamma với tham số a = 2 và b = 1/3 ($\theta \sim Gamma(a = 2,b = 1/3)$), khi đó xác suất là:

Với $I_{[0,5]}(\theta) = 1$ nếu $\theta$ thuộc khoảng [0,5]. Ý tưởng của xấp xỉ Monte Carlo là đếm số quan sát thuộc khoảng [0,5], rồi chia cho tổng số dữ liệu mô phỏng.

set.seed(6972)
#       number of iterations
n = 10^2
#       simulated uniform data
x= rgamma(n,shape = 2,1/3)
#       MonteCarlo approximation
MCa= mean(x <= 5)
#       Approximation error
e = error(pgamma(5,2,1/3),MCa)

rest =  data.frame(n = n,MCapprox = MCa,error = e)

for(k in 3:6){
  n = 10^k
  x= rgamma(n,shape = 2,1/3)
  mca= mean(x <= 5)
 rest = rbind(rest,c(n,mca,error(pgamma(5,2,1/3),mca) ) )
}

kable(rest,digits = 5,align = 'c',caption = "Probability Monte Carlo approximation results",
      col.names =c("Number of simulations","Monte Carlo approximation","Error approximation"))Được hỗ trợ bởi

Ví dụ 3: Xấp xỉ lượng phân vị của phân phối chuẩn

Giả sử ta muốn tính lượng phân vị 0,95 của biến ngẫu nhiên $\theta$ có phân phối chuẩn với các tham số $\mu = 20$ và $\sigma = 3$ ($\theta \sim normal(\mu = 20,\sigma^{2} = 9)$), khi đó lượng phân vị 0,95 là:

Ý tưởng chính là tìm giá trị mẫu lớn nhất cho xác suất bằng hoặc nhỏ hơn 0,95; xấp xỉ lượng phân vị Monte Carlo được ước lượng bằng hàm quantile() áp dụng cho dữ liệu mô phỏng.

set.seed(6973)
#       number of iterations
n = 10^2
#       simulated uniform data
  x= rnorm(n,20,3)
#       MonteCarlo approximation
MCa= quantile(x,0.95)
#       Approximation error
e = error(qnorm(0.95,20,3),MCa)

rest =  data.frame(n = n,MCapprox = MCa,error = e)

for(k in 3:6){
  n = 10^k
  x= rnorm(n,20,3)
  mca= quantile(x,0.95)
 rest = rbind(rest,c(n,mca,error(qnorm(0.95,20,3),mca) ) )
}

kable(rest,digits = 5,align = 'c',caption = "Quantile Monte Carlo approximation results",
      col.names =c("Number of simulations","Monte Carlo approximation","Error approximation"),row.names = FALSE)Được hỗ trợ bởi

Thảo luận và Kết luận

Các phương pháp xấp xỉ Monte Carlo cung cấp một công cụ thay thế để xấp xỉ tích phân và là công cụ thiết yếu trong suy luận Bayes, đặc biệt khi làm việc với các mô hình tinh vi và phức tạp. Như thấy ở cả ba ví dụ, phương pháp Monte Carlo cho độ xấp xỉ rất tốt, nhưng đòi hỏi số lượng mô phỏng rất lớn để sai số xấp xỉ tiến gần về 0.

Tài liệu tham khảo

1) Introducing Monte Carlo methods with R, Springer 2004, Christian P. Robert and George Casella.

2) Handbook of Markov Chain Monte Carlo, Chapman and Hall, Steve Brooks, Andrew Gelman, Galin L. Jones, and Xiao-Li Meng.

3) Introduction to mathematical Statistics, Pearson, Robert V. Hogg, Joseph W. Mckean, and Allen T. Craig.

4) Statistical Inference An Integrated Approach, Chapman and Hall, Helio S. Migon, Dani Gamerman, Francisco Louzada.