13 分层抽样法

用平均值法计算\(\int_C h(\boldsymbol x) \,d\boldsymbol x\)，若\(h(\boldsymbol x)\)在\(C\)内取值变化范围大则估计方差较大。重要抽样法选取了与\(f(x)\)形状相似但是容易抽样的密度\(g(\boldsymbol x)\)作为试投密度，大大提高了精度，但是这样的\(g(\boldsymbol x)\)有时难以找到。

如果把\(C\)上的积分分解为若干个子集上的积分，使得\(h(\boldsymbol x)\)在每个子集上变化不大，分别计算各个子集上的积分再求和，可以提高估计精度。这种方法叫做分层抽样法。这也是抽样调查中的重要技术。

例13.1 对函数 \[\begin{aligned} h(x) = \begin{cases} 1 + \frac{x}{10}, & 0 \leq x \leq 0.5 \\ -1 + \frac{x}{10}, & 0.5 < x \leq 1 \end{cases} \end{aligned}\] 求定积分 \[\begin{aligned} I = \int_0^1 h(x)\,dx, \end{aligned}\] 可以得\(I\)的精确值为\(I = 0.05\)。我们用平均值法和分层抽样法来估计\(I\)并比较精度。

在\((0,1)\)区间随机抽取\(N\)点用平均值法得\(\hat I_2\)，其渐近方差为 \[\begin{aligned} \text{Var}(\hat I_2) = \frac{\text{Var}(h(U))}{N} = \frac{143}{150N} \approx \frac{0.9533}{N} . \end{aligned}\]

模拟的R程序:

h <- function(x) ifelse(x <= 0.5, 1 + 0.1*x, -1 + 0.1*x)
I.true <- 0.05
set.seed(1)
N <- 10000
eta <- h(runif(N))
I2 <- mean(eta)
sd2 <- sd(eta)

取\(N=10000\)，一次平均值法得到的估计和相对误差为

I2
## [1] 0.0594168
(I2 - I.true)/I.true
## [1] 0.1883359

相对误差为0.2，仅有一位有效数字。估计值可报告为0.06。估计的\(N\text{Var}(\hat I_2)\)为

sd2^2
## [1] 0.9504117

理论值为0.9533。

把\(I\)拆分为\([0,0.5]\)和\([0.5,1]\)上的积分，即 \[\begin{aligned} I = a + b = \int_0^{0.5} h(x)\,dx + \int_{0.5}^1 h(x)\,dx , \end{aligned}\] 对\(a\)和\(b\)分别用平均值法，得 \[\begin{aligned} \hat a =& \frac{0.5}{N/2} \sum_{i=1}^{N/2} h(0.5 U_i) = \frac{0.5}{N/2} \sum_{i=1}^{N/2} (1 + 0.05 U_i), \\ \hat b =& \frac{0.5}{N/2} \sum_{i=(N/2)+1}^{N} h(0.5 + 0.5 U_i) = \frac{0.5}{N/2} \sum_{i=(N/2)+1}^{N} (-1 + 0.05 + 0.05 U_i), \\ \hat I_5 =& \hat a + \hat b, \end{aligned}\] 则分层抽样法结果\(\hat I_5\)的渐近方差为 \[\begin{aligned} \text{Var}(\hat I_5) =& \text{Var}(\hat a + \hat b) = \text{Var}(\hat a) + \text{Var}(\hat b) \\ =& 0.25 \frac{\text{Var}(1 + 0.05U)}{N/2} + 0.25 \frac{\text{Var}(-0.95 + 0.05U)}{N/2} = \frac{1/4800}{N}, \end{aligned}\] 分层后的估计方差远小于不分层的结果, 可以节省样本量约4500倍。

分层抽样法的一次模拟计算的R程序如下：

N <- 10000
set.seed(1)
eta1 <- h(runif(N/2, 0, 0.5))
a <- 0.5*mean(eta1)
eta2 <- h(runif(N/2, 0.5, 1))
b <- 0.5*mean(eta2)
I5 <- a+b

一次模拟的估计为：

I5
## [1] 0.0500084

\(N \text{Var}(\hat I_5)\)的估计为

Nvar <- 0.5*var(eta1) + 0.5*var(eta2); Nvar
## [1] 0.0002117651

理论值是0.002083。

※※※※※

一般地，设积分\(I=\int_C h(\boldsymbol x)\,d\boldsymbol x\)可以分解为\(m\)个不交的子集\(C_j\)上的积分, 即 \[\begin{aligned} I = \int_C h(\boldsymbol x)\,d\boldsymbol x = \int_{C_1} h(\boldsymbol x)\,d\boldsymbol x + \int_{C_2} h(\boldsymbol x)\,d\boldsymbol x + \dots + \int_{C_m} h(\boldsymbol x)\,d\boldsymbol x \end{aligned}\] 在\(C_j\)投\(n_j\)个随机点\(X_{ji} \sim \text{U}(C_j)\), \(i=1,\ldots,n_j\), 则\(I\)的\(m\)个部分可以分别用平均值法估计，由此得\(I\)的分层估计为 \[\begin{aligned} \hat I_5 = \sum_{j=1}^m \frac{V(C_j)}{n_j} \sum_{i=1}^{n_j} h(X_{ji}) \end{aligned}\] 记\(\sigma_j^2 = \text{Var}(h(X_{j1}))\), 划分子集时应使每一子集内\(h(\cdot)\)变化不大，即\(\sigma_j^2\)较小。这时 \[\begin{aligned} \text{Var}(\hat I_5) = \sum_{j=1}^m \frac{V^2(C_j) \sigma_j^2}{n_j} \end{aligned}\] 若\(\sigma_j^2\)可估计，应取\(n_j\)使 \[\begin{align} n_j \propto V(C_j) \sigma_j, \tag{13.1} \end{align}\] 即 \[\begin{aligned} n_j = N \frac{V(C_j) \sigma_j}{\sum_{k=1}^m V(C_k) \sigma_k}, \ j=1,2,\dots,m \end{aligned}\] 这样取的样本量\((n_1, n_2, \dots, n_m)\) 在所有满足\(n_1 + n_2 + \dots + n_m=N\)的取法中使得渐近方差最小。

在分层抽样法中，划分了子集后，每一子集上的积分也可用重要抽样法计算。

分层抽样法也可以用在求随机变量函数期望的问题中。设\(X\)为随机变量，要求\(X\)的函数\(h(X)\)的数学期望\(\theta = E h(X)\)。假设存在离散型随机变量\(Y\)，\(p_j = P(Y=y_j), j=1,2,\dots,m\), 在\(Y=y_j\)条件下可以从\(X\)的条件分布抽样，则 \[\begin{align} E[h(X)] = E\{E[h(X)|Y]\} = \sum_{j=1}^m E[h(X) | Y=y_j] p_j, \end{align}\] 如果在\(Y = y_j\)条件下生成\(X\)的\(N_j = N p_j\)个抽样值，设为\(X_i^{(j)}, i=1,2,\dots, N_j\), 则可以用\(\frac{1}{N_j} \sum_{i=1}^{N_j} h(X_i^{(j)})\)估计 \(E[h(X) | Y=y_j]\)，估计\(\theta\)为 \[\begin{align} \hat\theta = \sum_{j=1}^m \frac{1}{N_j} \sum_{i=1}^{N_j} h(X_i^{(j)}) p_j = \frac{1}{N} \sum_{j=1}^m \sum_{i=1}^{N p_j} h(X_i^{(j)}), \tag{13.2} \end{align}\] 这是\(\theta\)的无偏和强相合估计，且估计方差 \[\begin{align} \text{Var}(\hat\theta) =& \frac{1}{N^2} \sum_{j=1}^m N p_j \text{Var}[h(X) | Y=y_j] \nonumber \\ =& \frac{1}{N} \sum_{j=1}^m \text{Var}[h(X) | Y=y_j] p_j \tag{13.3} \\ =& \frac{1}{N} E \{ \text{Var}[h(X) | Y] \} \leq \frac{1}{N} \text{Var}[h(X)], \end{align}\] 比直接用平均值法估计\(E h(X)\)的方差小。这里用到了条件方差的性质 \[\begin{align} \text{Var}(X) = E \left[ \text{Var}(X|Y) \right] + \text{Var} \left[ E(X|Y) \right] \geq E \left[ \text{Var}(X|Y) \right], \end{align}\] 如果\(Y\)与\(X\)独立则\(E(X|Y)=EX\), \(\text{Var} \left[ E(X|Y) \right] = 0\), 这时分层抽样法比平均值法没有改进。从(13.3)可以看出，如果第\(j\)层样本的函数\(h(X_i^{(j)}), i=1,2,\dots,Np_j\)的样本方差为\(S_j^2\), 则\(\text{Var}(\hat\theta)\)的一个无偏估计是 \[\begin{align*} \widehat{\text{Var}(\hat\theta)} = \frac{1}{N} \sum_{j=1}^m S_j^2 p_j. \end{align*}\]

公式(13.2)取第\(j\)层样本数\(N_j = N p_j\)，仅考虑了\(Y\)的取值分布，而未考虑\(X | Y=y_j\)的条件分布情况。类似于(13.1)，应该对\(\text{Var}[h(X) | Y=y_j]\)较大的层取更多的样本。使得估计方差最小的分层样本量分配满足\(N_j \propto p_j \sqrt{\text{Var}[h(X)|Y=y_j]}\), 即 \[\begin{align} N_j = N \frac{ p_j \sqrt{\text{Var}[h(X)|Y=y_j]}}{\sum_{k=1}^m p_k \sqrt{\text{Var}[h(X)|Y=y_k]}}. \tag{13.4} \end{align}\] 在\(\text{Var}[h(X)|Y=y_j]\)未知的时候，可以预先抽取一个小的样本估计\(\text{Var}[h(X)|Y=y_j]\), 然后按估计的最优\(N_j\)分配各层的样本量。采用(13.4)的分层样本量后， \[\begin{align*} \hat\theta =& \sum_{j=1}^m \frac{1}{N_j} \sum_{i=1}^{N_j} h(X_i^{(j)}) p_j, \\ \text{Var}(\hat\theta) =& \sum_{j=1}^m \frac{p_j^2 \text{Var}[h(X)|Y=y_j]}{N_j}, \end{align*}\] 于是\(\text{Var}(\hat\theta)\)的估计为 \[\begin{aligned} \widehat{\text{Var}(\hat\theta)} =& \sum_{j=1}^m \frac{p_j^2 S_j^2}{N_j}, \end{aligned} \] 其中\(S_j^2\)是第\(j\)层样本函数\(\{ h(X_i^{(j)}), i=1,2,\dots,N_j \}\)的样本方差。

分层抽样法的本质是把\(X\)的值相近的抽样分入一层，使得同层的\(X\)条件方差较小，从而减小估计方差。

例13.2 设\(U \sim \text{U}(0,1)\), 用分层抽样法估计\(\theta = E h(U) = \int_0^1 h(x) \,dx\)。

令\(Y = \text{ceiling}(mU)\), 即当且仅当\(\frac{j-1}{m} < U \leq \frac{j}{m}\)时\(Y=j\)， \(j=1,2,\dots,m\), 可以按照\(Y\)分层抽样估计\(\theta\): \[\begin{aligned} \theta =& E[h(U)] = \sum_{j=1}^m E[h(U) | Y=j] \, P(Y=j) \\ =& \frac{1}{m} \sum_{j=1}^m E[h(U) | Y=j], \end{aligned}\] 易见\(Y=j\)条件下\(U\)服从\((\frac{j-1}{m}, \frac{j}{m})\)上的均匀分布，设\(U_1, U_2, \dots, U_n\)是U(0,1)的独立抽样，则用分层抽样法取每层\(N_j=1\)估计\(\theta = E h(U)\)为 \[\begin{aligned} \hat\theta =& \frac{1}{m} \sum_{j=1}^m h\left( \frac{j-1 + U_j}{m} \right). \end{aligned}\]

※※※※※

习题

习题1

设\(\sigma_j, j=1,2,\dots, m\)为\(m\)个正实数, \[\begin{aligned} f(\alpha_1, \alpha_2, \dots, \alpha_m) =& \frac{\sigma_1^2}{\alpha_1} + \frac{\sigma_2^2}{\alpha_2} + \dots + \frac{\sigma_m^2}{\alpha_m}, \ (\alpha_1, \dots, \alpha_m) \in (0,1)^m, \end{aligned}\] 则在\(\alpha_1 + \alpha_2 + \dots + \alpha_m = 1\)条件下 \(f(\alpha_1, \alpha_2, \dots, \alpha_m)\)的最小值点为 \[\begin{aligned} \alpha_j = \frac{\sigma_j}{\sum_{k=1}^m \sigma_k}, \ j=1,2,\dots,m \end{aligned}\] 最小值为\((\sigma_1 + \dots + \sigma_m)^2\)。

习题2

考虑定积分 \[\begin{aligned} I = \int_{-1}^1 e^x dx = e - e^{-1}. \end{aligned}\]

用随机模拟方法计算定积分\(I\), 分别用随机投点法、平均值法、重要抽样法和分层抽样法计算。
设估计结果为\(\hat I\), 如果需要以95%置信度保证计算结果精度在小数点后三位小数，这四种方法分别需要计算多少次被积函数值？
用不同的随机数种子重复以上的估计\(B\)次，得到\(\hat I_j, j=1,2,\dots,B\), 由此估计\(\hat I\)的抽样分布方差，与2的结果进行验证。
称 \[\begin{aligned} \mbox{MAE}(\hat I) = E | \hat I - I| \end{aligned}\] 为\(\hat I\)的平均绝对误差。从3得到的\(\hat I_j, j=1,2,\dots,B\)中估计\(\mbox{MAE}(\hat I)\)。比较这四种积分方法的平均绝对误差大小。

习题3

设\(h(x) = \frac{e^{-x}}{1 + x^2}\), \(x \in (0, 1)\), 用重要抽样法计算积分\(I = \int_0^1 h(x) \,dx\)，分别采用如下的试抽样密度: \[\begin{aligned} f_1(x) =& 1, \ x \in (0, 1), \\ f_2(x) =& e^{-x}, \ x \in (0, \infty), \\ f_3(x) =& \frac{1}{\pi (1 + x^2)}, \ x \in (-\infty, \infty), \\ f_4(x) =& (1 - e^{-1})^{-1} e^{-x}, \ x \in (0, 1), \\ f_5(x) =& \frac{4}{\pi (1 + x^2)}, \ x \in (0, 1). \end{aligned}\]

作\(h(x)\)和各试抽样密度的图形，比较其形状。
取样本点个数\(N=10000\), 分别给出对应于不同试抽样密度的估计\(\hat I_k\), \(k=1,2,3,4,5\), 以及\(\text{Var}(\hat I_k)\)的估计。
分析\(\text{Var}(\hat I_k)\)的大小差别的原因。
把\((0,1)\)区间均分为10段，在每一段内取\(N=1000\)个样本点用平均值法计算积分值，把各段的估计求和得到\(I\)的估计\(\hat I_6\)，估计其方差。
用例13.2的分层抽样方法计算积分的估计\(\hat I_7\)，估计\(\text{Var}(\hat I_7)\)并与前面的结果进行比较。