舍选法抽样matlab,12 重要抽样法 | 统计计算
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 重要抽样法 | 统计计算的全部內容，希望文章能夠幫你解決所遇到的問題。
- 上一篇： php 模板引擎 优点,Smarty模板
- 下一篇： matlab人脸追踪,求大神帮助我这个菜