12.2 帶有舍選控制的重要抽樣法

在重要抽樣法和標準化重要抽樣法的實際應用中，

好的試抽樣分布很難獲得，

所以權重\(\{ W_i = f(\boldsymbol X_i)/g(\boldsymbol X_i) \}\)經常會差別很大，

使得抽樣樣本主要集中在少數幾個權重最大的樣本點上。

為此，可以舍棄權重太小的樣本點，

重新抽樣替換這樣的樣本點,

這種方法稱為帶有舍選控制的重要抽樣法。

需要預先選定權重的一個閾值\(c\)，并計算

\[\begin{align*}

p_c = \int \min\left\{1, \frac{f(\boldsymbol x)}{c g(\boldsymbol x)} \right\} g(\boldsymbol x) \,d\boldsymbol x

= \int \min\left\{g(\boldsymbol x), \frac{f(\boldsymbol x)}{c} \right\} \,d\boldsymbol x .

\end{align*}\]

在產生每個抽樣點\(\boldsymbol X_i\)時，計算權重\(W_i\)，

當權重\(W_i \leq c\)時接受此抽樣點，但權重改為\(W_i^* = p_c W_i\)；

當\(W_i < c\)時僅以\(W_i/c\)的概率接受此抽樣點，并調整權重為\(p_c c\)，

如果沒有被接受則重新從\(g(\cdot)\)中抽取。

算法如下：

帶有舍選控制的重要抽樣法

for(\(i\) in \(1:N\)) {

\(\quad\) repeat {

\(\qquad\) 從\(g(\cdot)\)抽取\(\boldsymbol\xi_i\)

\(\qquad\) 計算\(W_i \leftarrow f(\boldsymbol\xi_i)/g(\boldsymbol\xi_i)\)

\(\qquad\) if(\(W_i \geq c\)) {

\(\qquad\quad\) 令\(\boldsymbol X_i \leftarrow \boldsymbol\xi_i\), \(W_i^* \leftarrow p_c W_i\)

\(\qquad\quad\) break

\(\qquad\) } else {

\(\qquad\quad\) 從U(0,1)抽取\(U\)

\(\qquad\quad\) if(\(U < W_i / c\)) {

\(\qquad\qquad\) 令\(\boldsymbol X_i \leftarrow \boldsymbol\xi_i\), \(W_i^* \leftarrow p_c c\)

\(\qquad\qquad\) break

\(\qquad\quad\) }

\(\qquad\) }

\(\quad\) } # repeat

} # for

輸出\((\boldsymbol X_i^*, W_i^*)\), \(i=1,2,\dots,N\)

在實際應用帶有舍選控制的重要抽樣法時，

一般先選略小的抽樣量\(N_0\)進行原始的重要抽樣，

分析權重\(W_i\)的分布，

據此選擇權重\(c\)，

比如可以選\(\{ W_i \}\)的某個分位數。

舍選控制算法中的\(p_c\)是歸一化常數，

如果使用標準化重要抽樣法，

\(p_c\)可以省略，

否則，

\(p_c\)可以從原始的重要抽樣法結果中估計為

\[\begin{align}

\hat p_c = \frac{1}{N_0} \sum_{i=1}^{N_0} \min\left\{1, \; \frac{W_i}{c} \right\}.

\end{align}\]

帶有舍選控制的重要抽樣法得到的\(\{(\boldsymbol X_i, W_i^*), \ i=1,\dots,N \}\)是關于\(f(\cdot)\)適當加權的,

被接受的\(\boldsymbol X_i^*\)的分布密度不再是試投密度\(g(\boldsymbol x)\)，而是改成了

\[\begin{align}

g^*(\boldsymbol x) = \frac{1}{p_c} \min\left\{g(\boldsymbol x), \frac{f(\boldsymbol x)}{c} \right\}.

(\#sim-rq-rjgstar)

\end{align}\]

例12.1用MC積分法計算

\(I = \int_0^1 e^x \,dx = e-1 \approx 1.718\)。

對被積函數\(h(x)=e^x\)做泰勒展開得

\[\begin{aligned}

e^x = 1 + x + \frac{x^2}{2!} + \cdots

\end{aligned}\]

取

\[\begin{aligned}

g(x)=c(1+x) = \frac{2}{3}(1+x)

\end{aligned}\]

要產生\(g(x)\)的隨機數可以用逆變換法，

密度\(g(x)\)的分布函數\(G(x)\)的反函數為

\[\begin{aligned}

G^{-1}(y) = \sqrt{1+3y}-1, \ 0

\end{aligned}\]

因此，取\(U_i\) iid U(0,1),

令\(X_i = \sqrt{1+3U_i}-1\), \(i=1,2,\ldots,N\),

則重要抽樣法的積分公式為

\[\begin{aligned}

\hat I_3 = \frac{1}{N}

\sum_{i=1}^N \frac{e^{X_i}}{\frac{2}{3}(1+X_i)}

\end{aligned}\]

漸近方差為

\[\begin{aligned}

\text{Var}(\hat I_3) = \frac{1}{N} \left(

\frac{3}{2}\int_0^1 \frac{e^{2x}}{1+x} dx - I^2 \right)

\approx 0.02691/N.

\end{aligned}\]

積分真值：

重要抽樣法的程序實現：

一次模擬的估計值和絕對誤差、相對誤差：

相對誤差很小，說明有效位數在三位以上。

用一次模擬的結果估計平均相對誤差為：

說明估計精度大約有三位有效數字。

從一次模擬中對\(\text{Var}(\hat I_3)*N\)進行估計：

理論值是0.02691。

如果用平均值法， 估計公式為

\[\begin{aligned}

\hat I_2 = \frac{1}{N} \sum_{i=1}^N e^{U_i},

\end{aligned}\]

漸近方差為

\[\begin{aligned}

\text{Var}(\hat I_2) = \frac{1}{N} \left(\int_0^1 e^{2x} \,dx - I^2 \right)

\approx 0.2420/N \approx 9.0 \times \text{Var}(\hat I_3)

\end{aligned}\]

是重要抽樣法方差的9倍。

因為隨機模擬誤差的方差通常是\(\text{Var}(Y_i)/N\)，

其中\(Y_i\)是重復\(N\)次抽樣中的一次的估計，

所以方差的倍數可以換算成計算量的倍數。

在這個問題中，

為了達到相同的精度，

重要抽樣法只需要平均值法的\(1/9\)的樣本量。

平均值法的R程序：

一次模擬的估計值和絕對誤差、相對誤差：

相對誤差很小，說明有效位數約有三位。

用一次模擬的結果估計平均相對誤差為：

說明估計精度大約有二位有效數字。

從一次模擬中對\(\text{Var}(\hat I_2)*N\)進行估計：

理論值是0.2420。

如果用隨機投點法，\(h(x)=e^x \leq e(0 < x < 1)\), 取上界\(M=e\),

向\([0,1] \times [0,M]\)隨機投點，落到\(f(x)\)下方的概率為

\[\begin{aligned}

p = I/(M(b-a)) = (e-1)/e,

\end{aligned}\]

設投\(N\)點落到\(h(x)\)下方的頻率為\(\hat p\), 用隨機投點法估計\(I\)的公式為

\[\begin{aligned}

\hat I_1 = \hat p \cdot M (b-a) = e \hat p,

\end{aligned}\]

漸近方差為

\[\begin{aligned}

\text{Var}(\hat I_1) = e^2 p(1-p) / N = (e-1)/N \approx 1.718 / N

\approx 7.1 \times \text{Var}(\hat I_2)

\approx 64.8 \times \text{Var}(\hat I_3)

\end{aligned}\]

即重要抽樣法只需要隨機投點法的\(1/65\)的樣本量。

可見選擇合適的抽樣算法對減少計算量、提高精度是十分重要的。

隨機投點法的R程序:

一次模擬的估計值和絕對誤差、相對誤差：

相對誤差很小，說明有效位數將近三位。

用一次模擬的結果估計平均相對誤差為：

說明估計精度大約有二位有效數字。

從一次模擬中對\(\text{Var}(\hat I_1)*N\)進行估計：

理論值是1.718。

※※※※※

例12.2

設二元函數\(f(x,y)\)定義如下

\[\begin{aligned}

f(x,y) =& \exp\{ -45(x + 0.4)^2 - 60(y-0.5)^2 \}

+ 0.5 \exp\{ -90(x-0.5)^2 - 45(y+0.1)^4 \}

\end{aligned}\]

求如下二重定積分

\[\begin{aligned}

I =& \int_{-1}^1 \int_{-1}^1 f(x,y) \,dx\,dy

\end{aligned}\]

\(f(x,y)\)有兩個分別以\((-0.4, 0.5)\)和\((0.5, -0.1)\)為中心的峰，

對積分有貢獻的區域主要集中在\((-0.4,0.5)\)和\((0.5, -0.1)\)附近，

在其他地方函數值很小，對積分貢獻很小。

\(f(x,y)\)寫成R函數：

用平均值法， 取點數\(N=10000\)。

一次模擬的估計值為：

估計的平均相對誤差為:

\(\hat I_2\)的一個估計值為\(\hat I_2 = 0.126\),

從這一次模擬估計的平均相對誤差為0.03,

僅有一位有效數字精度，

只能保證\(\hat I_2=0.1\)基本可信。

平均值法的缺點是在整個正方形區域\((-1,1) \times (-1,1)\)上均勻布點，

而被積函數\(f(x,y)\)僅在兩個峰附近有較大正值，

其它區域基本是零。

用重要抽樣法改進。取試投密度為

\[\begin{aligned}

g(x,y) \propto& \tilde g(x,y) \\

=& \exp \{ -45(x+0.4)^2 - 60(y-0.5)^2 \}

+ 0.5 \exp\{ -90(x-0.5)^2 - 10(y+0.1)^2 \}, \\

& -\infty < x < \infty, -\infty < y < \infty,

\end{aligned}\]

這樣抽取到\([-1,1]\times [-1,1]\)范圍外的點對積分沒有貢獻，

因為構成\(g(x,y)\)的兩個密度都很集中，所以效率損失不大。

需要求使得\(\tilde g(x,y)\)化為密度的比例常數。

記N(\(\mu, \sigma^2)\)的分布密度為

\(f(x;\mu, \sigma^2)\)，對\(\tilde g(x,y)\)積分得

\[\begin{aligned}

\int_{-\infty}^\infty \int_{-\infty}^\infty \tilde g(x,y) =&

\sqrt{2\pi/90} \int_{-\infty}^\infty f(x;-0.4, 90^{-1}) dx

\cdot \sqrt{2\pi/120} \int_{-\infty}^\infty f(y;0.5,120^{-1}) dy \\

& + 0.5 \sqrt{2\pi/180} \int_{-\infty}^\infty f(x;0.5,180^{-1}) dx

\cdot \sqrt{2\pi/20} \int_{-\infty}^\infty f(y;-0.1,20^{-1}) dy \\

=& \sqrt{2\pi/90}\sqrt{2\pi/120} + 0.5 \sqrt{2\pi/180}\sqrt{2\pi/20} \\

\approx& 0.1128199

\end{aligned}\]

于是令

\[\begin{aligned}

g(x,y) =& \tilde g(x,y) / 0.1128199 \\

=& 0.5358984 f(x;-0.4,90^{-1}) f(y;0.5,120^{-1}) \\

& + 0.4641016 f(x;0.5,180^{-1}) f(y;-0.1,20^{-1}), \\

& -\infty < x < \infty, -\infty < y < \infty,

\end{aligned}\]

用復合抽樣法對\(g(x,y)\)抽樣，然后用重要抽樣法得到估計值\(\hat I_3\)。

注意用重要抽樣法計算時需要\(X_i\)獨立同分布，

但是其次序并不重要，

當\(N\)很大時，

可以固定取混合分布中兩個分布的樣本個數嚴格等于混合比例，

而不需要從兩個分布中隨機抽取。

一次模擬的估計值為：

估計的平均相對誤差為:

\(N=10000\)時一次模擬得到的\(\hat I_3 = 0.1253\),

從這一次模擬估計的平均相對誤差為0.002, 估計精度有兩位有效數字以上，

估計結果可報告為0.13。

重要抽樣法與平均值法的方差相比：

兩種方法的方差差距有200倍以上，

說明重要抽樣法節省了200倍的計算量。

※※※※※

例12.3標準化的重要抽樣法在貝葉斯統計推斷中有重要作用。

例如，設獨立的觀測樣本\(Y_j\)服從如下的貝塔—二項分布：

\[\begin{align}

& f(y_j | K, \eta)

= P(Y_j = y_j) \\

=& \binom{n_j}{y_j}

\frac{B(K\eta + y_j,\ K(1-\eta) + n_j - y_j)}{B(K\eta,\, K(1-\eta))},

\ y_j=0,1,\dots, n_j,

\tag{12.8}

\end{align}\]

其中\(B(\cdot, \cdot)\)是貝塔函數，

\(n_j\)為已知的正整數，\(K>0\)和\(0

貝塔—二項分布用于描述比二項分布更為分散的隨機變量分布。

按照貝葉斯統計的做法， 假設參數\((K, \eta)\)也是隨機變量，

具有所謂的“先驗分布”， 假設\((K, \eta)\)有如下的“無信息”先驗分布密度:

\[\begin{align}

\pi(K, \eta) \propto \frac{1}{(1+K)^2} \frac{1}{\eta(1-\eta)},

\tag{12.9}

\end{align}\]

則\((K, \eta)\)有如下的“后驗密度”:

\[\begin{align}

\tilde p(K, \eta | \boldsymbol Y) \propto& \pi(K, \eta) \prod_{j=1}^n f(y_j | K, \eta) \nonumber \\

\propto& \frac{1}{(1+K)^2} \frac{1}{\eta(1-\eta)}

\prod_{j=1}^n

\frac{B(K\eta + y_j,\ K(1-\eta) + n_j - y_j)}{B(K\eta,\, K(1-\eta))}.

\tag{12.10}

\end{align}\]

要求

\(E (\log K | \boldsymbol Y) = \int_0^\infty \log K \, \tilde p(K, \eta | \boldsymbol Y) \,dK\)的值。

如果可以從后驗密度\(\tilde p(K, \eta | \boldsymbol Y)\)直接抽樣，

可以用平均值法估計\(E (\log K | \boldsymbol Y)\)，

但從(12.10)來看很難直接抽樣。

為此，使用標準化的重要抽樣法。 為了解除\((K, \eta)\)的取值限制，

作變換\(\alpha = \log K\), \(\beta = \log \frac{\eta}{1-\eta}\),

則\(\alpha, \beta \in (-\infty, \infty)\),

而(12.10)對應的\((\alpha, \beta)\)的后驗密度為:

\[\begin{align}

p(\alpha, \beta | \boldsymbol Y)

\propto

\frac{e^\alpha}{(1 + e^\alpha)^2}

\prod_{j=1}^n

\frac{B(\frac{e^{\alpha}}{1 + e^{-\beta}} + y_j,

\ \frac{e^{\alpha}}{1 + e^{\beta}} + n_j - y_j)}

{B(\frac{e^{\alpha}}{1 + e^{-\beta}},

\ \frac{e^{\alpha}}{1 + e^{\beta}})}.

\tag{12.11}

\end{align}\]

取值無限制的隨機變量試抽樣密度經常使用自由度較小的t分布，

比如t(4)分布，設t(4)分布密度函數為\(g(\cdot)\),

用獨立的t(4)分布生成\((\alpha, \beta)\)

的試抽樣樣本\((\alpha_i, \beta_i), i=1,2,\dots,N\)，

可以估計\(E(\log K | \boldsymbol Y)\)為

\[\begin{aligned}

\hat \alpha = \frac{\sum_{i=1}^N \alpha_i

\frac{p(\alpha_i, \beta_i | \boldsymbol Y)}{g(\alpha_i) g(\beta_i)}}

{\sum_{i=1}^N

\frac{p(\alpha_i, \beta_i | \boldsymbol Y)}{g(\alpha_i) g(\beta_i)}}.

\end{aligned}

\tag{3.39}

\]

其中的\(p(\alpha_i, \beta_i | \boldsymbol Y)\)只要用(12.11)的右側計算，

因為分子和分母的歸一化常數可以消掉。

下面用R語言產生一組模擬的觀測數據\(\boldsymbol Y\)，然后對這組數據估計\(E(\log K | \boldsymbol Y)\)。

先定義貝塔-二項分布的概率質量函數:

其中\(n\)是正整數，\(y\)是非負整數且\(0 \leq y \leq n\)，

K和eta是上述分布參數\(K\)和\(\eta\)。

給定\(n\)個獨立觀測，

設這組觀測的\(n_1, \dots, n_n\)保存在向量narr中，

\(y_1, \dots, y_n\)保存在向量yarr中，

于是以\((\alpha, \beta)\)為自變量的后驗密度函數(差一個未知的常數倍)，

將\(\alpha\)和\(\beta\)記為自變量a, b，R函數為：

注意其中的narr和yarr是還沒有賦值的全局變量。

在R的函數定義中，

允許使用尚未定義的變量，

只要調用該函數時變量已經在該函數定義的環境中賦值就可以。

下面用模擬方法生成觀測樣本，然后認為觀測樣本已知。

將這些模擬生成的數據固定下來。

下面用標準化的重要抽樣法估計\(E(\log K | \boldsymbol Y)=E(\alpha | \boldsymbol Y)\)。

取\(\alpha\)和\(\beta\)的試抽樣分布為t(4)。

\(\log(K)\)真值為：

用10000個抽樣的重要抽樣法給出的\(\log(K)\)的后驗估計為：

重要性權重最好取值比較均勻，

否則重要性權重很小的抽樣點基本不起作用，

效率較低。

上述模擬的歸一化重要性權重的分布情況：

可見絕大多數抽樣點都貢獻很小。

這也是求解比較復雜的問題時重要抽樣法應用的普遍現象。

舍選控制是解決這樣問題的一種技術，

下面嘗試采用這一技術。先看權重的分為點：

取舍選控制閾值\(c\)為90%分位數。

對\(N\)個權重，超過閾值的都保留；

小于閾值的，

以概率\(W_i/c\)保留，\(c\)是閾值，

其它的丟棄，這樣減少了樣本量。

原來的舍選控制法是預先選好閾值\(c\)，

對每個\(i=1,2,\dots, N\)的每個，

從\(g(\boldsymbol x)\)中抽取后\(X_i\)后，

都進行舍選控制；如果第\(i\)個被丟棄，就重新從\(g(\boldsymbol)\)中抽取，直到被接受。

在R的向量化做法中，可以多抽取一些，

然后丟棄的就不再重復抽取。

舍選控制：

原來的10000個抽樣僅保留了1000多個。

權重修改為：

在新的權重下保留的抽樣是關于后驗密度適當加權的。

新的估計為:

※※※※※

總結

以上是生活随笔為你收集整理的舍选法抽样matlab,12 重要抽样法 | 统计计算的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。