数值分析 $9 随机数与应用

§9 随机数与应用

C1随机数

1）随机数：独立分布的数集

• 伪随机数不独立，使用纯数学方法生成的是伪随机数

2）线性同余生成器(LCG)：$x_i=(a{x}_{i-1}+b)\%m$，a为乘子，b为偏移，m为模数，x₀称为种子

• 最小标准随机数生成器：$m=2^{31}-1,a=7^5=16807,b=0$

  形如$2^p-1$的数称为梅森素数

• randu生成器：$a=65539=2^{16}+3,m=2^{31}$

  $(x_{i+2}-6x_{i+1}+9x_i)\%m=0$,系数太小导致独立性不足

• MATLAB不使用LCG，使用的是延迟斐波那契生成器，重复周期超过$2^{1400}$

3）指数随机数：分布满足$p(x)=a{e}^{-ax},a>0$的随机数

• 均匀随机数$u$转指数随机数$x=\frac{-\ln (1-u)}{a}$
• 由$u\sim U[0,1]$生成满足$p(x)$分布的随即变量：$P(x)={\int }_0^{x}p(x)dx,P^{-1}(u)$为所求

4）标准正态随机数：分布满足$p(x)=\frac{1}{\sqrt{2\pi }}{e}^{-\frac{x^2}{2}}$的随机数

• Box-Muller生成法：$n=\sqrt{-2\ln u_1}\sin (2\pi u_2)或\sqrt{-2\ln u_1}\cos (2\pi u_2),u_1,u_2\sim U[0,1]$
• 避免使用三角函数的修正：$x_1,x_2\sim U[0,1],u_1=x_1^2+x_2^2$，且当$u_1<1$时重新选取。$u_2=\arctan (\frac{x_2}{x_1})$
• MATLAB：randn()

C2 蒙特卡罗模拟

1）蒙特卡洛模拟

• Ⅰ型蒙特卡罗：计算函数均值$\frac{1}{b-a}{\int }_a^{b}f(x)dx$

• Ⅱ型蒙特卡罗：计算满足条件的点集面积。使得特征函数取1的概率

• $\mathrm{E}\mathrm{r}\mathrm{r}\mathrm{o}\mathrm{r}\propto n^{\frac{1}{2}}$

2）拟随机数：数据间不独立，但自避(不聚集)

• 使用拟随机数的蒙特卡罗算法有更快的收敛速度，令$d$拟随机数的维度

  • Ⅰ型：$\mathrm{E}\mathrm{r}\mathrm{r}\mathrm{o}\mathrm{r}\propto (\ln n{)}^d{n}^{-1}$
  • Ⅱ型：$\mathrm{E}\mathrm{r}\mathrm{r}\mathrm{o}\mathrm{r}\propto n^{-\frac{1}{2}-\frac{1}{2d}}$

  以Ⅱ型为例，误差的主要来源是边界上点的取值。

  使用$n$个点能够将$d$维空间分为$n^{\frac{1}{d}}$个网格点，在边界上的点数正比于$n^{\frac{d-1}{d}}$

  因此估计方差为$n^{\frac{d-1}{d}}÷n=n^{-\frac{1}{d}}$，标准差$n^{-\frac{1}{2d}}$，取均值后$n^{-\frac{1}{2}-\frac{1}{2d}}$

• p进制低差异序列法：将数$1-n$按$p$(素数)进制写出，记$(b_k\dots b_0{)}_p$，则生成的拟随机数为$(0.b_0\dots b_k{)}_p$

C3 布朗运动

1）一维随机游动：

• $E(x)=0,D(x)=n$

• 逃逸：首次到达$[-b,a]$边缘。$E(n)=ab$。在$a$端逃逸的(后验)概率为$\frac{b}{a+b}$

2）连续布朗运动：

• 每个周期随机游动的步数变为$k$倍，步长缩减为$\frac{1}{\sqrt{k}}$，记为$W_t^k$。$E(W_t^k)=0,D(W_t^k)=t$

  令$k\to +\infty$，得到连续布朗运动$B_t$

• $\forall t_1,t_2,t_1<t_2,(B_{t_1}-B_{t_2})\sim N(0,t_2-t_1)$，且与一切$B_t$独立

C4 随机微分方程

1）随机微分方程（SDE)：解为随机过程的微分方程。如含有噪声项的微分方程

2）Ito积分：噪声项。${\int }_a^{b}f(s)d{B}_s=\underset{\Delta t\to 0}{\lim }f(t_{i-1})\Delta B_i,\Delta B_i=B_{t_i}-B_{t_{i-1}}$

• 布朗运动的噪声$d{B}_t$称白噪声

3）Ito公式：用于解随机微分方程。$y=f(t,x)$，则$dy=\frac{\partial f}{\partial t}dt+\frac{\partial f}{\partial x}dx+\frac{1}{2}\frac{{\partial }^{2}f}{\partial x^2}dxdx$

• 通过$d{B}_{t}d{B}_t=dt,dtdt=0,dtd{B}_t=d{B}_{t}dt=0$可计算出$dxdx$

4）数值积分法：

• 解法阶数：$E(|\hat{y}-y|)=O((\Delta t{)}^m)$

• Euler-Maruyama法:

  $\{\begin{array}{l}dy(t)=f(t,y)dt+g(t,y)d{B}_t\\ y(0)=y_0\end{array}$

  $w_0=y_0$
  for i = 0, 1, 2, $\cdots$
  $w_{i+1}=w_i+f(t_i,w_i)\Delta t_i+g(t_i,w_i)\Delta B_i$
  end

  其中，生成$\Delta B_i=z_i\sqrt{\Delta t_i},z_i\sim N(0,1)$是正态随机数

  阶数为$\frac{1}{2}$

• Milstein法：

  $w_0=y_0$
  for i = 0, 1, 2, $\cdots$
  $w_{i+1}=w_i+f(t_i,w_i)\Delta t_i+g(t_i,w_i)\Delta B_i+\frac{1}{2}\frac{\partial g}{\partial y}(\Delta B_i^2-\Delta t_i)$
  end

  $\frac{1}{2}\frac{\partial g}{\partial y}(\Delta B_i^2-\Delta t_i)$即随机泰勒序列

  阶数为1，收敛更快。但是需要用户提供偏导数

• 一阶随机龙格-库塔法：

  $w_0=y_0$
  for i = 0, 1, 2, $\cdots$
  $w_{i+1}=w_i+f(t_i,w_i)\Delta t_i+g(t_i,w_i)\Delta B_i+\frac{1}{2\sqrt{\Delta t_i}}[g(t_i,w_i+g(t_i,w_i)\sqrt{\Delta }{t}_i)-g(t_i,w_i)](\Delta B_i^2-\Delta t_i)$

  end