模拟

用实验的方法来确定概率值的方式称之为模拟。

而我们最常见的模拟即用计算机来模拟对应的概率分布，或者一个实验过程等。通常来说计算机都有一内置程序，成为随机数发生器，它产生一个伪随机数序列。就实用目的来说，这个伪随机数序列，与来自$(0,1)$均匀分布的样本没有区别。也就是相当于我们有一个均匀分布的随机变量。这个随机变量我们通常用$U$来表示。

模拟连续型随机变量的一般方法

本节主要介绍两种利用随机数模拟连续型随机变量的一般方法。

逆变换法

这个定理实际上就对应第5章问题26，29。证明方式一致，从分布函数出发。

$$ F_Y(a)=P\{Y\le a\}=P\{F^{-1}(U)\le a\} $$

注意$F$为分布函数一定是单调递增函数，所以逆函数也是单调递增函数。所以$F^{-1}(U)\le a$成立当且仅当$U\le F(a)$。于是上是变成：

$$ F_Y(a)=P\{U\le F(a)\}=F(a) $$

所以通常可以令$X=F^{-1}(U)$来模拟具有连续分布函数$F$得随机变量$X$。

舍取法

舍取法，其实类似于重复采样筛选出符合条件的值。这里更数学的表示这种方式的正确性。

假设有一种方式可以模拟密度函数为$g(x)$的随机变量，我们先模拟出一个密度函数为$g$的随机变量$Y$。然后以正比于$f(Y)/g(Y)$的概率采用$Y$的值，这样就能够模拟一个密度为$f(x)$的随机变量。

具体来说，如果存在一个常数$c$有

$$ \frac{f(y)}{g(y)}\le c \quad \text{for all }y $$

则可以通过如下方式来模拟具有密度$f$的随机变量

• 模拟具有密度$g$的$Y$，同时产生一随机数$U$。
• 若$U\le f(Y)/[cg(Y)]$，则$X=Y$，否则回到第一步。

证明

设$X$为舍取法模拟的随机变量，其由上述算法经过$N$次循环来模拟。

$$ P\{X\le x\}=P\{Y_N\le x\}=P\{Y\le x|U\le \frac{f(Y)}{cg(Y)}\} $$

根据算法步骤，最终$X\lt x$的情况，需在$U\le f(Y)/[cg(Y)]$条件下，使得最后生成的$Y$也小于$x$。即最后一部分条件概率。根据条件概率有：

$$ P\{Y\le x|U\le \frac{f(Y)}{cg(Y)}\}=P\{Y\le x,U\le \frac{f(Y)}{cg(Y)}\}/P\{U\le \frac{f(Y)}{cg(Y)}\} $$

现在考虑$Y,U$两个随机变量，其为独立随机变量。$Y$符合概率密度为$g$的分布，$U$为$(0,1)$上的均匀分布。则$Y,U$的联合概率密度为两个分布的乘积，即有：

$$ f(y,u)=g(y)\quad 0\lt u\lt 1 $$

令$K=P\{U\le f(Y)/cg(Y)\}$则原概率表达式有

$$ \begin{aligned} \frac{1}{K}P\{Y\le x,U\le \frac{f(Y)}{cg(Y)}\}&=\frac{1}{K}\int_{-\infty}^x\int_0^{f(y)/cg(y)}g(y)dudy\\ &=\frac{1}{cK}\int_{-\infty}^xf(y)dy \end{aligned} $$

注意其中$\int_{-\infty}^xf(y)dy$实际上就是$X$的概率分布表达式。现在令$x\to \infty$。可得

$$ 1=\frac{1}{cK}\int_{-\infty}^\infty f(y)dy=\frac{1}{cK} $$

这说明$cK$值必然为1。进而说明此时的$X$即具有概率密度为$f$的随机变量。

其实这个算法过程就相当于模拟一个量，随机一个指定位置，然后用$f(Y)/[cg(Y)]$比例来筛选是否合格。合格则接受这个值，否则拒绝。因为存在失败可能性，所以存在多次模拟可能性。

根据算法可见，每次循环独立地接受$Y$的概率为$P\{U\le f(Y)/[cg(Y)]\}=K=1/c$。可知循环次数为以$c$为均值的几何分布。

几个模拟例子

这里需要记录几个常用的模拟方式。

正态分布

在随机变量的联合分布章节中，其实有提过模拟正态分布的方式。若$X,Y$是相互独立的标准正态随机变量，则他们的极坐标$R=\sqrt{X^2+Y^2}$和$\Theta=\arctan(Y/X)$各自独立。$R^2$服从均值为2的指数随机变量。$\Theta$服从$(0,2\pi)$上的均匀分布。

所以设$U_1,U_2$为随机数则

$$ X=R\cos \Theta=(-2\ln U_1)^\frac{1}{2}\cos(2\pi U_2)\\ Y=R\sin \Theta=(-2\ln U_1)^\frac{1}{2}\sin(2\pi U_2)\\ $$

为两个独立的标准正态分布。这个方法也成为Box_Muller方法。

由于计算正弦值和余弦值，算法效率会收到影响。有一种方法可以避免潜在的费时的困难。这个方式就是通过舍取法模拟圆内的均匀分布，给出坐标$(V_1,V_2)$。因为是圆内的均匀分布，可以证明，其极坐标表示$\bar{R},\bar{\Theta}$相互独立，$\bar{\Theta}$服从$(0,2\pi)$上的均匀分布，$\bar{R}^2$服从$(0,1)$上的均匀分布。同时有

$$ \sin\Theta=\frac{V_2}{\sqrt{V_1^2+V_2^2}}=\frac{V_2}{\bar{R}} \quad \cos\Theta=\frac{V_1}{\sqrt{V_1^2+V_2^2}}=\frac{V_1}{\bar{R}} $$

进而可以化简原来模拟函数中的三角函数计算部分，即有：

$$ X=R\cos \Theta=(-2\ln U_1)^\frac{1}{2}V_2/\bar{R}\\ Y=R\sin \Theta=(-2\ln U_1)^\frac{1}{2}V_1/\bar{R} $$

还可以注意到$\bar{R}^2$服从$(0,1)$上的均匀分布这正好是均匀变量$U$的分布所以可以直接用$\bar{R}^2$来替换$U_1$即有：

$$ X=R\cos \Theta=(-2\ln \bar{R}^2)^\frac{1}{2}V_2/\bar{R}\\ Y=R\sin \Theta=(-2\ln \bar{R}^2)^\frac{1}{2}V_1/\bar{R} $$

这个方式是避免了三角函数的计算。但是因为舍取法的问题，会有重复计算的问题。

圆内均匀分布

游戏常见随机在圆内生成物体。这要求生成一个坐标$(x,y)$其服从半径为$1$内的均匀分布。简单来说其要生成两个联合随机变量$(X,Y)$满足如下的概率密度函数：

$$ f_{X,Y}(x,y)= \begin{cases} \frac{1}{\pi} \quad &\text{if }x^2+y^2\le 1\\ 0\quad &\text{other} \end{cases} $$

模拟圆分布的一个问题是，模拟的随机变量要限制在圆内，满足圆约束。或者说$X,Y$不是相互独立的。

一种方式也是通过舍取法来获得。即生成两个$(-1,1)$上的均匀随机变量$U_1,U_2$。其相当于在$(-1,1)$的正方形内均匀分布。如果$(x,y)$满足在圆内，则使用。否则舍弃。 这个方式是可以生成圆内均匀分布的随机变量。也是前面正态分布中用来优化三角函数选取的方式。可以知道其存在循环求取过程，循环次数等于$4/\pi$。

还有一种方式是尝试通过联合分布函数将$X,Y$拆解成两个可以独立的随机变量。再各自模拟对应的随机变量来求取。最直接的就是极坐标表示。

令$R=(X^2+Y^2)^{1/2}$和$\Theta=\arctan(Y/X)$。我们现在来看$R,\Theta$的分布。

根据随机变量函数的联合分布有

$$ f_{R,\Theta}(r,\theta)=f_{X,Y}(x,y)\left|J(x,y)\right|^{-1} $$

其中

$$ \left|J(x,y)\right|^{-1}=\begin{vmatrix} \frac{x}{\sqrt{x^2+y^2}}&\frac{y}{\sqrt{x^2+y^2}}\\ -\frac{y}{x^2+y^2}&\frac{x}{x^2+y^2}\\ \end{vmatrix}^{-1}=\sqrt{x^2+y^2} $$

可得

$$ f_{R,\Theta}(r,\theta)= \begin{cases} \frac{1}{2\pi}2r\quad &\text{if }0\le r\le 1\space 0\le\theta\le 2\pi\\ 0\quad &\text{other} \end{cases} $$

所以可以看到$R,\Theta$独立，$R$服从密度函数为$r/2$的分布函数。而$\Theta$为$(0,2\pi)$上的均匀分布。所以现在通过逆变法用随机数模拟$R,\Theta$再转换成$X,Y$即可。

由$f_R(r)=2r$可得$F_R(a)=a^2$其中$0\le a \le 1$。所以要用$U_1,U_2$模拟$R,\Theta$得分布可令

$$ R=F^{-1}(U_1)=\sqrt{U_1} \\ \Theta=2\pi U_2 $$

最后转化成$X,Y$即可：

$$ X=R\cos\Theta=\sqrt{U_1}\cos2\pi U_2\\ y=R\sin\Theta=\sqrt{U_1}\sin2\pi U_2\\ $$