统计分析基于R—随机模拟（附代码）

模拟是指把现实或抽象系统的某些特征或部分状态用模拟模型来替代。并根据具体问题使其变成一个概率模型求解。随机模拟方法就是对产生符合模型的大量随机数进行分析而求解问题。该方法又称为蒙特卡洛方法。最经典的蒙特卡洛方法就是通过计算所有落在正方形的点中落在圆形区域内的个数来估计 $\pi$ 的值。随机模拟中的随机性来自于模型本身的随机变量或者非随机问题转换为概率模型后特征量的估计问题，从而引入所谓的随机模拟误差。一般地，假设随机变量 $X$ 的期望为 $\theta$ ，方差为 $\sigma^2$ 。由中心极限定理知样本平均值 $\bar{X_N}=\frac{1}{N}\sum_{i=1}^NX_i$ 近似服从 $N(\theta,\frac{\sigma^2}{N})$ ，因此随机模拟误差幅度为。渐进标准差可用样本方差代替计算。要使标准差小于 $\sigma_0$ ，可以先取小样本计算出 $S^2_{N_0}=\frac{1}{N_0}\sum^{N}_{i=0}(X_i-\bar{X_N})^2$ ，然后根据不等式 $\frac{S_{N_0}}{\sqrt{N}}<\sigma_0,N>\frac{S_{N_0}^2}{\sigma^2_0}$ 计算需要的 $N$ 的大小。同时也可以根据中心极限定理计算 $\theta$ 的近似 $95$ % 置信区间： $\bar{X_N}\pm \frac{2S_N}{\sqrt{N}}$ 。

随机模拟方法避免了繁复的理论分析，但是其结果的随机性需要考虑在最终的分析之中。因为其算法简单条件易改变等特点，随机模拟用在科学研究中，常常作为探索性试验使用。

二、随机模拟积分

随机模拟解决非随机问题的典型问题是计算定积分，随机模拟积分也称蒙特卡洛积分。

随机投点法

几何直观上来说，计算函数的积分就是计算函数曲线下和坐标轴围成的面积。如名所示，随即投点法就是模拟产生平面上的个随机点，用这个点种落入曲线下方区域的百分比来估计概率，从而得到定积分的近似值。由强大数律可知估计值的精度可以随着样本个数的增加而无限提高，但与此同时付出的代价是计算量的成倍增加。

平均值法

相比于随机投点法，期望值的估计法效率更高。具体过程如下：取 $\left \{ U_i,i=1,2,... \right \}$ 独立同分布于 $U(0,1)$ ，则 $h(U_i)$ 也是独立同分布随机变量，且 $E\left [ h(U) \right ]=\int _a^bh(u)\frac{1}{b-a}du=\frac{I}{b-a}$ 。取 $\tilde{I}=\frac{b-a}{N}\sum^N_ih(U_i)$ 为 $I$ 的强相合估计，由强大数律 $\bar{Y}=\frac{1}{N}\sum^N_ih(U_i)\rightarrow E(h(U))=\frac{I}{b-a}$ 。平均值法本质上是将定积分问题转换为求函数期望问题。对随机变量 $X$ 的取值范围也不必局限于有限空间。 $\hat{Y}=\frac{1}{N}\sum^N_{i=1}h(U_i)$ 是函数期望的无偏估计和强相合估计，若随机变量方差存在，则 $\hat{I}$ 有渐进正态分布 $N(Eh(X),\frac{1}{N}Var(h(X)))$ 。

一维积分用数值方法均匀分布点计算一般更加有效，比如估计 $I$ 的复合梯形求积公式

$I_0=\frac{b-a}{N}\left [ \frac{1}{2}h(a)+\sum^{N-1}_{i=1}h(x_i)+\frac{1}{2}h(b)] \right ]$

上述两种计算定积分的方法均可以推广到高维情形。称用式子 $\hat{I_1}=\hat{p}\cdot M\prod ^d_{j=1}(b_i-a_i)$ 计算方法为随机投点法，其中 $\hat{p}$ 为落入下方区域的概率。称用式子 $\hat{I_2}=\frac{1}{N}\sum^N_{i=1}h(U_i)\cdot\prod^d_{j=1}(b_j-a_j)$ 计算方法为平均值法。通过计算两种方法的方差之差不难得出，平均值法的精度更高，这是因为随即投点法多用的随机数会导致抽样随机误差的增加。

维数诅咒问题

当维数增加时计算量迅猛增加以至于无法计算的问题称为维数诅咒。Monte Carlo积分方法之所以在高维问题中有重要应用，就是因为其精度和维度关系不大。相比之下，网格法就显得效率很低。

重要抽样方法

结合数学知识不难发现，只有在 $\left | h(x) \right |$ 较大的地方密集投点，在 $\left | h(x) \right |$ 较小的地方稀疏投点才可以有效利用样本避免浪费。称 $g(x)$ 为试投密度或重要抽样密度若其满足以下特征：其形状与 $\left | h(x) \right |$ 相似，且当 $g(x)=0$ 时 $h(x)=0$ ，当 $\left \| x \right \|$ 趋于无穷时 $h(x)=o(g(x))$ 。令 $\eta _i=\frac{h(X_i)}{g(X_i)},i=1,2,...$ ，则 $E\eta_i=\int \frac{h(x)}{g(x)}g(x)dx=\int h(x)dx=I$ ，因此可用此样本均值来估计 $I$ ，即 $\hat{I_3}=\frac{1}{N}\sum^N_{i=1}\frac{h(x_i)}{g(x_i)}$ 。这种方法称为重要抽样方法。事实上，计算 $\hat{I_3}$ 的方差后不难得出，当 $g(X)$ 和 $\left | h(x) \right |$ 形状接近时 $Var(\hat{I_3})$ 很小。

重要抽样方法适用于随机变量的随机数难以生成或者其分布集中于 $h(y)$ 接近零的位置。特别地，称 $W_i=\frac{f(X_i)}{g(X_i)}$ 为重要性权重。

进一步，当 $f(\cdot)$ 本身未知时，重新定义重要权重为 $\frac{\hat{f(x)}}{g(x)}$ ，其中 $\hat{f(x)}=cf(x)$ ，常数 $c$ 未知，可以将公式改成 $\hat{I_4}=\frac{\sum^N_{i=1}W_ih(X_i)}{\sum^N_{i=1}W_i}$ 。这称为标准化重要方法。注意，此时 $\hat{I_4}$ 是 $Eh(Y)$ 的强相合估计，但不是无偏估计。该方法在贝叶斯统计推断中有重要的作用。

在重要抽样方法的实际应用中，好的抽样分布很难获得，所以权重经常相差甚大，使得抽样样本主要集中在几个权重最大的样本上，为此可以舍弃权重太小的样本点，重新抽样替换以上样本点，这种方法被称为舍选控制。其具体过程如下：确定权重的阈值 $c$ ，然后对逐个样本点进行以下操作：若其权重大于阈值，则接受该样本点；否则以 $\frac{W_i}{c}$ 的概率接受。一轮操作后所有样本点的权重调整为 $W_i^*=p_cmax(W_i,c)$ ，其中 $p_c$ 为归一化常数。阈值的选取可以根据权重选择，如取为 $W_i$ 的某个分位数。

分层抽样

重要抽样方法选取了和形状相似的密度 $g(x)$ 作为试投密度，此举提高了精度，但是 $g(x)$ 有时难以找到。分层抽样法的思想就是把定义域上的积分分解为若干个子集上的积分，且 $h(x)$ 在每个子集上的振幅不大，分别计算积分后在求和。在每一子集上的积分也可用重要抽样法计算。一般地，设 $I=\sum^N_{i=1}\int_{C_i}h(x)dx$ ，在 $C_i$ 投随机点 $X_{ij}\sim U(C_i),j=1,...,n_i$ ，由此得 $I$ 的分层估计为 $\hat{I_5}=\sum^N_{i=1}\frac{V(C_i)}{n_i}\sum^{N_i}_{j=1}h(X_{ij})$ 。每一子集内振幅不大意味着 $\sigma_i^2$ 较小，此时可以通过随机点的个数选择使得估计量的方差在所有取法中取到最小。由此可以看出，分层抽样法的本质是把 $X$ 的值近似的抽样分入一层，从而减小估计方差。

三、方差缩减方法

随机模拟方法往往具有精度低、计算量大的缺点，减小随机模拟误差方差就可以有效地节省随机模拟时间。接下来以随机变量 $X$ 的期望 $\theta=EX$ 的估计为例探讨方差缩减方法。

控制变量法

由期望和方差的性质可知，若能找到零均值随机变量 $Y$ 且满足 $\rho_{xy}\neq 0$ ，就有 $Var(X+bY)=Var(X)+b^2Var(Y)+2bCov(XY)$ ，因此可以适当选取 $b$ 的值从而减小 $\theta$ 的估计方差。且 $X$ 和 $Y$ 的相关性越强，其改善幅度越大，这种方法被称之为控制变量法。实际中， $EY$ 的值可能是未知，可以先模拟一个小样本对上数值进行估计。另一种可能是 $EX=EY$ 且均未知，这时令 $Z=\alpha X+(1-\alpha)Y$ 同样可以对期望进行合理估计。

对立变量法

模拟中尝试用均匀分布随机数经函数 $g(U)$ 变换产生新的随机数。由数学知识可知，当该变换是单调的，则随机变量 $Y=g(1-U)$ 与负相关，且前后者同分布。因此可以通过求 $Y$ 的期望来估计。最常见的逆变换是 $F(X)$ ，其中设为连续分布函数， $\boldsymbol{U}\sim U(0,1)$ ， $X=F^{-1}(U)$ ， $Y=F^{-1}(1-U)$ ， $Z=\frac{X+Y}{2}$ ，则 $X$ 和 $Y$ 同分布 $F(X)$ 且 $Cov(X,Y)<0$ ， $Var(Z)<\frac{1}{2}Var(X)$ 。这种提高精度的方法叫做对立变量法，相比之下这种方法更简便易行。

条件期望法

上述两种方法均没有考虑到额外的相关信息对最终结果的影响。由重期望公式 $E\left \{ E{(Y|X)} \right \}=EY$ 和Rao-Blackwell不等式 $Var\left \{ E(Y|X) \right \}\leq Var(Y)$ 可知用 $E(Y|X)$ 的样本平均值来估计 $EY$ 的精度更高，这种方法称为条件期望法。一般来说 $Y$ 的抽样分为两步：首先对 $X$ 抽样，第二步利用 $Y|X$ 的条件分布对 $Y$ 抽样，用条件期望法则避免了第二步抽样引起的随机误差。

随机数复用

统计方法的性能指标包含偏差、方差、覆盖率等等，这些指标可以对不同的统计方法进行定量比较。为了减少计算量、提高结果精度，可以重复生成 $N$ 组随机变量，对每个样本同时用不同统计方法计算结果。

四、随机服务系统

随机服务系统模型应该包含以下三个要素（以食堂排饭为例）：

输入过程：比如已知同学进入食堂的规律

排队规则：比如同学选择哪个窗口的规则，如排队长短等

服务机构：有多少个窗口几个阿姨，做好一份要多久等等

按时间顺序记录发生的事件的模拟数据服务系统也叫做离散事件模拟，该算法主要可分为三类：活动模拟、事件模拟、过程模拟。顾名思义，这三类算法就是将连续时间离散化的不同思路。以事件模拟算法为例：当前时刻、随时间而变化的计数变量和系统状态这些变量仅在有事件发生时才更新，这样记录下来的就是时间点而不是一段时间。

五、MCMC

重要抽样方法虽然有较良好的高位推广，但是对于维数过高且分布复杂的情况则束手无策，MCMC（马氏链蒙特卡洛）方法模拟一个马氏链，使马氏链的平稳分布为目标分布，由此产生大量的近似服从目标分布但并非相互独立的样本。简单来说，马氏链的特点是在随机过程中某时刻的状态仅和上一时刻有关，设 $p_{ij}$ 为从上一时刻状态转移到状态的转移概率， $P=(P_{ij})_{m*m}$ 为转移概率矩阵。若有 $P(X_{t+k}=j|X_t=i)=p^{(k)}_{ij}$ 不依赖 $t$ 成立，则称为 $k$ 步转移概率。不可约马氏链中状态相互连通，且总能经过若干步相互转移。对只有有限个状态的非周期不可约马氏链有 $\lim_{n\rightarrow \infty }P(X_n=j)=\pi_j$ ，其中 $\pi_j$ 为常数， $\left \{ \pi_j \right \}$ 称为极限分布。

若可能的状态有可列个时，对状态 $i$ ，如果从该状态出发总能返回原态且所用时间期望是有限的，则称该状态是正态返的。所有状态均为非周期正常返状态的马氏链成为遍历马氏链，它存在唯一的极限分布。如果存在 $\left \{ \pi_j,j=1,2,...,m \right \},\pi_j\geq ,\sum^m_{j=1}\pi_j=1$ ，且满足 $\pi_ip_{ij}=\pi_jp_{ji},i\neq j$ ，则这样的马氏链为细致平衡的且 $\left \{ \pi_j \right \}$ 是该链的平稳分布。

遍历的不可约马氏链有平稳分布 $\left \{ \pi_j \right \}$ ，令该分布为目标分布，从高维或取值空间结构复杂的随机向量 $X$ 取值构造马氏链。从任意初值欻模拟产生的序列 $X_t$ 在 $t$ 足够大时其分布近似服从。抛弃开始一段老化期的部分抽样值，把剩余部分作为 $X$ 的非独立抽样。对于 $\theta = Eh(Y)$ ，可用估计量 $\hat{\theta}=\frac{1}{n-k}\sum^n_{k-1}h(X_t)$ 来估计。以上过程就是MCMC方法，注意到该估计是强相合估计但估计效率比独立抽样低。常见的MCMC算法有MH方法和Gibbs等，算法改进的方向是如何快速收敛到平稳分布且不会长时间在取值空间的局部区域内滞留。

Metropolis-Hastings抽样

具体过程如下：从 $\chi$ 中任意取初值，经过 $t$ 步后记当前状态为 $X^{(t)}$ ，下一步试转移概率分布为 $T(X|Y)$ ， $\mathit{U}\sim U(0,1)$ 按以下规则转移： $X^{(t+1)}=\left\{\begin{matrix} Y,U \leq r(X^{(t)},Y)\\ X^{(t)},else \end{matrix}\right.,r(x,y)=min\left \{ 1,\frac{\pi(y)T(x|y)}{\pi(x)T(y|x)} \right \}$

经验证，转移概率 $A(x|y)=P(X^{(t+1)}=y|X^{t}=x)$ 满足 MH 的细致平稳条件 $\pi(x)A(y|x)=\pi(y)A(x|y)$ ，大多数情况下 MH 构造的马氏链以 $\left \{ \pi(x) \right \}$ 为极限分布，注意目标分布可能只是确定到差一个常数倍。

如果，则称该算法为 $Metropolis$ 抽样法；若 $T(y|x)$ 不依赖于 $x$ 取值，则称相应算法为 $Metropolis$ 独立抽样法；对于连续型目标分布，设 $\pi(x)$ 为目标分布的密度， $T(y|x)$ 为给定 $x$ 条件下的试抽样密度，算法和离散型目标分布一致；由于 $T(y|x)$ 较难确定，随机游动 $Metropolis$ 抽样就是从作随机游动的试转移方法。令 $y=x^{(t)}+\varepsilon_{t},\varepsilon _t\sim g(x;\sigma )$ ， $\varepsilon _t$ 独立同分。注意随机游动的步幅 $\sigma$ 取值会影响结果，取值过大会导致拒绝率大，过小会使收敛到平衡态的速度过慢。