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

计算机模拟备考纲要

作者:互联网

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

大纲如下:

  1. 随机数的生成

    • 均匀分布随机数
    • 离散随机数生成
    • 逆变换法
    • 拒绝接受法
    • 正态分布抽样
  2. MCMC

    • Metropolis-Hastings抽样法
    • 概率转移矩阵的选择
    • Ising模型
    • Metropolis算法的原理
  3. MC积分

    • 随机投点法
    • 重要性抽样法
  4. 最优化算法

    • 模拟退火法
    • 遗传算法
  5. 实例模拟

    • 服务系统
    • 随机游走

使用的库:

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)\)中抽样,步骤如下:

  1. 积分\(p(x)\)得到\(F(x)\),求出反函数\(F^{-1}(x)\);
  2. 生成01均匀分布随机数\(u\sim U(0,1)\);
  3. 得到随机数\(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

逆变换法的优点是简单快捷,缺点有:

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

如果生成的均匀分布含边界值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)\)。拒绝接受法的算法步骤如下:

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

    \[\delta(z)=\frac{f(z)}{Mh(z)}. \]

  2. 生成\(U\sim U(0,1)\)的随机数\(u\);

  3. 如果\(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抽样法的优点是:不需要分布函数的反函数。缺点是:

  1. 需要舍弃许多随机数\(z\),尤其是\(z\)发生在\(\delta(z)\)很小的地方时;
  2. 要能找到能完全覆盖\(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})\),其步骤为

  1. 获得初始状态\(x_0\);
  2. 基于\(x_t\)有条件分布\(g(\cdot|x_t)\),从\(g(\cdot|x_t)\)中产生建议随机数\(y\);
  3. 基于\(y\)产生一个接受概率\(h(x_t,y)\),产生01随机数\(r\);
  4. 如果\(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链。

关于对称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\)存在大量的局部最优点。

典型非线性函数: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\)的最小值点。

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

  1. 初始温度\(T_0\)。初始函数\(T_0\)要足够大,否则容易落入局部最优。

  2. 温度\(T\)的衰减函数,即每一次降温过程是如何发生的,下面是两种常用衰减函数:

    \[T_{k+1}=\alpha T_{k}; \\ T_{k+1}=\frac{c}{a+\ln T_k}. \]

    衰减函数的衰减速率不宜过大。

  3. 每一温度下的迭代次数\(L_k\),这相当于Markov链的长度;如果温度衰减足够快,有时可以省略这个参数。

  4. 终止温度\(T_f\),要足够小。

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

  1. 确定冷却进度表,生成初始解。

  2. 在每一个温度下进行\(L_K\)次迭代:

    1. 基于原解的邻域产生一个新解;

    2. 计算接受概率,决定是否接受新解。

      \[h = \min\{1, \exp[(E(S')-E(S))/T] \}. \]

  3. 达到\(L_k\)次迭代后降温。

  4. 返回最优解。

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

4-3 遗传算法

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

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

编码方式:

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

全排列编码:如TSP问题。

综合以上要点,给出遗传算法的步骤:

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

Part 5:实例模拟

5-1 随机服务系统

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

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

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

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

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

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

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

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

构造一种较复杂的算例:服务台数量为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. \]

这里:

股价模拟: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