當我們遇到比較複雜的機率分配，難以直接抽樣時，可以利用以下兩種方法處理，分別為拒絕抽樣（Rejection Sampling）和重要性抽樣（Importance Sampling）。

1 拒絕抽樣

1.1 從矩形中抽樣

假設一機率分配\(g(x)\)，不屬於任何一種常見的機率分配，其中變數值\(x\)的範圍有限。我們可利用一個矩形包覆住機率分配其密度不等於\(0\)的範圍，從矩形中均勻地抽出點\((U_1,U_2)\)，再篩選出落在機率分配內的點，這些點的橫軸座標\(U_1\)即為來自該機率分配的亂數。此方法適用在隨機變數值是有限範圍的情況。

上圖三角形為一特殊的機率分配，變數值範圍介於0到2。為了生成此分配的亂數， 我們可以從一個包覆著分配的矩形中均勻的抽樣，再篩選出落在三角形內的樣本。

1.2 從機率分配抽樣

假設我們希望從機率分配\(g(x)\)抽樣。可利用機率分配\(f(x)\)，滿足\(kf(x) \geq g(x), \; \forall \;x\)。從\(kf(x)\)中抽出樣本點\((X_1,U_2)\)，其中\(X_1 \rightarrow f(x))\)且\(U_2 \rightarrow U(0,k)\)，篩選出落在\(g(x)\)內的樣本點，這些點的橫座標\(F_1\)即為來自\(g(x)\)的亂數。\(k\)應該盡可能選用較小的值，使得\(kf(x)\)更貼近\(g(x)\)，增加樣本點落在\(g(x)\)內的機率。

1.3 練習

1.3.1

假設我們希望從以下機率分配\(g(x)\)抽樣： \[ g(x) = \left\{\begin{array}{ll} 1- |1-x|, & if \;\; 0\leq x\leq 2 \\ 0, & \mbox{otherwise} \end{array} \right. \] 觀察可知\(g(x)\)在\([0,2]\)的範圍內，因此可用一個長度為1、寬度為2的矩形包覆\(g(x)\)。首先從矩形均勻地抽出樣本點\((U_1,U_2)\)，其中\(U_1 \rightarrow U(0,2)\)且\(U_2 \rightarrow U(0,1)\)，再挑選出落在\(g(x)\)底下的點。

kN <- 1000000
kU1 <- runif(kN, 0, 2)
kU2 <- runif(kN, 0, 1)

# 篩選出落於g(x)下的點
r.gx <- kU1[kU2 < (1 - abs(1 - kU1))]

# 以直方圖觀察樣本的分配
hist(r.gx, freq = F, breaks = "Scott")

1.3.2

使用拒絕抽樣方法從標準常態分配的區間[-2,2]內內抽樣

kN <- 1000000


# 從矩形中均勻的抽樣
r.unif.x <- runif(kN, min=-2, max=2)  # 橫座標
r.unif.y <- runif(kN, min=0, max=1)  # 縱座標


# 挑選落在標準常態分配下的亂數
r.norm.x <- r.unif.x[r.unif.y < dnorm(r.unif.x)]

# 以直方圖觀察亂數的分配
hist(r.norm.x, freq=F, breaks="Scott", xlim=c(-3,3))

1.3.3

假設有一機率分配\(Cg(x) = Ce^{-x^{1.5}} \; , \forall \; x>0\)，\(C\)是未知常數。已知\(e^{-x^{1.5}} \leq 2e^{-x}\)，可令\(f(x) = e^{-x}\)且$ k = 2\(。從\)kf(x)抽出點(X_1,U_2)\(，再挑出落在\)g(x)$底下的點。

set.seed(7)
kN <- 1000000
kK <- 2

# 設置g(x)
f.gx <- function(x) exp(-x^(1.5))

# 生成來自kf(x)的點
r.exp <- rexp(kN, rate = 1)
r.unif <- runif(kN)
r.kfx <- kK * r.unif * dexp(r.exp)

# 取落在g(x)/C面積下的點
r.gx <- r.exp[r.kfx  < f.gx(r.exp)]


# 繪圖：直方圖為Cg(x)的亂數，虛線為kf(x)函數
hist(r.gx, freq = F, breaks = "Scott")  # breaks 指定bin個數或計算方式
f.fx <- function(x) kK * dexp(x)
curve(f.fx, from = 0, to = max(r.gx), add = TRUE, lty = 2)

2 重要性抽樣

重要性抽樣是一蒙地卡羅積分的一個延伸方法，藉由使用便於抽樣的分配取代原本的分配生成亂數，以估計原分配的參數（例如期望值）。

2.1 方法

\[\mathbb{E}_f(h(X)) = \int_x h(x)f(x)dx = \int_x [h(x) \frac{f(x)}{g(x)}] g(x) dx = \mathbb{E}_g(h(x) \frac{f(x)}{g(x)}) \]

其中\(g(x)\)的支撐集（support）需包含隨機變數函數在原分配下的支撐集 \[supp(g) ⊃ supp(h × f)\] \[\frac{1}{n} \sum^n_{j=1} h(x_j) \frac{f(x_j)}{g(x_j)} \rightarrow \mathbb{E}_f(h(X)) \]

重要性抽樣的優點

  1. 可以使用便於抽樣的分配代替複雜的原分配生成樣本
  2. 藉由權重調整每個樣本點對估計的影響，可以保留所有樣本，比拒絕抽樣有效率，對參數的估計也能更精確
  3. 計算平均值時，若除以權重總和會調整樣本的尺度，不用求出原分配的未知常數

上圖紅線代表任意一種我們估計期望值的分配，而藍線則是為了便於抽樣而引進的分配，其生成的亂數為圖中藍點，亂數乘上權重後轉變為紅點。雖然加權後的亂數分布不會符合原本的分配，但能夠正確的估計出原分配的期望值，如黑色虛線所示。

2.2 步驟

  1. 從一個便於抽樣的機率分配\(g(x)\)抽樣
  2. 計算每個樣本點的權重\(f(x)/g(x)\)
  3. 計算樣本點乘上權重後的平均值

2.3 練習

2.3.1

估計機率分配\(g(x)\)的期望值和變異數 \[g(x) = Ce^{-x^{1.5}}, \; x>0\]

kN <- 10000000

f.gx <- function(x) exp(-x^(1.5))
r.x <- rexp(kN)

# 計算亂數的權重
r.w <- f.gx(r.x) / dexp(r.x)

# 估計g(x)期望值
k.Avg <- sum(r.x * r.w) / sum(r.w) 
k.Avg
## [1] 0.6592922

# 估計g(x)變異數
k.Var <- (r.x - k.Avg)^2
sum(k.Var * r.w) / sum(r.w)
## [1] 0.3033673

2.3.2

假設\(X \rightarrow N(0,1)\)，欲估計\(P(X>4.5)\)。由於\(X>4.5\)的機率非常低，若直接從常態分配中抽樣，需要非常大量的樣本才能精確的估計出機率值。因此我們利用重要性抽樣方法，引入位移的指數分配取代原始分配做抽樣。由下列推導可知，我們可用權重的平均值估計\(P(X>4.5)\)。

\[ \begin{align} P(X>4.5) &= \int^{\infty}_{4.5}f(x)dx\\ &= \int_xI_{4.5,\infty}(x)f(x)dx \\ &= \int_xI_{4.5,\infty}(x) \frac{f(x)}{g(x)}g(x)dx \\ &= \mathbb{E}_g(\frac{f(x)}{g(x)}) \end{align} \]

\[ \begin{align} \frac{1}{n} \sum^n_{i = 1} \frac{f(x_i)}{g(x_i)} \rightarrow \mathbb{E}_g(\frac{f(x)}{g(x)}) = P(X>4.5) \end{align} \]

set.seed(7)
kN <- 10000000
r.norm <- rnorm(kN)
# 從位移指數分配抽樣
r.exp <- rexp(kN) + 4.5


# 計算權重
r.weight <- dnorm(r.exp) / dexp(r.exp - 4.5)

# 計算平均值
sum(1 * r.weight) / kN
## [1] 3.396683e-06

par(mfrow=c(1,2))
# 逐一增加樣本的數量並計算平均值，觀察收斂速度
plot( cumsum(r.weight)/ (1 : kN),
      type="l",
      xlim = c(0,1000),
      ylim = c(0,10^(-5)),
      xlab = 'numbers of sample',
      ylab = 'prob',
      main = 'Importance Sampling')
abline(a=pnorm(-4.5),b=0,col="red")  # P(X>4.5)理論值線

# 使用常態分配直接抽樣，觀察收斂速度
plot(cumsum(r.norm > 4.5) / (1 :kN),
     type="l",
     xlim = c(0,10000000),
     ylim = c(0,10^(-5)),
     xlab = 'numbers of sample',
     ylab = 'prob',
     main = 'Sampling from N(0,1)')
abline(a=pnorm(-4.5),b=0,col="red")

從圖左可看出使用重要性抽樣估計\(P(X>4.5)\)，樣本數到400後就已經接近理論值\(3.397673 \times 10^{-06}\)；圖右則顯示直接從標準常態分配抽樣，需要將近1千萬個樣本估計值才能貼近理論值。

本文是蒙地卡羅方法的筆記與實作，實作的題目出自以下參考來源：

"A First Course in Statistical Programming with R"
"Monte Carlo Methods with R"