首页 > 其他分享> > 计算机模拟备考纲要

计算机模拟备考纲要

2021-01-26 16:36:16 作者：互联网

本文基于《随机模拟方法与应用》，适用于备考浙江大学计算机模拟（3学分，课程号82120010），参考了较多王何宇老师讲义的内容。基于不同教材的计算机模拟侧重点不同，考点也不同，请勿完全参考。由于内容是一个人整理的，缺少审阅，难免出现错误，如有发现，请在评论区指正。所有的代码都可以在导入所需的库后直接运行，但如果是画图代码，则需要自行添加plt.show()。

大纲如下：

随机数的生成
- 均匀分布随机数
- 离散随机数生成
- 逆变换法
- 拒绝接受法
- 正态分布抽样
MCMC
- Metropolis-Hastings抽样法
- 概率转移矩阵的选择
- Ising模型
- Metropolis算法的原理
MC积分
- 随机投点法
- 重要性抽样法
最优化算法
- 模拟退火法
- 遗传算法
实例模拟
- 服务系统
- 随机游走

使用的库：

import numpy as np
import numpy.random as npr
import matplotlib.pyplot as plt
from collections import Counter
from scipy import stats

方便起见，笔记中不含图，可以结合课件或讲义插图配合理解。分布由于是复习笔记，所以内容较为简略。

Part 1：简单随机数的生成

1-1 \(U(0,1)\)随机数

平方取中法：取位数\(s\)，种子设为\(2s\)位的整数\(X_0\)，基于以下算法生成随机数：

\[X_{i+1}=\left[\frac{X_i^2}{10^s}\right]\mod{10^{2s}}. \]

def rand_mid_square(X_0, s):
    X = np.zeros(N)
    X[0] = X_0
    for i in range(1, N):
        X[i] = X[i-1]**2 // 10**s % 10**(2*s)
    return X/10**(2*s)

生成的随机数效果极差，会陷入循环。

Fibonacci产生器：生成算法为

\[X_i=(X_{i-2}+X_{i+1})\mod{M}. \]

def rand_add_mod(X_0, X_1, M, N):
    """取自讲义"""
    X = np.zeros(N)
    X[0] = X_0
    X[1] = X_1
    for i in range(2, N):
        X[i] = (X[i-2]+X_[i-1]) % M
    return X/(M-1)

线性同余产生器：基于以下算法生成随机数：

\[X_i=(AX_{i-1}+C)\mod{M}. \]

需要初值\(X_0\)与模型参数\(A,C,M\)，指定产生的随机数数量。

def rand_mul_mod(X_0, M, A, C, N):
    """取自讲义"""
    X = np.zeros(N)
    X[0] = X_0
    for i in range(1, N):
        X[i] = (A*X[i-1]+C) % M
    return X/(M-1)

一般将\(M\)取得很大如\(2^{16},2^{31}\)，这样产生的随机数具有足够的间距。

这两种方法具有的问题是随机数是在低维空间上的随机，可能在映射到高维空间上时具有规律性。IBM360的经典线性同余发生器模型参数是\(M=2^{31},A=65539,C=0\)，这样产生的随机数映射到高维空间具有规律性。

现在使用的均匀分布随机数生成算法主要为梅森旋转。均匀分布是实现下列所有随机数生成的基础。

Python中，常用于生成\(U(0,1)\)的函数有

npr.rand()  # 生成一个U(0, 1)
npr.rand(N)  # 生成N个U(0, 1)

1-2 离散分布随机数

离散分布的生成可以参考其阶梯型分布函数，基于\(U(0,1)\)生成随机数，将\([0,1]\)区间按照离散分布的概率分段即可。

def rand_discrete(pmf, N):
    cdf = np.cumsum(pmf) / np.sum(pmf)
    X = np.zeros(N)
    for i in range(N):
        h = npr.rand()  # 由于已经讨论过01均匀分布随机数的生成，直接用内置函数替代01均匀分布随机数
        index = 0
        while h > cdf[index]:
            index += 1
        X[i] = index + 1
    return X

这只适用于有限离散分布，而不适用于可列离散分布。

1-3 逆变换法

逆变换法希望从一个已知分布\(p(x)\)中抽样，步骤如下：

积分\(p(x)\)得到\(F(x)\)，求出反函数\(F^{-1}(x)\)；
生成01均匀分布随机数\(u\sim U(0,1)\)；
得到随机数\(x=F^{-1}(u)\)。

举例：生成指数分布随机数。

\[p(x)=\lambda e^{-\lambda x},\\ F(x)=1-e^{-\lambda x},\\ F^{-1}(u)=-\frac{1}{\lambda}\ln(1-u). \]

容易验证\(F^{-1}(F(x))=x\)，实现逆变换抽样：

def inv_exp(l, N):
    X = np.zeros(N)
    for i in range(N):
        u = npr.rand()
        x = -1/l*np.log(1-u)  # 这里l是均值的倒数
        X[i] = x
    return X

逆变换法的优点是简单快捷，缺点有：

并不是所有的分布函数逆变换在计算上都是高效的；
分布函数必须严格递增，否则反函数不存在；
有的分布不存在解析分布函数。

如果生成的均匀分布含边界值0和1，则需要格外注意。

逆变换法有效性的证明：设生成的随机数是\(X_i\)，则\(X_i=F^{-1}(U_i)\)，\(U_i=F(X_i)\)且\(U_i\sim U(0, 1)\)。所以

\[\mathbb{P}(X_i\le x)=\mathbb{P}(F^{-1}(U_i)\le x)=\mathbb P(U_i\le F(x))=F(x). \]

1-4 拒绝接受法(ARM)

拒绝接受法基于一个已知且容易抽样的分布\(h(x)\)，从另一个复杂分布\(f(x)\)抽样。\(f(x)\)可以作如下分解：

\[f(x)=Mg(x)h(x)，\\ 0\le g(x)\le 1. \]

这自然要求\(f\)有显然的上界，并且总有\(Mh(x)\le f(x)\)。拒绝接受法的算法步骤如下：

生成\(Z\sim h(x)\)的随机数\(z\)，它应该被接受的比例是

\[\delta(z)=\frac{f(z)}{Mh(z)}. \]
生成\(U\sim U(0,1)\)的随机数\(u\)；
如果\(u<\delta(z)\)，则接受随机数\(z\)，否则舍弃这个随机数。

讲义中给到的例子是从标准半正态分布中抽样：

\[p(x)=\sqrt{\frac{2}{\pi}}e^{-\frac{x^2}{2}}I_{0\le x<\infty}. \]

以指数分布为简单分布抽样，进行如下的分解：

\[p(x)=\sqrt{\frac{2e}{\pi}}\cdot e^{\frac{2x-x^2}{2}} \cdot e^{-x},\\ M=\sqrt{\frac{2e}{\pi}},h(x)=e^{-x}. \]

实现ARM抽样：

def arm_half_norm(N):
    num = 0
    M = np.sqrt(2*np.e/np.pi)
    X = np.zeros(N)
    while num < N:
        Z = npr.exponential(1)
        U = npr.rand()
        accept = np.sqrt(2/np.pi)*np.exp(-Z**2/2) / (M*np.exp(-Z))
        if U < accept:
            X[num] = Z
            num += 1
    return X

AR抽样法有效性的证明：将\((U,Z)\)看成一个二维独立随机变量，且

\[U\sim N(0, 1),Z\sim h(x).\\ p_{U,Z}(u,v)=h(v)I_{0\le u\le 1}. \]

我们实际采取的随机变量是\(Z\)，但只有在\(U\le \delta(Z)\)的条件下接受，所以要证明

\[\mathbb{P}\left(Z\le z\bigg|U\le \frac{f(Z)}{Mh(Z)} \right)=F(z), \]

实际上

\[\begin{aligned} & \mathbb{P}\left(Z\le z\bigg|U\le \frac{f(Z)}{Mh(Z)} \right) \\ =& \frac{\mathbb{P}\left(Z\le z,U\le \frac{f(Z)}{Mh(Z)}\right)}{\mathbb{P}\left(U\le \frac{f(Z)}{Mh(Z)} \right)} \end{aligned} \]

分别求分子分母的概率，注意到\(U,V\)是独立的：

\[\begin{aligned} & \mathbb{P}\left(U\le \frac{f(Z)}{Mh(Z)} \right)\\ =& \int_{-\infty}^\infty\mathbb{P}\left(U\le \frac{f(z)}{Mh(z)}\bigg|Z=z \right)h(z){\rm d}z \\ =&\int_{-\infty}^\infty\frac{f(z)}{Mh(z)}h(z){\rm d}z\\ =&\frac{1}{M},\\ & \mathbb{P}\left(U\le \frac{f(Z)}{Mh(Z)}, Z\le z \right)\\ =& \int_{-\infty}^z\int_{0}^{\frac{f(v)}{Mh(v)}}h(v){\rm d}u{\rm d}v \\ =&\frac{F(z)}{M}, \end{aligned} \]

代入得到

\[\mathbb{P}\left(Z\le z\bigg|U\le \frac{f(Z)}{Mh(Z)} \right)=\frac{\frac{F(z)}{M}}{\frac{1}{M}}=F(z). \]

AR抽样法的优点是：不需要分布函数的反函数。缺点是：

需要舍弃许多随机数\(z\)，尤其是\(z\)发生在\(\delta(z)\)很小的地方时；
要能找到能完全覆盖\(f(z)\)的简单抽样函数\(Mh(z)\)，这一点不容易做到。

如果\(Mh(z)\)不能很好地拟合\(f(z)\)，AR法的效果就会比较差。

1-5 正态分布\(N(0,1)\)抽样

正态分布的简单抽取方式：中心极限定理，找到期望为0方差为1的随机变量\(U_i\sim U(-\frac{1}{2},\frac{1}{2})\)，有

\[\mathbb{E}(U_i)=0,\mathbb{D}(U_i)=\frac{1}{12}. \]

由中心极限定理，对独立同分布的\(U_i\)，有

\[\sqrt{\frac{12}{n}}\sum_{i=1}^nU_i\stackrel {d}\to N(0,1). \]

如果\(n\)很大，可以用这种方式近似生成标准正态随机数。

def norm_simple(N, n):
    X = np.zeros(N)
    C = np.sqrt(12/n)
    for i in range(N):
        U = npr.rand(n)
        X[i] = C*np.sum(U)
    return X

正态分布精准抽样：如果\(U\sim U(0,1)\)，\(V\sim E(1)\)且相互独立，则如下定义的\(X,Y\)都服从\(N(0,1)\)：

\[X=\sqrt{2V}\cos(2\pi U),Y=\sqrt{2V}\sin(2\pi U). \]

证明：

\[p(u,v)=e^{-v}I_{0\le u\le 1}I_{v\ge 0},\\ \frac{J(x,y)}{J(u,v)}=\begin{bmatrix} -2\pi\sqrt{2V}\sin(2\pi U) & \frac{1}{\sqrt{2V}}\cos(2\pi U) \\ 2\pi\sqrt{2V}\cos(2\pi U) & \frac{1}{\sqrt{2V}}\sin(2\pi U) \end{bmatrix},\\ \left|\frac{J(x,y)}{J(u,v)} \right|=2\pi, \left|\frac{J(u,v)}{J(x,y)} \right|=\frac{1}{2\pi}. \]

而

\[V=\frac{X^2+Y^2}{2},\\ U = \frac{1}{2\pi}\arctan(\frac{Y}{X}), \]

所以

\[\begin{aligned} f(x,y)=&p(u,v)\cdot\frac{1}{2\pi}\\ =&\frac{1}{2\pi}e^{-\frac{x^2+y^2}{2}}, \end{aligned} \]

这是\(N_2(0,0,1,1,0)\)的密度函数，所以\(X,Y\stackrel{i.i.d}\sim N(0,1)\)。

多元正态分布的抽取：由Cholesky分解，相关系数阵\(R\)可以分解为下三角矩阵\(L\)与其转置相乘，即

\[R=LL', \]

所以抽取\(n\)元正态分布向量\(X\sim N_p(0,I_p)\)，有

\[Y=LX+\mu\sim N_p(\mu,\Sigma). \]

Part 2：MCMC

2-1 Metropolis Hastings抽样步骤

Metropolis-Hastings抽样需要一个条件分布：\(g(\cdot|x_{t})\)，其步骤为

获得初始状态\(x_0\)；
基于\(x_t\)有条件分布\(g(\cdot|x_t)\)，从\(g(\cdot|x_t)\)中产生建议随机数\(y\)；
基于\(y\)产生一个接受概率\(h(x_t,y)\)，产生01随机数\(r\)；
如果\(r<h(x,y)\)，则\(x_{t+1}=y\)，否则\(x_{t+1}=x_t\)。

要解决的问题是条件分布\(g(\cdot |x_t)\)与接受概率\(h(x_t, y)\)的寻找。如果将\(g(\cdot |x_t)\)看成状态转移矩阵，则抽取的样本来自的总体，是\(g(\cdot|x_t)\)的平稳分布。

首先是，如何确定平稳分布\(f(y)\)，接下来我们假设将从\(f(y)\)中抽样。在M-S算法抽样中，平稳分布实际上依赖于\(h(x_t,y)\)。如果取

\[h(x_t,y)=\min\left\{1,\frac{f(y)g(x_t|y)}{f(x)g(y|x_t)} \right\}, \]

则产生的平稳分布就是\(f(y)\)。

接下来考虑\(g(\cdot |x_t)\)，为了使\(h(x_t,y)\)易于计算，常常选择\(g(\cdot |x_t)\)为对称分布，即

\[g(x|y)=g(y|x), \]

这样，就可以简化接受概率为\(h(x_t,y)\)为

\[h(x_t,y)=\min\left\{1,\frac{f(y)}{f(x_t)} \right\}. \]

此时的M-H抽样算法称为对称M-H抽样算法。

2-2 \(g(x|y)\)的选择与简单实例

基于对称MH抽样法，\(g(x|y)\)有两种典型的选择：极大邻域法与极小邻域法。

对于离散型分布，极小邻域法的状态转移矩阵为

\[\begin{pmatrix} \frac12 & \frac12 & 0 &\cdots & 0 & 0 \\ \frac12 & 0 & \frac12 & \cdots & 0 & 0 \\ 0 & \frac12 & 0 & \cdots & 0 & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & 0 & \cdots & 0 & \frac12 \\ 0 & 0 & 0 & \cdots & \frac12 & \frac12 \end{pmatrix}, \]

极大邻域法的状态转移矩阵为（总状态数为\(n\)）

\[\begin{pmatrix} \frac1n & \frac1n & \frac1n &\cdots & \frac1n & \frac1n \\ \frac1n & \frac1n & \frac1n & \cdots & \frac1n & \frac1n \\ \frac1n & \frac1n & \frac1n & \cdots & \frac1n & \frac1n \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ \frac1n & \frac1n & \frac1n & \cdots & \frac1n & \frac1n \\ \frac1n & \frac1n & \frac1n & \cdots & \frac1n & \frac1n \end{pmatrix}. \]

可以验证这两种状态转移矩阵\(g(x|y)\)都满足对称条件（其实就是对称矩阵）。在Part 1中提供的离散型轮盘赌抽样法中，只能对有限分布列进行抽样，而M-H抽样可以推广到可列离散型分布列。

由于M-H抽样中密度函数（分布列）的出现总是上下齐次的，所以不需要归一化，只要密度函数（分布列）能表征概率的相对大小即可。接下来对于掷双骰子实验，使用对称M-H抽样。

\(x\)	2	3	4	5	6	7	8	9	10	11	12
\(p\)	1	2	3	4	5	6	5	4	3	2	1

这是最小邻域法抽样的实现，注意状态转移的过程。

def mh_doubledice_1(N):
    f = np.array([0, 0, 1, 2, 3, 4, 5, 6, 5, 4, 3, 2, 1], dtype=np.int8)
    X = np.zeros(10000+N, dtype=np.int8)  # 需要一个趋于平稳的过程
    y = 0
    X[0] = npr.randint(2, 13, 1)
    for i in range(1, 10000+N):
        h = npr.rand()  # 状态转移
        if X[i-1] == 2:
            y = 2 if h < 0.5 else 3
        elif X[i-1] == 12:
            y = 12 if h < 0.5 else 11
        else:
            y = X[i-1] - 1 if h < 0.5 else X[i-1] + 1
        accept = np.min([1, f[y]/f[X[i-1]]])  # 接受概率的计算
        X[i] = y if npr.rand() < accept else X[i-1]
    return X[10000:]

Rand6 = mh_doubledice_1(10000)
plt.hist(Rand6, bins=11)

这是最大邻域法抽样的实现：

def mh_doubledice_2(N):
    f = np.array([0, 0, 1, 2, 3, 4, 5, 6, 5, 4, 3, 2, 1], dtype=np.int8)
    X = np.zeros(10000+N, dtype=np.int8)
    y = 0
    X[0] = npr.randint(2, 13, 1)
    for i in range(1, 10000+N):
        y = npr.randint(2, 13, 1)  # 均等产生新的随机数
        accept = np.min([1, f[y]/f[X[i-1]]])
        X[i] = y if npr.rand() < accept else X[i-1]
    return X[10000:]

Rand7 = mh_doubledice_2(10000)
plt.hist(Rand6, bins=11)

注意到以上两个抽样方式都从10000开始，这是因为Markov链到达平稳分布需要一定的时间，开头的有限项应当作为收敛前项，不被考虑作为随机数。

对于连续型随机变量的分布，如果随机变量的取值域是有限的，则最大邻域法就是在取值域上均匀取值。以Beta分布为例，Beta分布的概率密度函数为

\[f(x;\alpha,\beta)=\frac{1}{\beta(\alpha,\beta)}x^{\alpha-1}(1-x)^{\beta-1}I_{0<x<1}. \]

正如之前讨论的，不需要知道\(\beta(\alpha,\beta)\)的具体值。运用最大邻域法给出MH抽样算法，可以自行绘图查看效果。

def mh_beta(N, alpha, beta):
    f = lambda x:x**(alpha-1)*(1-x)**(beta-1)
    X = np.zeros(10000+N)
    X[0] = 0.5
    for i in range(1, 10000+N):
        y = npr.rand()
        accept = np.min([1, f(y)/f(X[i-1])])
        X[i] = y if npr.rand() < accept else X[i-1]
    return X[10000:]

Rand8 = mh_beta(1000000, 2, 3)
plt.hist(Rand8, bins=100, density=True, stacked=True)
x = np.linspace(0, 1, 10000)
plt.plot(x, stats.beta(2, 3).pdf(x))

如果连续型随机变量的取值域包含无穷值，则适合采用极小邻域法，一般是从当前取值的周围均值。这里采取柯西分布的绝对值进行抽样，即

\[f(x)=\frac{1}{2\pi(1+x^2)}I_{x\ge 0}, \]

现在\(x\)的取值域是无穷的但不是\(\mathbb{R}\)。最小邻域一般取当前值\(x_t\)周围的区间\([x_t-1,x_t+1]\)进行均匀抽样\(y\sim U(x_{t}-1,x_{t}+1)\)，但如果\(x_t-1<0\)，就令\(y=x_t\)即可。这种做法类似离散时候的边界值。

def mh_cauchy(N):
    f = lambda x:1/(1+x**2)
    X = np.zeros(10000+N)
    X[0] = 1
    for i in range(1, 10000+N):
        y = 2 * npr.rand() - 1 + X[i-1]
        if y < 0:
            X[i] = X[i-1]
            continue
        accept = np.min([1, f(y)/f(X[i-1])])
        X[i] = y if npr.rand() < accept else X[i-1]
    return X[10000:]

Rand9 = mh_cauchy(10000000)
plt.hist(Rand9, bins=1000, density=True, stacked=True)
x = np.linspace(0, 12, 10000)
plt.plot(x, stats.cauchy().pdf(x))

2-3 复杂实例：Ising模型

本质：从玻尔兹曼分布中抽样（铁磁体的状态），并计算随机变量函数的数学期望（磁矩）。如果\(N\)很大，则用轮盘赌法要计算很复杂的期望，但是用MCMC抽样就可以忽略正则化因子，将接受概率设置为一个比较好计算的值。

本文对教材中提到的二维\(N\times N\)周期边界二维网格来模拟一个铁磁体，其中每一个网格的取值为1或-1。显然，\(N\times N\)一共含有\(2^{(N^2)}\)种状态，每一种状态（样本点）记作\(S\)，\(\Omega\)是所讨论的样本空间。

对于每一个样本点\(S\)，要定义其能量函数\(E(S)\)，这依赖于邻域的选取。现在我们选择的邻域是十字邻域，即

\[\Lambda(i, j)=\{(i+1,j),(i-1,j),(i,j+1),(i,j-1) \}. \]

此时定义的能量函数为

\[E(S)=-\frac{1}{2}J\sum_{j=1}^N\sum_{i=1}^N \sum_{\mu\in\Lambda(i,j)}s(\mu)s(i,j), \]

也就是说，对于每一个格点，其能量值是它与周围四个点状态的对比：每有一个同号状态就加一，每有一个异号状态就减一；整个系统的能量与所有格点的能量加和成正比。

系统的能量总是越低越稳定的，但是铁磁体不会总处于最低能量，所以铁磁体的状态存在一个概率分布，即\(\Omega\)中的每一个样本点都有一定的存在概率。我们的目标其实是求磁体总磁矩量，对于磁体的每一个状态\(S\)，其磁矩为

\[M(S)=\frac{1}{N^2}\sum_{j=1}^N\sum_{i=1 }^Ns(i,j), \]

所以磁体总磁矩量是其期望，即每一个状态的磁矩乘以状态的产生概率。由于MCMC用于计算随机变量的期望，所以此问题可以用随机模拟来解决。

Ising模型的模拟关键在于从玻尔兹曼模型抽样，状态为\(S\)的概率服从玻尔兹曼分布，即

\[\mathbb{P}(S)=\frac{1}{Z}e^{-\frac{E(S)}{k_BT}},\\ \]

这里\(Z\)是正则化因子，\(k_B\)是一个物理常量。采用极小邻域法，即每一次只对一个格点进行扰动（相对而言，极大邻域法相当于重新初始化了一个解），那么从\(S\)扰动到\(S'\)，其接受概率为

\[h(S, S')=\min\left\{1,\frac{\mathbb{P}(S')}{\mathbb{P}(S)} \right\}=\min\left\{1,e^{\frac{E(S)-E(S')}{k_BT}} \right\}. \]

这样我们就能够从玻尔兹曼中抽样，并计算玻尔兹曼分布的数学期望了。

def initState(N):
    S = 2 * npr.randint(0, 2, [N, N])-1
    return S

def Energy(S, N):
    E = 0
    for i in range(N):
        for j in range(N):
            E += S[i,j]*S[i-1,j] + S[i,j]*S[(i+1)%N,j] + S[i,j]*S[i,j-1] + S[i,j]*S[i,(j+1)%N]
    E *= -0.5
    return E

def trans(S, N):
    S_trans = S.copy()
    i, j = npr.randint(0, N, 2)
    S_trans[i,j] = -S_trans[i,j]
    return S_trans

def Magnetic(S, N):
    return np.sum(S) / N**2

def Ising(N, trials, T):
    beta = 1 / T
    M = np.zeros(10000+trials)
    S = initState(N)
    M[0] = Magnetic(S, N)
    for i in range(1, 10000+trials):
        S2 = trans(S, N)
        accept = np.min([1, np.exp(beta*(Energy(S, N)-Energy(S2, N)))])
        S = S2 if npr.rand() < accept else S
        M[i] = Magnetic(S, N)
    return M[10000:]

Rand10 = Ising(3, 100000, 100)
plt.hist(Rand10, bins=19)

2-4 Metropolis的原理

基础：Markov链。

（齐次的）状态转移矩阵\(P\)。
平稳分布：如果某个分布\(\pi\)满足\(\pi P=\pi\)，则称之为平稳分布。
极限分布：对于初始分布\(p_0\)，极限分布为

\[p_\pi=\lim_{n\to \infty}p_0P^n. \]
对于齐次的非周期不可约Markov链，极限分布收敛于平稳分布。不可约指的是从任意一个状态出发都有一定概率转移到其他所有状态，非周期指的是不存在某最小正整数\(d>1\)使得\(P^d=P\)。

关于对称Metropolis有效性的证明，书上与讲义中给出的都是直观理解。

假设有一系列独立的个体位于状态空间的所有位置，每一个时间步下，这些个体都按照状态转移矩阵独立行走。经过\(t\)步行走后，处在\(x\)状态的个体数目记作\(m_t(x)\)，处在\(y\)状态的个体数目记作\(m_t(y)\)，则动态平衡要求

\[\forall x,y,\quad m_t(x)\mathbb{P}(x\to y)=m_t(y)\mathbb{P}(y\to x). \]

也就是说，在平衡态（记作\(e\)时刻），应该有

\[\forall x,y,\quad \frac{m_e(x)}{m_e(y)}=\frac{\mathbb{P}(y\to x)}{\mathbb{P}(x\to y)}. \]

因为\(\mathbb{P}(y\to x)=g(y|x)h(x,y)\)，即当前位置在\(x\)下一步在\(y\)，且被接受的概率。基于状态转移矩阵的对称性，有

\[\frac{m_e(x)}{m_e(y)}=\frac{\mathbb{P}(y\to x)}{\mathbb{P}(x\to y)}=\frac{h(y,x)}{h(x,y)}. \]

如果\(f(x)\ge f(y)\)，则\(h(y,x)=1\)，\(h(x,y)=f(y)/f(x)\)，所以

\[\frac{m_e(x)}{m_e(y)}=\frac{f(x)}{f(y)}; \]

如果\(f(x)<f(y)\)，则\(h(x,y)=1\)，\(h(y,x)=f(x)/f(y)\)，所以

\[\frac{m_e(x)}{m_e(y)}=\frac{f(x)}{f(y)}. \]

所以\(m_e(x)\propto f(x)\)，也就证明了Metropolis算法的有效性。

遍历性定理：齐次不可约非周期Markov链的一次实现平均值趋近于平稳分布的期望，这对于任意映射\(\xi:\Omega\to \mathbb{R}\)成立，即

\[\lim_{N\to \infty}\frac{1}{N}\sum_{i=1}^N\xi(x_i)=\mathbb{E}_\pi(\xi). \]

Part 3：MC积分法

3-1 随机投点法

计算高维定积分，这里\(D=[a_1,b_1]\times\cdots\times[a_d,b_d]\)，\(f(\boldsymbol{x})\in[0,M]\)：

\[I=\int_Df(\boldsymbol{x}){\rm d}{\boldsymbol{x}}=\int_{a_1}^{b_1}\cdots\int_{a_d}^{b_d}f(x_1,\cdots,x_d){\rm d}x_1\cdots{\rm d}x_d. \]

随机投点法产生若干个区域\(\Omega=D\times M\)上的均匀分布随机向量（点），如果一共产生了\(N\)个点，有\(V\)个落在\(f(\boldsymbol{x})\)下方，则积分值为

\[\hat I=M\prod_{i=1}^d(b_i-a_i)\cdot\frac{V}{N}\xlongequal{def}M|D|\cdot\frac{V}{N}. \]

使用随机投点法计算平面\(x+y+z=1\)与三个坐标平面所围成的封闭区域，也就是

\[I=\int_{0}^{1}\int_{0}^{1}(1-x-y){\rm d}x{\rm d}y. \]

这个积分的实际值是\(\frac 16\)，实现方式：

def random_point(N):
    Volumn = 0
    for i in range(N):
        x, y, z = npr.rand(3)
        if z < 1 - x - y:
            Volumn += 1
    return Volumn / N

随机投点法的精度十分有限，其原理可以通过\(\hat I\)的分布来计算。令

\[p=\frac{I}{M|D|} \]

则有\(V\sim B(N,p)\)，\(\hat I=\frac{V}{N}\cdot M|D|\)。所以

\[\mathbb{E}(\hat I)=\frac{\mathbb{E}(V)}{N}\cdot M|D|=pM|D|=I; \\ \mathbb{D}(\hat I)=\frac{(M|D|)^2\mathbb{D}(V)}{N^2}=\frac{M|D|-I}{N}. \]

这表明\(\hat I\)是无偏的，但标准差正比于\(\sqrt{1/N}\)，由切比雪夫不等式，

\[\mathbb{P}(|\hat I-I|\ge\varepsilon)\le\frac{\mathbb{D}(\hat I)}{\varepsilon^2}, \]

想要将精度增加一位，其方差就需要增加两位精度，也就是\(N\)要增加一百倍。进一步，如果要求\(\hat I-I\le \epsilon\)的概率不小于\(1-\alpha\)，则切比雪夫不等式过于粗糙，依据中心极限定理有\((V-Np)/\sqrt{Np(1-p)}\to N(0,1)\)，再结合精度求分位数构造区间估计。

3-2 重要性抽样

重要性抽样法的考虑方向是，不从均匀分布出发，而是采用一个更有效率的分布，这可能从直观上变得不那么具有几何意义，但从数学角度出发是完全可行的。

假设这个更有效率的分布是\(X\sim g_X(x)\)，要计算\(f(x)\)在\([a,b]\)上的积分，则\(X\)也应该定义在\([a,b]\)上，此时有

\[I=\int_a^b f(x){\rm d}x=\int_a^b \frac{f(x)}{g_X(x)}g_X(x){\rm d}x=\mathbb{E}\left(\frac{f(X)}{g_X(X)}\right). \]

也就将\(I\)表示成了

\[Y=\frac{f(X)}{g_X(X)} \]

的总体均值。由辛钦大数定律，\(Y\)的样本均值依概率收敛于总体均值，所以现在只要从总体\(X\)中取值，经过函数变化变成\(Y\)的观测值\(y\)，\(\hat I\)就是\(Y\)的样本均值。

如果对有限区间上的函数求积分，可以使用均匀分布密度作积分，这样就变成了求\(f(X)\)的样本均值，因此称为样本平均法，也就是

\[I=\int_{a}^{b} f(x){\rm d}x=(b-a)\mathbb{E}f(X),\quad X\sim U(a,b). \]

例如计算：

\[I=\int_0^{\infty}e^{-\frac{x^2}2}{\rm d}x, \]

这个积分的精确值是\(\sqrt{\pi/2}=1.2533141\)。取指数分布为抽样分布，则将其改写为

\[I=\int_{0}^{\infty} \frac{e^{-\frac{x^2}{2}}}{e^{-x}}e^{-x}{\rm d}x\xlongequal{def}\int_0^{\infty}h(x)e^{-x}{\rm d}x, \]

此时的\(h(x)=e^{-\frac{x(x-2)}{2}}\)，\(X\sim E(1)\)。

def sample_mean(N):
    V = npr.exponential(1, N)
    Y = np.exp(-V*(V-2)/2)
    return np.mean(Y)

事实上，重要性抽样法并没有很大改善精度问题，但解决了无界积分问题。积分估计量为

\[\hat I=\frac{1}{N}\sum_{i=1}^N\frac{f(X_i)}{g(X_i)},\quad X\sim g_X(x) \]

由数理统计的知识，样本均值的期望为总体期望，方差为总体方差的1/N，所以在精度上，依然是100倍样本点增加一位精度。但重要性抽样法的方差取决于重要性抽样总体方差，即

\[\mathbb{D}\left(\frac{f(X)}{g_X(X)}\right)\xlongequal{def}\mathbb{D}(h(X)). \]

如果能选择合适的\(g_X(X)\)使得\(h(X)\)的方差最小化，则可以一定程度上增加抽样效率。直观上，则需要\(h(x)\)的变化幅度尽可能小，这要求\(g(x)\)与\(f(x)\)的拟合程度比较好。

相对于随机投点法，重要性抽样在缩小误差方面作出的努力是，它缩小了总体的方差。如果总体标准差缩小了2两倍，则相当于节省了4倍的抽样努力。

误差控制还有一种方法：控制变量法。它与重要性抽样法的区别在于，它的原理基于积分的线性性质：

\[\int_{a}^{b} f(x){\rm d}x=\int_{a}^{b}[f(x)-g(x)]{\rm d}x+\int_{a}^{b} g(x){\rm d}x. \]

如果希望以上的式子能达到方差缩减的目的，则需要满足两条：第一是\(g(x)\)的积分是方便计算的，第二是\(f(x)-g(x)\)的震荡度要低于\(f(x)\)。如果这样，我们就可以用样本平均法求\(f(x)-g(x)\)的期望。当然，这要求\(f(x)\)是在有限区间上积分的，否则就要更换重要性函数。

Part 4：最优化算法

4-1 最优化问题

最优化问题指求某函数\(f(x)\)的最小值，不同的是\(x\)往往存在约束，这是实际情况导致的，也就是

\[\arg\min_{x\in\Omega}f(x). \]

如果\(f\)是一个关于\(x\)的显函数并且连续可导，一般通过求导的方式求最小值，但实际的最优化问题往往达不到这一点。另外，实际上存在着一些非线性程度很高的问题，这时的\(f\)存在大量的局部最优点。

全局最优：称\(x^*\in \Omega\)为全局最优解，如果\(\forall x\)都有\(f(x^*)\le f(x)\)。
局部最优：称\(\tilde x\)为局部最优解，如果\(\forall x\in U(\tilde x)\)都有\(f(\tilde x)\le f(x)\)。这里\(U(\tilde x)\)称为\(\tilde x\)的邻域，实际问题中往往根据实际情况设定。

典型非线性函数：Shekel函数。

\[f(x)=\sum_{i=1}^m\frac{1}{a_i(x-b_i)^2+c_i}. \]

它的各个峰值分别位于\(b_i\)，衰变率为\(a_i\)。

二维Shekel函数：

\[f(x_1, x_2)=\sum_{i=1}^m\frac{1}{a_{1i}(x_1-b_{1i})^2+a_{2i}(x_2-b_{2i})^2+c_i}. \]

为了找到问题的全局最优，基于蒙特卡洛方法，有模拟退火法与遗传算法两种常用的随机算法。

4-2 模拟退火(SA)

对模拟退火讨论可以从Ising模型中玻尔兹曼分布入手。在MCMC抽样时，提到对称MH抽样法的接受概率为

\[h(x_t, y)=\min\left\{1,\frac{f(x)}{f(y)}\right\}, \]

如果取\(f(x)\)为玻尔兹曼分布即

\[P(S)=\frac{1}{Z}\exp\left(-\frac{E(S)}{k_B T} \right), \]

那么新状态的接受概率为

\[h(S,S')=\min\left\{1,\exp\left(\frac{E(S')-E(S)}{k_B T}\right) \right\}, \]

这个概率将在接下来的讨论中发挥重要作用。

接下来进入模拟退火法的解释。假定组合优化问题具有能量函数\(f(x)\)，其状态概率服从玻尔兹曼分布，在固定温度下，每一次迭代都基于某种方式产生一个新解\(x'\)，被接受的概率就是上面所列的接受概率。为了让此状态转移过程成为Markov链，在固定温度下需要具有一定的迭代次数，同时也要满足对称MH算法的不可约遍历与对称性。

观察玻尔兹曼分布，每一个状态\(f(x)\)的产生概率与温度有关：如果\(T\)很高，则不同能量\(f(x)\)的产生概率差不多；如果\(T\)很低，则不同能量\(f(x)\)的产生概率就会有很大差异。因此，当\(T\)接近于0时，由概率分布就能导出\(f\)的最小值点。

接下来的问题是降温，降温方案称为冷却进度表，它包含以下几个方面：

初始温度\(T_0\)。初始函数\(T_0\)要足够大，否则容易落入局部最优。
温度\(T\)的衰减函数，即每一次降温过程是如何发生的，下面是两种常用衰减函数：

\[T_{k+1}=\alpha T_{k}; \\ T_{k+1}=\frac{c}{a+\ln T_k}. \]
衰减函数的衰减速率不宜过大。
每一温度下的迭代次数\(L_k\)，这相当于Markov链的长度；如果温度衰减足够快，有时可以省略这个参数。
终止温度\(T_f\)，要足够小。

综合以上要点，给出模拟退火法的算法步骤：

确定冷却进度表，生成初始解。
在每一个温度下进行\(L_K\)次迭代：
1. 基于原解的邻域产生一个新解；
2. 计算接受概率，决定是否接受新解。
  
  \[h = \min\{1, \exp[(E(S')-E(S))/T] \}. \]
达到\(L_k\)次迭代后降温。
返回最优解。

定理的有效性：如果在每一个温度下都能达到当前温度下的平稳分布，则算法能收敛到全局最优解。

4-3 遗传算法

遗传算法基于生物进化过程中的自适应演化，其核心在于生成一个初始种群，然后经过交叉、变异、选择，相当于进行了一轮生物进化。

染色体：用来代指一个解。本课程中常用的是二进制编码与全排列编码。
种群：由一定规模的个体（染色体）构成的一个群组。
适应度：染色体在当前背景下的优越程度，可以视为能量函数取反。
交叉：选择两条染色体，通过某种方式获得一条子代染色体，包含父代两条染色体的部分信息。
变异：选择一条染色体，按一定概率进行扰动。
选择：按照轮盘赌的方式对种群进行优胜劣汰。

编码方式：

二进制编码：用{0，1}字符产生的字符串表示问题的候选解。连续问题需要离散化，离散程度根据问题所需要的精度确定。
全排列编码：用一个序列的重排代表问题的候选解。

二进制编码：如函数求极值问题，将求值区间\([a, b]\)按照精度\(\epsilon\)划分成\((b-a)/\epsilon\)个区间，取一个最小\(k\)使得\(2^k>(b-a)/\epsilon\)，这就是二进制编码的位数。

单点交叉：找到一个交叉点和两个交叉序列，交换交叉点后面的部分即可。
变异：找到一个变异点，将01互换。

全排列编码：如TSP问题。

变异：随机选择两个点，交换这两个点上的数字。
部分映射交叉（PMX）：选择一段进行交叉，然后将非交叉部分的重复元素进行部分映射。映射由交叉部分决定。

如12345678和45687213的交叉，选择第四位以后的段落交叉，则初始交叉结果为

\[C_1 : 12347213,\\ C_2 : 45685678. \]
交叉部分建立了如下的映射：

\[5 \leftrightarrow (7)\leftrightarrow 1,\\ 2\leftrightarrow 6,\\ 3\leftrightarrow 8. \]
这里7被括号起来代表7是映射的中间元素。因此，第一条染色体中，重复的部分是123，按照映射换成568；第二条染色体中，重复的部分则对应为568，按照应设换成123。所以，得到的交叉结果是

\[C_1 :56847213, \\ C_2 :41235678. \]
顺序交叉（OC）：选择一段进行交叉，然后将原序列中交叉段以外的数字按顺序填补到剩余的地方。

如12345678和45687213的交叉，选择第四位以后的段落交叉，则初始交叉的结果为

\[C_1: ****7213,\\ C_2: ****5678. \]
此时\(C_1\)交叉前剩余数码的顺序是4568，\(C_2\)交叉前剩余数码的顺序是4123，所以得到的交叉结果是

\[C_1: 45687213,\\ C_2: 41235678. \]
基于位置的交叉（PBC）：随机选择一半（向下取整）的位置保持不变，然后将另一条染色体中不重复的元素按顺序填入，这时另一条染色体中没有被选中的元素保持不变，把上一条染色体中不重复的元素按顺序填入。

如12345678和45687213的交叉，假设第一步选中了2368四个位置，则这四个位置对应的元素（也是2368）保持不变，另一条染色体中这四个元素也保持不变。

\[C_1:*23**6*8,\\ C_2:**68*2*3. \]
将C1、C2中要移动的位置按顺序填入C2、C1即可。

\[C_1:42357618,\\ C_2:14685237. \]

综合以上要点，给出遗传算法的步骤：

确定参数，即种群大小、变异概率、最大迭代次数。
每一次迭代，对种群进行交叉、变异、选择，注意每次生成的种群是固定规模的。
直到达到最大迭代次数，返回剩下的最优解。

Part 5：实例模拟

5-1 随机服务系统

随机服务系统的三个组成部分：

输入：即顾客的到达过程。
排队过程：顾客是否会离开、排队机制是LCFS还是FCFS。
服务：即服务过程。

符号约定：排队系统用X/Y/Z/A/B/C表示，这里

X：顾客到达时间间隔的分布，主要是M，代表指数分布（满足无记忆性、无并发性，构成泊松过程）。
Y：服务时间分布，主要是M。
Z：并行服务台数量。
A：排队系统的最大容量，默认为\(\infty\)。
B：顾客源数量，默认为\(\infty\)。
C：排队规则，默认为FCFS。

最典型的排队模型是M/M/1的，即顾客到达服从指数分布、服务时间服从指数分布，且只有一个服务台。

衡量指标：用来评判一个服务系统的好坏。

平均顾客数\(L_s\)：指系统中的顾客数，包括在服务的与排队中的。
平均队列长\(L_q\)：指系统中排队中的顾客数。
平均逗留时间\(W_s\)：指一个顾客在系统中停留时间的平均值。
平均等待时间\(W_q\)：指一个顾客在系统中排队等待时间的平均值。
平均忙期\(T_b\)：指服务机构连续繁忙时间的平均值。

算例模拟：一般是给定一个时间，用一个二维数组guests存储相关的信息，各行作用如下：

顾客的到来时间，可以先产生到来时间再用cumsum累加。
对每个顾客的（理论）服务时间。
顾客在队列中的等待时间。
顾客的离开时间。
系统的附加信息（Option），如属于哪一个服务台服务、是否接受服务等等。

另外，对每一个顾客i，都要计算到来时系统中存在的人数n，同时要将已经运算完毕的顾客进行存储。各个指标之间的计算：

始终存在guests[3, i] = guests[0, i] + guests[1, i] + guests[2, i]。
系统中已经存在的人数：guests[0, i]比多少个guests[3, j]大，这里j in range(i)。
如果到来时系统中存在的人数小于服务台数，则guests[2, i] = 0。
如果服务台数量不是1而是k，要求等待时间，要将目前为止的离开时间进行降序排列，然后取出第k大的离开时间作为服务时间。

构造一种较复杂的算例：服务台数量为4，每一次服务时间服从\(U(3, 5)\)，客户的到达间隔服从\(E(1)\)。如果顾客看到队列长度\(q>2\)，即系统总人数\(n>6\)则会以\(p=0.3\)的概率直接离开。

def serve(T):
    # T为设定的服务时间，假设为720。
    guests = np.zeros([5, 2*T])  # 构造一个5行2T列的数组，一般认为到来的人数不超过2*lambda*T
    
    # %% 设定前两行的参数
    guests[0,] = npr.exponential(1, 2*T)
    guests[0,] = np.cumsum(guests[0,])
    guests[1,] = npr.uniform(3, 5, 2*T)
    total_ser = np.max(np.where(guests[0,] < T))
    
    # %% 对前4个用户0,1,2,3进行初始化，因为有4个服务台
    for j in range(4):
        guests[2, j] = 0  # 不需要等待
        guests[3, j] = guests[0, j] + guests[1, j]  # 计算离开时间
        guests[4, j] = j+1  # 代表为其服务的服务台
        
    ser_list = [0, 1, 2, 3]  # 这代表已经计算完毕的顾客
    
    # %% 对剩下的人进行入队模拟
    for j in range(4, total_ser+1):
        system_num = np.sum(guests[3, :j] > guests[0, j])  # 计算当前队列中的人数
        if system_num < 4:  # 系统中人数比服务台总数小，代表有可用服务台
            guests[2, j] = 0
            guests[3, j] = guests[0, j] + guests[1, j]
        else:
            if system_num > 6:  # 系统中人数大于6，触发离开判定
                h = npr.rand()
                if h < 0.3:  # 马上离开
                    guests[2, j] = 0
                    guests[3, j] = guests[0, j]
                    continue
              # 留在队列中
            temp_finish_list = np.sort(guests[3, ser_list])  # 将目前为止选择进入服务的人的离开时间升序
            guests[3, j] = temp_finish_list[-4] + guests[1, j]  # 倒数第4个完成服务的时间作为其服务时间
            guests[2, j] = guests[3, j] - guests[0, j] - guests[1, j]  # 等待时间
        
        for server in [1,2,3,4]:
            index = np.max(np.where(guests[4,] == server))
            if guests[3, index] <= guests[0, j] + guests[2, j] + 0.000001:  # 加入修正项，防止浮点数运算问题
                guests[4, j] = server
                break
        
        ser_list.append(j)  # 代表这个人选择接受服务
    return guests[:, :total_ser+1]

5-2 随机游走

一维随机游走的模拟：注意如果将一秒划分成\(n\)刻，则每一步行走的距离是\(\sqrt{n}\)刻。

def rw_1d(t, N, seg):
    # 这里seg表示每一秒分成多少刻
    RW = np.zeros(N)
    number = int(t*seg)
    for i in range(N):
        step = 2*(npr.rand(number)<0.5)-1
        RW[i] = 1/np.sqrt(seg)*np.sum(step)
    return RW

二维随机游走轨迹的模拟：假设每秒粒子在\(X\)方向与\(Y\)方向都按照\(N(0,1)\)行走。

def rw_2d_trace(N):
    xlst = np.zeros(N)
    ylst = np.zeros(N)
    for i in range(N):
        xlst[i] = xlst[i-1] + npr.normal()
        ylst[i] = ylst[i-1] + npr.normal()
    return xlst, ylst

二维随机游走终点分布的模拟：假设同上。

def rw_2d(N, time):
    xlst = np.zeros(N)
    ylst = np.zeros(N)
    for i in range(N):
        for j in range(time):
            xlst[i] += npr.normal()
            ylst[i] += npr.normal()
    return xlst, ylst

布朗运动的标度不变性：在大部分的时间步中，其位移长度大致等于标准差，看上去像聚集在一个区域里；但偶尔会出现一个较大的位移，造成聚集的分区现象。

布朗运动的常返性：一维与二维的对称随机游走具有常返性，而对应布朗运动也具有常返性，这指的是\(\forall\varepsilon > 0\)，布朗运动总会以概率1返回原点的\(\varepsilon\)邻域内。

股票价格模拟：股价服从对数正态分布（伊藤过程），即股价的回报率服从正态分布：

\[\frac{S_t-S_0}{S_0}\sim N(\mu_t,\sigma_t),\quad \forall t. \]

这表明

\[\frac{\Delta S_t}{S_t}=\mu\Delta t+\sigma\Delta W_t. \]

这里：

\(\mu\)和\(\sigma\)是股价变化的特征参数，分别称为漂移率和波动率。
\(\mu\)虽然难以估计，但是在期权定价中不起作用，所以通常以无风险利率替代。
波动率\(\sigma\)反映了股价的波动幅度。
\(\Delta W_t\)是布朗运动增量，独立服从\(N(0, \Delta t)\)。
需要注意\(\mu\)和\(\sigma\)的单位，在金融上常常按年给出，但在模拟时要用日参数，转换公式为

\[\mu_{day} = \frac{\mu_{year}}{365},\quad \sigma_{day}=\frac{\sigma_{year}}{\sqrt{365}}. \]

股价模拟：90天，当前价格20元，波动率为每年0.6，无风险利率为每年0.031。

def stock(times, S0=20, mu=0.031, sigma=0.6, days=90):
    muDay = mu / 365
    sigmaDay = sigma / np.sqrt(365)
    stockPrice = np.zeros(times)
    for i in range(times):
        S = S0
        Wt = npr.normal(0, 1, days)
        for j in range(days):
            dS = S*(muDay + sigmaDay*Wt[j])
            S += dS
        stockPrice[i] = S
    return stockPrice

在此基础上可以获得股票到期日时的价格分布，因而可以获得股价期望\(S_T\)。如果欧式看涨期权的执行价格为\(K\)，则期权的预期收益为\(S_T-K\)。如果无风险利率为\(r_f\)，则期权预期收益的贴现额是

\[(S_T-K)e^{-r_fT}. \]

这就是期权的定价，可以用模拟得到的结果与Black-Scholes期权定价公式进行对比。

凯利问题：在简单赌博\(G\)中，每一次获胜的概率是\(p\)，能赢下下注额\(\gamma\)倍的钱，输的概率为\(q\)，其数学期望

\[\mathbb{E}(G)=p\gamma + q>0, \]

因此是一个长期来看必胜的游戏。凯利游戏想要最快地增长财富，求每一次下注的投注比例。数学推导得到的最佳投注比例是

\[f=\frac{E(G)}{\gamma}=\frac{p\gamma -q}{\gamma}. \]

模拟凯利游戏：胜率为\(p=0.6\)的赌博，赔率为\(1:1\)，模拟最佳下注比例与其他比例的区别。

def Kary(times, p):
    S0 = 10
    win_rate = 0.6
    Cash = np.zeros(times)
    Cash[0] = S0
    result = 2*(npr.rand(times) < win_rate) - 1
    for i in range(1, times):
        Cash[i] = Cash[i-1] + Cash[i-1]*p*result[i]
    return Cash

标签：mathbb,抽样,frac,备考,纲要,np,npr,模拟,guests
来源： https://www.cnblogs.com/jy333/p/14331059.html