重要性抽樣中的方差縮減技術—

榣山遺韻：重要性抽樣中降低方差的幾種方法

這一篇相比上一篇非常繞，我儘量寫的簡明有邏輯一點。上一篇中已經介紹了一種重要的方差縮減技術，重要性抽樣，並介紹了幾種簡單的提議重要性測度的方法。但是，我們還是希望有一個系統的提議方法。很多時候，由於對統計量的估計需要用到大量的樣本，所以，我們更關注樣本數n很大時估計量的漸進行為

一些定義

為什麼需要重要性抽樣

速率分析與大偏差近似

方差速率函式與重要性測度的有效性

例子——馬爾可夫鏈

一、一些定義

考慮下面的抽樣系統

$\Omega \stackrel{Z_{p(q), n} \sim P(Q)}{\longrightarrow} S^n \stackrel{f_n}{\rightarrow} \mathbb{R}^d$

，即n個樣本構成的隨機向量

把機率空間映到一個完備可分的度量空間，再由Borel函式

將其映到實空間。舉個例子，對於均值的估計，

$\begin{split} f_n : & S^n \rightarrow \mathbb{R}\\ &Z = (Z_1, Z_2, \cdots, Z_n) \mapsto \sum_{i=1}^n Z_i \end{split}$

，估計量

$\hat{m} = \frac{f_n(Z)}{n}$

。當然，如果是對連續隨機過程，則

應該改成

即從時間指標集映到狀態空間的全體函式，即樣本軌道的全體。

我們希望瞭解估計量的分佈情況，尤其是某些稀有事件（rare event）的機率，即

$\rho_n = P(\frac{f_n(Z_{p, n})}{n} \in E)$

。再舉個例子，比如非平衡統計物理中感興趣的是流的漲落，而流是樣本軌道的函式，常由隨機積分或者均值來定義，比如

$J_T = \frac{1}{T} \int_0^T d X_t$

，而想知道流的漲落則需要知道

$P(J_T \in E)$

。這一機率常常蘊涵著非平衡系統的響應等等性質。

二、為什麼需要重要性抽樣

Prop.

在置信度y下，要使估計量的相對誤差被控制在x以下至少需要的樣本數是

$k = (\frac{t_y \sigma_k}{x m_k})^2$

，其中，

$m, \sigma$

是被估計的量的均值、標準差，

是正態分佈的雙側分位數。

Pf.

根據中心極限定理

$P(|\frac{\hat{\rho}^k - m}{\sigma/\sqrt{k}}| \leq t_y) = P(|\frac{\hat{\rho}^k - m}{m}|\frac{m}{\sigma/\sqrt{k}} \leq t_y) = P(|\frac{\hat{\rho}^k - m}{m}| \leq x)= y$

如果用直接蒙特卡洛進行抽樣，

$\hat{\rho}_n^k = \frac{1}{k} \sum_{i=1}^k \mathbb{I}_{\frac{f_n(Z_{p, n}^{(i)})}{n} \in E}$

，由於稀有事件很難被擊中，所以導致直接蒙卡估計量的方差很大，從上面的命題可以看出，當估計量的方差很大時，所需要的樣本數也是極其大的。這時候，就需要引入重要性測度Q使得在P下的稀有事件在Q下不再稀有，從而達到降低方差，提高抽樣效率的目的。由於不管是在連續還是離散的情況下，我們所關心的n往往都比較大，所以，只需要考慮這一機率的漸進行為即可。下面介紹一種處理思路，即大偏差近似。有一個不太嚴格但是比較好說明大偏差近似是什麼的推理。

三、速率分析與大偏差近似

根據中心極限定理，對於n個零均值、方差為1的iid rv

，

$P_n(\bar{X} \geq l) = P(\frac{\sum_{i=1}^n X_i / n}{\sqrt{1/n}} \geq \sqrt{n} l )\stackrel{n\rightarrow \infty}{\longrightarrow} \frac{1}{\sqrt{2\pi}}\int_{\sqrt{n} l}^{\infty}e^{-\frac{1}{2}x^2}\mathrm{d}x \sim e^{-n \frac{1}{2} l^2+o(n)} \sim e^{-nI(l)+o(n)}$

，即對於偏離均值較大的估計量的機率漸進的正比於

$e^{-nI(l)}$

，其中

表示指數衰減的速率，我們把它叫速率函式。速率函式越大，當樣本數增大時，其衰減的也就越快。當速率函式為0時，在n趨於無窮的過程中，其機率也不會衰減到0，實際上也就表示估計量依機率收斂到速率函式的零點。更嚴謹的定義可以Google，這裡就不去花太多功夫說了。而在樣本數n逐漸變大的過程中，估計量在一個區域內的機率的衰減速率是由速率函式的值最小的點所控制的，因此，對於任何區域E，定義速率函式

$I(E) = \lim_{n \rightarrow \infty} - \frac{1}{n} \ln P(\frac{f_n(Z)}{n} \in E) = \inf_{x \in E} I(x)$

，在上面的例子中，由於在邊界l處衰減最慢，所以速率函式

$I([l, +\infty)) = I(l)$

。對於速率函式的計算，給出一個重要的定理。

Thm. (Gartner-Ellis)

對於任意

$\theta \in \mathbb{R}^d$

，定義累積生成函式

$\phi(\theta) = \lim_{n \rightarrow \infty} \frac{1}{n} \ln \mathbb{E}[e^{\langle \theta, f_n(\cdot) \rangle}]$

，其中尖括號表示內積，則速率函式為累積生成函式的Legendre-Fenchel變換，

$I(x) = \sup_{\theta \in \mathbb{R}^d} \{\langle \theta, x\rangle - \phi(\theta)\}$

因此，區域的速率函式為

$I(E) = \inf_{x \in E} \sup_{\theta \in \mathbb{R}^d} \{\langle \theta, x\rangle - \phi(\theta)\}$

對於重要性抽樣的估計量，

$\hat{\rho}_n^k = \frac{1}{k} \sum_{i=1}^k \frac{d P}{d Q}(Z_{q, n}^{(i)}) \mathbb{I}_{\frac{f_n(Z_{q, n}^{(i)})}{n} \in E}$

，我們

Def.

定義函式

$c_n(\theta) = \frac{1}{n} \ln \int_{S^n} \frac{d P_n}{d Q_n}(z) \exp{(\langle \theta, f_n(z)\rangle)} P_n(dz)$

，如果定義機率

$\mu_n(B) = \exp{(-n c_n(0))}\int_B \frac{d P_n}{d Q_n}(z) P_n(dz)$

，則在

$\mu_n$

下的累積生成函式

$\phi(\theta) = \lim_{n \rightarrow \infty} \phi_n(\theta) = \lim_{n \rightarrow \infty} \frac{1}{n} \ln \mathbb{E}_{\mu_n}[e^{\langle \theta, f_n(\cdot)\rangle}] = \lim_{n \rightarrow \infty} c_n(\theta) - c_n(0) = c(\theta) - c(0)$

Rmk.

可以證明

是一個凸函式。

根據Gartner-Ellis定理，

$I(x) = \lim_{n \rightarrow +\infty} - \frac{1}{n} \ln \mu_n(\frac{f_n(Z_n)}{n} \in dx)$

可以由

$\phi$

的Legendre-Fenchel變換給出，即

$\begin{split} I(E) &= \inf_{x\in E} I(x) = \inf_{x\in E}\sup_{\theta \in \mathbb{R}^d} \{ \langle \theta, x\rangle - \phi(\theta)\}\\ &= \inf_{x\in E}\sup_{\theta \in \mathbb{R}^d} \{ \langle \theta, x\rangle - c(\theta)\} + c(0)\\ &= R(E) + c(0) \end{split}$

這裡，記c的Legendre-Fenchel變換為R，實際上和I值相差一個規範化常數，所以，R是一個測度的速率函式，但不是機率測度。那麼，R到底是什麼測度的速率函式，和重要性抽樣又有什麼樣的關係？

四、方差速率函式與重要性測度的有效性

上一篇已經寫了分析重要性抽樣方差的恆等式

$k Var(\rho_Q) = F_Q - \rho^2$

，其中，被估計的量是rv

$\eta$

，

$\rho = \mathbb{E}[\eta], ~ F_Q = \mathbb{E}[\eta^2 \frac{d P}{d Q}]$

為方差函式，也是評價重要性測度好壞的標準，方差函式越小，重要性測度Q所帶來的方差縮減效果也就越大。

在上面的情況中，

$\eta = \mathbb{I}_{\frac{f_n(\cdot)}{n} \in E}$

，示性函式的平方不改變，因此，方差函式

$F_Q(E) = \mathbb{E}[\mathbb{I}_{\frac{f_n(Z_q)}{n} \in E} \frac{d P}{d Q}(Z_q)]$

Prop.

方差速率函式

為方差函式

的速率函式，即

$\lim_{n \rightarrow \infty}\frac{1}{n} \ln F_n(E) = - R(E)$

Pf.

$\begin{split} F_n(E) &= \int_{S^n} \frac{d P_n}{d Q_n}(z) \mathbb{I}_{\frac{f_n(z)}{n} \in E} P_n(dz)\\ &= \exp{(nc_n(0))} \mu_n(\{z \in S^n : \frac{f_n(z)}{n} \in E\})\\ \lim_{n \rightarrow +\infty} \frac{1}{n} \ln F_n(E) &= c(0) - I(E)\\ &= - R(E) \end{split}$

Def.

當方差速率函式

等於被估計機率本身速率函式

的兩倍，即

，稱重要性測度是有效的

Rmk.

還是從恆等式來看

$k Var(\rho_Q) = F_Q - \rho^2 = F_n(E) - P(\frac{f_n(Z)}{n} \in E)^2 \sim e^{- n R(E)} - e^{- 2 n I(E)} = e^{- 2 n I(E)} (e^{- n (R(E) - 2 I(E))} - 1)$

當有效性條件成立時，重要性抽樣的估計量的方差漸進為0

五、例子——馬爾可夫鏈

1。馬爾可夫鏈MC

$\{X_n : n\in \mathbb{N}_+\}$

的狀態空間

$S = \{0, 1, \cdots, k-1\}$

，轉移機率

$p_{ij}$

，函式

$f : S \rightarrow \mathbb{R}^d$

，

$f_n(X_1, X_2, \cdots, X_n) = \sum_{i=1}^n f(X_i)$

。記狀態空間上

$S \rightarrow \mathbb{R}^d$

的函式全體為

$\mathcal{G}$

2。定義線性運算元

$T_\theta : \mathcal{G} \rightarrow \mathcal{G}$

使得

$T_\theta (g)(x) = \mathbb{E}_P^x[e^{\langle\theta, f(\cdot)\rangle}g(\cdot)] = \sum_{y=0}^{k-1} e^{\langle\theta, f(y)\rangle}g(y)p_{xy}$

。由於

$T_\theta$

是嚴格正運算元，由Perron-Frobenius定理，存在無簡併的主特徵值

$\lambda(\theta) = \sup_{g:\|g\|=1} \|T_\theta g\|$

，和對應的單位特徵函式

$r_\theta(x)$

。矩母函式

$\begin{split} M_{f_n}(\theta) &= \mathbb{E}[\exp{(\langle \theta, f_n(\cdot) \rangle)}]\\ &= \sum_{x_1} p(x_1)e^{\langle \theta, f(x_1)\rangle} \prod_{i=2}^n \sum_{x_i} p_{x_{i-1}x_i} e^{\langle \theta, f(x_i) \rangle}\\ &= \sum_{x_1} T_\theta^{(n-1)}(1)(x_1)p(x_1)e^{\langle \theta, f(x_1)\rangle} \end{split}$

對函式g用

$T_\theta$

的特徵函式展開，

$T_\theta^{(n-1)} (g)(x) = c(g)\lambda(\theta)^{n-1} r_\theta(x) + o(n)$

，則可以證明主特徵值的對數等於累積生成函式

$\phi(\theta) = \lim_{n\rightarrow +\infty} \phi_n(\theta) = \lim_{n\rightarrow +\infty} \frac{1}{n} \ln [\lambda(\theta)^{n-1}(\sum_{x_1} c(1) r_\theta(x_1) p(x_1) e^{\langle\theta, f(x_1)\rangle} + o(n))] = \ln \lambda(\theta)$

3。同樣的，對於線性運算元

$K_{q, \theta}(g)(x) = \mathbb{E}_Q^x[\frac{p_{x\cdot}^2}{q_{x\cdot}^2} e^{\langle\theta, f(\cdot)\rangle} g(\cdot)] = \sum_{y=0}^{k-1} \frac{p_{xy}^2}{q_{xy}} e^{\langle\theta, f(y)\rangle} g(y)$

，其主特徵值的對數等於c

$\begin{split} c(\theta) &= \lim_{n\rightarrow +\infty} \frac{1}{n} \ln \mathbb{E}_Q[\frac{p(x_1)^2}{q(x_1)^2}e^{\langle\theta, f(x_1)\rangle}\prod_{i=2}^n\frac{p_{x_{i-1}x_i}^2}{q_{x_{i-1}x_i}^2}e^{\langle\theta, f(x_i)\rangle}]\\ &= \lim_{n\rightarrow +\infty} \frac{1}{n} \ln \mathbb{E}_{q(x_1)}[K_{q, \theta}^{(n-1)}(1)(x_1) \frac{p(x_1)^2}{q(x_1)^2}e^{\langle\theta, f(x_1)\rangle}]\\ &= \ln \gamma_q(\theta) \end{split}$

Prop.

$\nu = \arg\min_{x \in E} I(x)$

為最小速率點，其Legendre對偶

$\theta_\nu$

使得

$\nabla \phi(\theta_\nu) = \nu$

，則

$q_{o,xy} = p_{xy}\frac{e^{\langle \theta_{\nu}, f(y) \rangle}}{\lambda(\theta_{\nu})} \frac{r_{\theta_\nu}(y)}{r_{\theta_\nu}(x)}$

是有效的重要性轉移機率

Pf.

$\begin{split} K_{q_o, \theta}(g)(x) &= r_{\theta_\nu}(x) \lambda(\theta_{\nu}) \sum_{y=0}^{k-1} p_{xy} e^{\langle\theta-\theta_\nu, f(y)\rangle} \frac{g(y)}{r_{\theta_\nu}(y)}\\ &= r_{\theta_\nu}(x) \lambda(\theta_{\nu}) T_{\theta-\theta_\nu}(\frac{g}{r_{\theta_\nu}})(x)\\ &= \lambda(\theta_\nu)\lambda(\theta-\theta_\nu)r_{\theta_\nu}(x)r_{\theta-\theta_\nu}(x) \end{split}$