2021-08-11

蒙特卡罗方法

​ 蒙特卡罗（Monte Carlo）方法，也称为计算机随机模拟方法，是一种基于“随机数”的计算方法。这一方法源于美国在第一次世界大战期间研制原子弹的“曼哈顿计划”。它使用随机数来进行场景的模拟或者过程的仿真，其思想核心就是通过模拟出来的大量样本集或者随机过程去近似我们想要研究的实际问题对象，这是一类非常重要的数值计算方法。

​ 在实际问题中，面对一些带随机因素的复杂系统，用分析方法建模常常需要作许多简化假设，与面临的实际问题可能相差甚远，以致解答根本无法应用。这时，计算机模拟几乎成为唯一的选择。

一、随机变量的模拟

1. 离散随机变量的模拟

​ 设随机变量 $X$ 的分布律为 P { X = x i } = p i , i = 1 , 2 , ⋯ P\left\{X=x_{i}\right\}=p_{i}, i=1,2, \cdots P{X=xi​}=pi​,i=1,2,⋯ , 令
$$p^{(0)}=0,p^{(i)}=\sum _{j=1}^i{p}_j,i=1,2,\cdots$$
​ 将 ${\mathbf{p}}^{(i)}$ 作为分点, 将区间 $(0,1)$ 分为一系列小区间 $\left({\mathbf{p}}^{(i-1)},{\mathbf{p}}^{(i)}\right)(\mathbf{i}=1,2,\cdots )$ 。 对于均匀分布的随机变量 $\mathbf{R}\sim \mathbf{U}(0,1)$, 则有
$$P\left\{p^{(i-1)}<R\le p^{(i)}\right\}=p^{(i)}-p^{(i-1)}=p_i,i=1,2,\cdots ,$$
​ 由此可知，事件 $\left\{{\mathbf{p}}^{(i-1)}<\mathbf{R}\le {\mathbf{p}}^{(i)}\right\}$ 和事件 { X = x i } \left\{\boldsymbol{X}=\boldsymbol{x}_{i}\right\} {X=xi​} 有相同的发生概率。因此我们可以用随机变量 $\mathbf{R}$ 落在小区间内的情况来模拟离散的随机变量 $\mathbf{X}$ 的取值情况。

---

​ 例： 已知在一次随机试验中，事件 $A,B,C$ 发生的概率分别为 $0.2,0.3,0.5$ 。模拟 100,000 次随机试验, 计算事件 $\mathbf{A},\mathbf{B},\mathbf{C}$ 发生的频率。在一次随机试验中，事件发生的概率如下表。
$$\begin{array}{cccc}\text{ 事件 }& \mathbf{A}& \mathbf{B}& \mathbf{C}\\ \text{ 概率 }& 0.2& 0.3& 0.5\\ \text{ 累积概率 }& 0.2& 0.5& 1\end{array}$$
​ 我们用产生[0,1]区间上均匀分布的随机数, 来模拟事件 $\mathbf{A},\mathbf{B},\mathbf{C}$ 的发生。由几何概率的知识，可以认为如果产生的随机数在区间[0,0.2]上, 事件 $A$ 发生了; 产生的随机数在区间 $(0.2,0.5]$ 上, 事件 $B$ 发生了; 产生的随机数在区间 $(0.5,1]$ 上, 事件 $C$ 发生了。产生 100000 个[0,1]区间上均匀分布的随机数, 统计随机数落在相应区间上的次数, 就是在 这 100000 次模拟中事件 $A,B,C$ 发生的次数,再除以总的试验次数 100000,即得到事件 $\mathbf{A},\mathbf{B},\mathbf{C}$ 发生的频率。

代码

import numpy as np
n = 100000
a = np.random.rand(n) #产生[0,1)的均匀分布随机数
n1 = np.sum(a<=0.2)
n2 = np.sum((a>0.2)&(a<=0.5))
n3 = np.sum(a>0.5)
f = np.array([n1,n2,n3])/100000
print(f)

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-in9ErdrV-1628690832671)(C:\Users\haz\AppData\Roaming\Typora\typora-user-images\image-20210811203815021.png)]

2. 连续随机变量的模拟

​ 利用[0,1]区间上的均匀分布随机数可以产生具有给定分布的随机变量数列。我们知道, 若随机变量 $\xi$ 的概率密度函数和分布函数分别为 $f(x),F(x)$ ，则随机变量 $\eta =F(\xi )$ 的分布就是区间 $[0,1]$ 上的均匀分布。因此, 若 $R_i$ 是 $[0,1]$中圴匀分布的随机数，那么方程
$${\int }_{-\infty }^{x_i}f(x)dx=R_i$$
的解 ${\mathbf{x}}_i$ 就是所求的具有概率密度函数为 $f(x)$ 的随机抽样。这可简单解释如下。

若某个连续型随机变量 $\xi$ 的分布函数为
$$F(x)={\int }_{-\infty }^{x}f(t)dt,$$
不失一般性, 设 $F(x)$ 是严格单调增函数, 存在反函数 $x={\mathbf{F}}^{-1}(\mathbf{y})$, 下面我们证明随机变量 $\eta =F(\xi )$ 服从 $[0,1]$ 上的均匀分布, 记 $\eta$ 的分布函数为 $G(y)$, 由于 $F(x)$ 是分布函数, 它的取值在 $[0,1]$ 上，从而当 $0<y<1$ 时
G ( y ) = P { η ≤ y } = P { F ( ξ ) ≤ y } = P { ξ ≤ F − 1 ( y ) } = F ( F − 1 ( y ) ) = y G(y)=P\{\eta \leq y\}=P\{F(\xi) \leq y\}=P\left\{\xi \leq F^{-1}(y)\right\}=F\left(F^{-1}(y)\right)=y G(y)=P{η≤y}=P{F(ξ)≤y}=P{ξ≤F−1(y)}=F(F−1(y))=y
因而 $\eta$ 的分布函数为
$$G(y)=\{\begin{array}{ll}0,& y\le 0\\ y,& 0<y<1\\ 1,& y\ge 1\end{array}$$
$\eta$ 服从[0,1]上的均匀分布。 $R$ 为 $[0,1]$ 区间上均匀分布的随机变量, 则随机变量 $\xi =F^{-1}(R)$ 的分布函数 为 $F(x)$, 概率密度函数为 $f(x)$, 这里 $F^{-1}(x)$ 是 $F(x)$ 的反函数。

所以, 只要分布函数 $F(x)$ 的反函数 $F^{-1}(x)$ 存在, 由[0,1]区间上均匀分布
的随机数 ${\mathbf{R}}_t$, 求 ${\mathbf{x}}_t={\mathbf{F}}^{-1}\left({\mathbf{R}}_t\right)$, 即解方程
$$F\left(x_t\right)=R_t,$$
就可得到分布函数为 $F(x)$ 的随机抽样 $x_t$ 。

---

例： 求具有指数分布
$$f(x)=\{\begin{array}{ll}\lambda e^{-\lambda x},& x>0\\ 0,& x\le 0\end{array}$$
的随机抽样。
设 ${\mathbf{R}}_i$ 是[0,1]区间上均匀分布的随机数, 则
$$R_i={\int }_{-\infty }^{x_i}f(x)dx={\int }_0^{x_i}\lambda e^{-\lambda x}dx=1-e^{-\lambda x_i}$$
所以
$$x_i=-\frac{1}{\lambda }\ln \left(1-R_i\right)$$
就是所求的随机抽样。

由于 $1-{\mathbf{R}}_i$ 也服从均匀分布，所以上式又可简化为
$$x_i=-\frac{1}{\lambda }\ln R_i$$
Python 中 $NumPy$ 库函数可以产生常用分布的随机数，实际上不需要 我们去生成各种分布的随机数。

二、$MonteCarlo$法的应用

1.定积分的计算

​ 在实际问题中，许多需要计算多重积分的复杂问题，用蒙特卡洛方法一般都能够很有效地予以解决，尽管蒙特卡洛方法计算结果的精度不很高，但它能很快提供出一个低精度的模拟结果也是很有价值的。而且, 在多重积分中蒙特卡洛方法的计算误差与积分重数无关。

例1： 求定积分

$$I={\iint }_{x^2+y^2\le 1}\sqrt{1-x^2}dxdy$$

根据积分的几何意义，它的值等于以$x^2+z^2=1(z>0)$为曲面顶，$x^2+y^2\le 1,z=0$为底的柱体$D$的体积，用以下思路求$I$的近似值：

$D$被包含在长方体$\Omega$：$|x|\le 1,|y|\le 1,0\le z\le 1$ 的内部，长方体 $\Omega$ 的体积为 4 。 而 $D$ 是 $x^2+y^2\le 1$, $0\le \mathbf{z}\le \sqrt{1-{\mathbf{x}}^2}$ 所围成的区域。若在 $\Omega$ 内产生均匀分布的 $N$ 个点, 有 $n$ 个点落 在D的内部。由频率近似于概率, 得到在 $\Omega$ 内任取一点, 落在 $D$ 内的概率 $p=\frac{D\text{ 的体积 }}{\Omega \text{ 的体积 }}=\frac{I}{4}\approx \frac{n}{N}$, 所以 $I\approx \frac{4n}{N}$ 。

代码

import numpy as np
N = 10000000
x = np.random.uniform(-1,1,size=N)
y = np.random.uniform(-1,1,N)
z = np.random.uniform(0,1,N)
n = np.sum((x**2+y**2<=1)&(0<=z)&(z<=np.sqrt(1-x**2)))
I = n/N*4
print('I近似为',I)

2.概率在计算中的应用

例1：求 $\pi$ 的近似值

在单位正方形$(x\in [0,1],y\in [0,1])$中取 10,000,000 个随机点 $\left({\mathbf{x}}_i,{\mathbf{y}}_i\right)$, $\mathbf{i}=1,2,\cdots ,10000000$, 统计点落在 ${\mathbf{x}}^2+{\mathbf{y}}^2\le 1$ 内的频数 $\mathbf{n}$ 。则由几何概率知, 任取单位正方形内一点, 落在单位圆内部（第一象限部分）的概率为 $p=\frac{\pi }{4}$,由于试验次数充分多,频率近似于概率,有 $\frac{n}{10000000}\approx \frac{\pi }{4}$, 所以 $\pi \approx \frac{4n}{10000000}$ 。

代码

import numpy as np
N = 10000000
x = np.random.rand(N)
y = np.random.rand(N)
n = np.sum(x**2+y**2<1)
PI = n/N*4
print(PI)

例2：炮弹射击

炮弹射击的目标为一椭圆 $\frac{x^2}{120^2}+\frac{y^2}{80^2}=1$ 所围成的区域的中心, 瞄准目标的中心发射时, 受到各种因素的影响, 炮弹着地点与目标中心有
随机偏差。设炮弹着地点围绕目标中心呈二维正态分布, 且偏差的标准差在$x$ 和 $y$ 方向均为 100 米，并相互独立, 用 Monte Carlo 法计算炮弹落在椭圆区域内的概率，并与数值积分计算的概率进行比较。

解: 炮弹的落点为二维随机变量, 记为 $(\mathbf{X},\mathbf{Y}),(\mathbf{X},\mathbf{Y})$ 的联合概率密度函数为
$$f(x,y)=\frac{1}{20000\pi }{e}^{-\frac{x^2+y^2}{20000}}$$
炮弹落在椭圆区域内的概率为
$$p={\iint }_{\frac{x^2}{120^2}+\frac{y^2}{80^2}\le 1}\frac{1}{20000\pi }{e}^{-\frac{x^2+y^2}{20000}}dxdy$$
利用 Python 数值解的命令，求得 $\mathbf{p}=0.3754$ 。

我们也可以使用 Monte Carlo 法求概率。模拟发射了 $N$ 发炮弹，统计炮弹落在椭圆 $\frac{x^2}{120^2}+\frac{y^2}{80^2}=1$ 内部的次数 $n$, 用炮弹落在椭圆内的频率近似所求的概率，模拟结果为所求的概率在 $0.3754$ 附近波动。

代码

import numpy as np
from scipy.integrate import dblquad

f = lambda x,y: 1/(20000*np.pi)*np.exp(-(x**2+y**2)/20000)
fy = lambda x: 80*np.sqrt(1-x**2/120**2)
p1 = dblquad(f,-120,120,lambda x: -fy(x),lambda x: fy(x))
print("概率的数值解为：",p1[0])

N = 1000000
mean = [0,0]
cov = np.eye(2)*10000
a = np.random.multivariate_normal(mean,cov,N)
n = (a[:,0]**2/120**2+a[:,1]**2/80**2<=1).sum()
p2 = n/N
print("概率的近似值为",p2)

例3：求全局最优解

求下列函数的最大值:
$$f(x)=\left(1-x^3\right)\sin (3x),-2\pi \le x\le 2\pi .$$
解 为了便于理解，先做出图象

import numpy as np
import matplotlib.pyplot as plt

x = np.linspace(-2*np.pi,2*np.pi,100)
y = (1-x**3)*np.sin(3*x)
plt.plot(x,y)
plt.show()

可见，函数在-6和 6 附近达到最大值,最大值接近 195。

如果使用优化命令 $fminbound$ 求解，只能求得局部极大点 x=-3.7505, 对应的局部极大值为52.0046。显然结果是错误的, 原因是 $fminbound$ 容易䧂入局部极值，这也是许多优化算法难以克服的一个困难。

from scipy import optimize
import numpy as np
f = lambda x: (1-x**3)*np.sin(3*x)
a = optimize.fminbound(lambda x: -f(x),-np.pi*2,np.pi*2) # 标量函数的有界最小化。
print(a,f(a))

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传
用随机模拟的方法，就是随机产生若干个自变量的值来搜索，求得最大值再194附近，效果要好得多。

import numpy as np
f = lambda x: (1-x**3)*np.sin(3*x)
x = np.random.uniform(-2*np.pi,2*np.pi,100)
y = f(x)
y_max = y.max()
x_ = x[f(x) == y_max]
print(x_[0],y_max)

例4：库存管理问题

某小贩每天以 $a=2$ 元/束的价格购进一种鲜花, 卖价为 $b=3$ 元/束，当天卖不出去的花全部损失. 顾客一天内对花的需求量 $X$ 是随机变量,$X$ 服从泊松分布
P { X = k } = e − λ λ k k ! , k = 0 , 1 , 2 , ⋯ P\{X=k\}=e^{-\lambda} \frac{\lambda^{k}}{k !}, k=0,1,2, \cdots P{X=k}=e−λk!λk​,k=0,1,2,⋯
其中参数 $\lambda =10$. 问小贩每天应购进多少束鲜花才能得到好收益?

解 这是一个随机决策问题，要确定每天应购进的鲜花数量以使收入最高。设小贩每天购进 $u$ 束鲜花. 如果这天需求量 $X\le \mathbf{u}$, 则其收入为 $b\mathbf{X}-\mathbf{a}\mathbf{u}$,如果需求量 $\mathbf{X}>\mathbf{u}$, 则其收入为 $\mathbf{b}\mathbf{u}-\mathbf{a}\mathbf{u}$, 因此小贩一天的期望收入为
$$J(u)=-au+\sum _{k=0}^{u}bk\cdot e^{-\lambda }\cdot \frac{{\lambda }^k}{k!}+\sum _{k=u+1}^{\infty }bu\cdot e^{-\lambda }\cdot \frac{{\lambda }^k}{k!}$$
问题归结为在 $a,b,\lambda$ 已知时，求 $u$ 使得 $J(u)$ 最大. 因而最佳购进量 $u^{\ast }$ 满足
$$J\left(u^{\ast }\right)\ge J\left(u^{\ast }+1\right),J\left(u^{\ast }\right)\ge J\left(u^{\ast }-1\right),$$
由于 $J(u+1)-J(u)=-a+b{e}^{-\lambda }{\sum }_{k=u+1}^{\infty }\frac{{\lambda }^k}{k!}=-a+b\left(1-{\sum }_{k=0}^u{e}^{-\lambda }\frac{{\lambda }^k}{k!}\right)$

最佳购进量 ${\mathbf{u}}^{\ast }$ 满足
$$\begin{array}{rl}& 1-\sum _{k=0}^{u^{\ast }}{e}^{-\lambda }\frac{{\lambda }^k}{k!}\le \frac{a}{b}\\ & 1-\sum _{k=0}^{u^{\ast }-1}{e}^{-\lambda }\frac{{\lambda }^k}{k!}\ge \frac{a}{b}\end{array}$$
记泊松分布的分布函数为 F ( i ) = P { X ≤ i } = ∑ k = 0 i e − λ λ k k ! F(i)=P\{X \leq i\}=\sum_{k=0}^{i} e^{-\lambda} \frac{\lambda^{k}}{k !} F(i)=P{X≤i}=∑k=0i​e−λk!λk​, 则最佳购进量 $u^{\ast }$ 满足
$$F\left(u^{\ast }-1\right)\le 1-\frac{a}{b}\le F\left(u^{\ast }\right)$$
查 Poisson 分布表, 或利用 Python 软件, 求得最佳购进量 ${\mathbf{u}}^{\ast }=9$ 。

from scipy.stats import poisson
a=2; b=3
lamda =10; p=1-a/b
u=poisson.ppf(1-a/b,lamda) # 求最佳订购量
p1=poisson.cdf(u-1,lamda) #p1和p2是为验证最佳购进量
p2=poisson.cdf(u,lamda)
print(u,p1,p,p2)

下面用计算机模拟进行检验。 对不同的 $\mathbf{a},\mathbf{b},\lambda$ ，用计算机模拟求最优决策 $u$ 的算法如下:

• 步骤 1 给定 $\mathbf{a},\mathbf{b},\lambda$, 记进货量为 $u$ 时, 收益为 ${\mathbf{M}}_u$, 当 $\mathbf{u}=0$ 时, ${\mathbf{M}}_0=0$;令 $\mathbf{u}=1$, 继续下一步。
• 步骤 2 对随机需求变量 $X$ 做模拟，求出收入，共做 $n$ 次模拟，求出收入的平均值 $M_u$ 。
• 步骤 3 若 ${\mathbf{M}}_u\ge {\mathbf{M}}_{u-1}$, 令 $\mathbf{u}=\mathbf{u}+1$, 转步骤 $2;$ 若 ${\mathbf{M}}_u<{\mathbf{M}}_{u-1}$, 输出 ${\mathbf{u}}^{\ast }=\mathbf{u}-1$, 停止。

用 Python 软件进行模拟，求得最佳进货量为 8 或 9 , 发现其与理论推导符合得很好。模拟的 Python 程序如下:

import numpy as npa=2; b=3; lamda=10; M1=0;u=1; n=10000;for i in range(1,2*lamda):    d=np.random.poisson(lamda,n)  #产生n个服从Poiss分布的需求量数据    M2=np.mean(((b-a)*u*(u<=d)+((b-a)*d-a*(u-d))*(u>d)))  #求平均利润    if M2>M1: M1=M2; u=u+1;    else:        print('最佳购进量:',u-1);        break

三、概率论相关知识补充

超几何分布

设有N件产品，其中有M件次品，现从中任取n件，用X表示其中的次品数，求其分布律。
$$P(X=k)=\frac{C_M^k{C}_{N-M}^{n-k}}{C_N^n},k=0,1,2,\cdots ,l.$$

$$l=min(M,n)$$

泊松分布

若
$$P(X=k)=e^{-\lambda }\frac{{\lambda }^k}{k!},k=0,1,2,\cdots$$
其中 $\lambda >0$ 是常数，则称 $X$ 服从参数为 $\lambda$ 的泊松 (Poisson) 分布. 记作 $X\sim \pi (\lambda )$ 或 $P(\lambda )$

应用场合 在某个时段内：

• 市级医院急诊病人数;
• 某地区拨错号的电话呼唤次数;
• 某地区发生的交通事故的次数. # 蒙特卡罗方法

​ 蒙特卡罗（Monte Carlo）方法，也称为计算机随机模拟方法，是一种基于“随机数”的计算方法。这一方法源于美国在第一次世界大战期间研制原子弹的“曼哈顿计划”。它使用随机数来进行场景的模拟或者过程的仿真，其思想核心就是通过模拟出来的大量样本集或者随机过程去近似我们想要研究的实际问题对象，这是一类非常重要的数值计算方法。

​ 在实际问题中，面对一些带随机因素的复杂系统，用分析方法建模常常需要作许多简化假设，与面临的实际问题可能相差甚远，以致解答根本无法应用。这时，计算机模拟几乎成为唯一的选择。

一、随机变量的模拟

1. 离散随机变量的模拟

​ 设随机变量 $X$ 的分布律为 P { X = x i } = p i , i = 1 , 2 , ⋯ P\left\{X=x_{i}\right\}=p_{i}, i=1,2, \cdots P{X=xi​}=pi​,i=1,2,⋯ , 令
$$p^{(0)}=0,p^{(i)}=\sum _{j=1}^i{p}_j,i=1,2,\cdots$$
​ 将 ${\mathbf{p}}^{(i)}$ 作为分点, 将区间 $(0,1)$ 分为一系列小区间 $\left({\mathbf{p}}^{(i-1)},{\mathbf{p}}^{(i)}\right)(\mathbf{i}=1,2,\cdots )$ 。 对于均匀分布的随机变量 $\mathbf{R}\sim \mathbf{U}(0,1)$, 则有
$$P\left\{p^{(i-1)}<R\le p^{(i)}\right\}=p^{(i)}-p^{(i-1)}=p_i,i=1,2,\cdots ,$$
​ 由此可知，事件 $\left\{{\mathbf{p}}^{(i-1)}<\mathbf{R}\le {\mathbf{p}}^{(i)}\right\}$ 和事件 { X = x i } \left\{\boldsymbol{X}=\boldsymbol{x}_{i}\right\} {X=xi​} 有相同的发生概率。因此我们可以用随机变量 $\mathbf{R}$ 落在小区间内的情况来模拟离散的随机变量 $\mathbf{X}$ 的取值情况。

---

​ 例： 已知在一次随机试验中，事件 $A,B,C$ 发生的概率分别为 $0.2,0.3,0.5$ 。模拟 100,000 次随机试验, 计算事件 $\mathbf{A},\mathbf{B},\mathbf{C}$ 发生的频率。在一次随机试验中，事件发生的概率如下表。
$$\begin{array}{cccc}\text{ 事件 }& \mathbf{A}& \mathbf{B}& \mathbf{C}\\ \text{ 概率 }& 0.2& 0.3& 0.5\\ \text{ 累积概率 }& 0.2& 0.5& 1\end{array}$$
​ 我们用产生[0,1]区间上均匀分布的随机数, 来模拟事件 $\mathbf{A},\mathbf{B},\mathbf{C}$ 的发生。由几何概率的知识，可以认为如果产生的随机数在区间[0,0.2]上, 事件 $A$ 发生了; 产生的随机数在区间 $(0.2,0.5]$ 上, 事件 $B$ 发生了; 产生的随机数在区间 $(0.5,1]$ 上, 事件 $C$ 发生了。产生 100000 个[0,1]区间上均匀分布的随机数, 统计随机数落在相应区间上的次数, 就是在 这 100000 次模拟中事件 $A,B,C$ 发生的次数,再除以总的试验次数 100000,即得到事件 $\mathbf{A},\mathbf{B},\mathbf{C}$ 发生的频率。

代码

import numpy as np
n = 100000
a = np.random.rand(n) #产生[0,1)的均匀分布随机数
n1 = np.sum(a<=0.2)
n2 = np.sum((a>0.2)&(a<=0.5))
n3 = np.sum(a>0.5)
f = np.array([n1,n2,n3])/100000
print(f)

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

2. 连续随机变量的模拟

​ 利用[0,1]区间上的均匀分布随机数可以产生具有给定分布的随机变量数列。我们知道, 若随机变量 $\xi$ 的概率密度函数和分布函数分别为 $f(x),F(x)$ ，则随机变量 $\eta =F(\xi )$ 的分布就是区间 $[0,1]$ 上的均匀分布。因此, 若 $R_i$ 是 $[0,1]$中圴匀分布的随机数，那么方程
$${\int }_{-\infty }^{x_i}f(x)dx=R_i$$
的解 ${\mathbf{x}}_i$ 就是所求的具有概率密度函数为 $f(x)$ 的随机抽样。这可简单解释如下。

若某个连续型随机变量 $\xi$ 的分布函数为
$$F(x)={\int }_{-\infty }^{x}f(t)dt,$$
不失一般性, 设 $F(x)$ 是严格单调增函数, 存在反函数 $x={\mathbf{F}}^{-1}(\mathbf{y})$, 下面我们证明随机变量 $\eta =F(\xi )$ 服从 $[0,1]$ 上的均匀分布, 记 $\eta$ 的分布函数为 $G(y)$, 由于 $F(x)$ 是分布函数, 它的取值在 $[0,1]$ 上，从而当 $0<y<1$ 时
G ( y ) = P { η ≤ y } = P { F ( ξ ) ≤ y } = P { ξ ≤ F − 1 ( y ) } = F ( F − 1 ( y ) ) = y G(y)=P\{\eta \leq y\}=P\{F(\xi) \leq y\}=P\left\{\xi \leq F^{-1}(y)\right\}=F\left(F^{-1}(y)\right)=y G(y)=P{η≤y}=P{F(ξ)≤y}=P{ξ≤F−1(y)}=F(F−1(y))=y
因而 $\eta$ 的分布函数为
$$G(y)=\{\begin{array}{ll}0,& y\le 0\\ y,& 0<y<1\\ 1,& y\ge 1\end{array}$$
$\eta$ 服从[0,1]上的均匀分布。 $R$ 为 $[0,1]$ 区间上均匀分布的随机变量, 则随机变量 $\xi =F^{-1}(R)$ 的分布函数 为 $F(x)$, 概率密度函数为 $f(x)$, 这里 $F^{-1}(x)$ 是 $F(x)$ 的反函数。

所以, 只要分布函数 $F(x)$ 的反函数 $F^{-1}(x)$ 存在, 由[0,1]区间上均匀分布
的随机数 ${\mathbf{R}}_t$, 求 ${\mathbf{x}}_t={\mathbf{F}}^{-1}\left({\mathbf{R}}_t\right)$, 即解方程
$$F\left(x_t\right)=R_t,$$
就可得到分布函数为 $F(x)$ 的随机抽样 $x_t$ 。

---

例： 求具有指数分布
$$f(x)=\{\begin{array}{ll}\lambda e^{-\lambda x},& x>0\\ 0,& x\le 0\end{array}$$
的随机抽样。
设 ${\mathbf{R}}_i$ 是[0,1]区间上均匀分布的随机数, 则
$$R_i={\int }_{-\infty }^{x_i}f(x)dx={\int }_0^{x_i}\lambda e^{-\lambda x}dx=1-e^{-\lambda x_i}$$
所以
$$x_i=-\frac{1}{\lambda }\ln \left(1-R_i\right)$$
就是所求的随机抽样。

由于 $1-{\mathbf{R}}_i$ 也服从均匀分布，所以上式又可简化为
$$x_i=-\frac{1}{\lambda }\ln R_i$$
Python 中 $NumPy$ 库函数可以产生常用分布的随机数，实际上不需要 我们去生成各种分布的随机数。

二、$MonteCarlo$法的应用

1.定积分的计算

​ 在实际问题中，许多需要计算多重积分的复杂问题，用蒙特卡洛方法一般都能够很有效地予以解决，尽管蒙特卡洛方法计算结果的精度不很高，但它能很快提供出一个低精度的模拟结果也是很有价值的。而且, 在多重积分中蒙特卡洛方法的计算误差与积分重数无关。

例1： 求定积分

$$I={\iint }_{x^2+y^2\le 1}\sqrt{1-x^2}dxdy$$

根据积分的几何意义，它的值等于以$x^2+z^2=1(z>0)$为曲面顶，$x^2+y^2\le 1,z=0$为底的柱体$D$的体积，用以下思路求$I$的近似值：

$D$被包含在长方体$\Omega$：$|x|\le 1,|y|\le 1,0\le z\le 1$ 的内部，长方体 $\Omega$ 的体积为 4 。 而 $D$ 是 $x^2+y^2\le 1$, $0\le \mathbf{z}\le \sqrt{1-{\mathbf{x}}^2}$ 所围成的区域。若在 $\Omega$ 内产生均匀分布的 $N$ 个点, 有 $n$ 个点落 在D的内部。由频率近似于概率, 得到在 $\Omega$ 内任取一点, 落在 $D$ 内的概率 $p=\frac{D\text{ 的体积 }}{\Omega \text{ 的体积 }}=\frac{I}{4}\approx \frac{n}{N}$, 所以 $I\approx \frac{4n}{N}$ 。

代码

import numpy as np
N = 10000000
x = np.random.uniform(-1,1,size=N)
y = np.random.uniform(-1,1,N)
z = np.random.uniform(0,1,N)
n = np.sum((x**2+y**2<=1)&(0<=z)&(z<=np.sqrt(1-x**2)))
I = n/N*4
print('I近似为',I)

2.概率在计算中的应用

例1：求 $\pi$ 的近似值

在单位正方形$(x\in [0,1],y\in [0,1])$中取 10,000,000 个随机点 $\left({\mathbf{x}}_i,{\mathbf{y}}_i\right)$, $\mathbf{i}=1,2,\cdots ,10000000$, 统计点落在 ${\mathbf{x}}^2+{\mathbf{y}}^2\le 1$ 内的频数 $\mathbf{n}$ 。则由几何概率知, 任取单位正方形内一点, 落在单位圆内部（第一象限部分）的概率为 $p=\frac{\pi }{4}$,由于试验次数充分多,频率近似于概率,有 $\frac{n}{10000000}\approx \frac{\pi }{4}$, 所以 $\pi \approx \frac{4n}{10000000}$ 。

代码

import numpy as np
N = 10000000
x = np.random.rand(N)
y = np.random.rand(N)
n = np.sum(x**2+y**2<1)
PI = n/N*4
print(PI)

例2：炮弹射击

炮弹射击的目标为一椭圆 $\frac{x^2}{120^2}+\frac{y^2}{80^2}=1$ 所围成的区域的中心, 瞄准目标的中心发射时, 受到各种因素的影响, 炮弹着地点与目标中心有
随机偏差。设炮弹着地点围绕目标中心呈二维正态分布, 且偏差的标准差在$x$ 和 $y$ 方向均为 100 米，并相互独立, 用 Monte Carlo 法计算炮弹落在椭圆区域内的概率，并与数值积分计算的概率进行比较。

解: 炮弹的落点为二维随机变量, 记为 $(\mathbf{X},\mathbf{Y}),(\mathbf{X},\mathbf{Y})$ 的联合概率密度函数为
$$f(x,y)=\frac{1}{20000\pi }{e}^{-\frac{x^2+y^2}{20000}}$$
炮弹落在椭圆区域内的概率为
$$p={\iint }_{\frac{x^2}{120^2}+\frac{y^2}{80^2}\le 1}\frac{1}{20000\pi }{e}^{-\frac{x^2+y^2}{20000}}dxdy$$
利用 Python 数值解的命令，求得 $\mathbf{p}=0.3754$ 。

我们也可以使用 Monte Carlo 法求概率。模拟发射了 $N$ 发炮弹，统计炮弹落在椭圆 $\frac{x^2}{120^2}+\frac{y^2}{80^2}=1$ 内部的次数 $n$, 用炮弹落在椭圆内的频率近似所求的概率，模拟结果为所求的概率在 $0.3754$ 附近波动。

代码

import numpy as np
from scipy.integrate import dblquad

f = lambda x,y: 1/(20000*np.pi)*np.exp(-(x**2+y**2)/20000)
fy = lambda x: 80*np.sqrt(1-x**2/120**2)
p1 = dblquad(f,-120,120,lambda x: -fy(x),lambda x: fy(x))
print("概率的数值解为：",p1[0])

N = 1000000
mean = [0,0]
cov = np.eye(2)*10000
a = np.random.multivariate_normal(mean,cov,N)
n = (a[:,0]**2/120**2+a[:,1]**2/80**2<=1).sum()
p2 = n/N
print("概率的近似值为",p2)

例3：求全局最优解

求下列函数的最大值:
$$f(x)=\left(1-x^3\right)\sin (3x),-2\pi \le x\le 2\pi .$$
解 为了便于理解，先做出图象

import numpy as np
import matplotlib.pyplot as plt

x = np.linspace(-2*np.pi,2*np.pi,100)
y = (1-x**3)*np.sin(3*x)
plt.plot(x,y)
plt.show()

可见，函数在-6和 6 附近达到最大值,最大值接近 195。

如果使用优化命令 $fminbound$ 求解，只能求得局部极大点 x=-3.7505, 对应的局部极大值为52.0046。显然结果是错误的, 原因是 $fminbound$ 容易䧂入局部极值，这也是许多优化算法难以克服的一个困难。

from scipy import optimize
import numpy as np
f = lambda x: (1-x**3)*np.sin(3*x)
a = optimize.fminbound(lambda x: -f(x),-np.pi*2,np.pi*2) # 标量函数的有界最小化。
print(a,f(a))

用随机模拟的方法，就是随机产生若干个自变量的值来搜索，求得最大值再194附近，效果要好得多。

import numpy as np
f = lambda x: (1-x**3)*np.sin(3*x)
x = np.random.uniform(-2*np.pi,2*np.pi,100)
y = f(x)
y_max = y.max()
x_ = x[f(x) == y_max]
print(x_[0],y_max)

例4：库存管理问题

某小贩每天以 $a=2$ 元/束的价格购进一种鲜花, 卖价为 $b=3$ 元/束，当天卖不出去的花全部损失. 顾客一天内对花的需求量 $X$ 是随机变量,$X$ 服从泊松分布
P { X = k } = e − λ λ k k ! , k = 0 , 1 , 2 , ⋯ P\{X=k\}=e^{-\lambda} \frac{\lambda^{k}}{k !}, k=0,1,2, \cdots P{X=k}=e−λk!λk​,k=0,1,2,⋯
其中参数 $\lambda =10$. 问小贩每天应购进多少束鲜花才能得到好收益?

解 这是一个随机决策问题，要确定每天应购进的鲜花数量以使收入最高。设小贩每天购进 $u$ 束鲜花. 如果这天需求量 $X\le \mathbf{u}$, 则其收入为 $b\mathbf{X}-\mathbf{a}\mathbf{u}$,如果需求量 $\mathbf{X}>\mathbf{u}$, 则其收入为 $\mathbf{b}\mathbf{u}-\mathbf{a}\mathbf{u}$, 因此小贩一天的期望收入为
$$J(u)=-au+\sum _{k=0}^{u}bk\cdot e^{-\lambda }\cdot \frac{{\lambda }^k}{k!}+\sum _{k=u+1}^{\infty }bu\cdot e^{-\lambda }\cdot \frac{{\lambda }^k}{k!}$$
问题归结为在 $a,b,\lambda$ 已知时，求 $u$ 使得 $J(u)$ 最大. 因而最佳购进量 $u^{\ast }$ 满足
$$J\left(u^{\ast }\right)\ge J\left(u^{\ast }+1\right),J\left(u^{\ast }\right)\ge J\left(u^{\ast }-1\right),$$
由于 $J(u+1)-J(u)=-a+b{e}^{-\lambda }{\sum }_{k=u+1}^{\infty }\frac{{\lambda }^k}{k!}=-a+b\left(1-{\sum }_{k=0}^u{e}^{-\lambda }\frac{{\lambda }^k}{k!}\right)$

最佳购进量 ${\mathbf{u}}^{\ast }$ 满足
$$\begin{array}{rl}& 1-\sum _{k=0}^{u^{\ast }}{e}^{-\lambda }\frac{{\lambda }^k}{k!}\le \frac{a}{b}\\ & 1-\sum _{k=0}^{u^{\ast }-1}{e}^{-\lambda }\frac{{\lambda }^k}{k!}\ge \frac{a}{b}\end{array}$$
记泊松分布的分布函数为 F ( i ) = P { X ≤ i } = ∑ k = 0 i e − λ λ k k ! F(i)=P\{X \leq i\}=\sum_{k=0}^{i} e^{-\lambda} \frac{\lambda^{k}}{k !} F(i)=P{X≤i}=∑k=0i​e−λk!λk​, 则最佳购进量 $u^{\ast }$ 满足
$$F\left(u^{\ast }-1\right)\le 1-\frac{a}{b}\le F\left(u^{\ast }\right)$$
查 Poisson 分布表, 或利用 Python 软件, 求得最佳购进量 ${\mathbf{u}}^{\ast }=9$ 。

from scipy.stats import poisson
a=2; b=3
lamda =10; p=1-a/b
u=poisson.ppf(1-a/b,lamda) # 求最佳订购量
p1=poisson.cdf(u-1,lamda) #p1和p2是为验证最佳购进量
p2=poisson.cdf(u,lamda)
print(u,p1,p,p2)

下面用计算机模拟进行检验。 对不同的 $\mathbf{a},\mathbf{b},\lambda$ ，用计算机模拟求最优决策 $u$ 的算法如下:

• 步骤 1 给定 $\mathbf{a},\mathbf{b},\lambda$, 记进货量为 $u$ 时, 收益为 ${\mathbf{M}}_u$, 当 $\mathbf{u}=0$ 时, ${\mathbf{M}}_0=0$;令 $\mathbf{u}=1$, 继续下一步。
• 步骤 2 对随机需求变量 $X$ 做模拟，求出收入，共做 $n$ 次模拟，求出收入的平均值 $M_u$ 。
• 步骤 3 若 ${\mathbf{M}}_u\ge {\mathbf{M}}_{u-1}$, 令 $\mathbf{u}=\mathbf{u}+1$, 转步骤 $2;$ 若 ${\mathbf{M}}_u<{\mathbf{M}}_{u-1}$, 输出 ${\mathbf{u}}^{\ast }=\mathbf{u}-1$, 停止。

用 Python 软件进行模拟，求得最佳进货量为 8 或 9 , 发现其与理论推导符合得很好。模拟的 Python 程序如下:

import numpy as npa=2; b=3; lamda=10; M1=0;u=1; n=10000;for i in range(1,2*lamda):    d=np.random.poisson(lamda,n)  #产生n个服从Poiss分布的需求量数据    M2=np.mean(((b-a)*u*(u<=d)+((b-a)*d-a*(u-d))*(u>d)))  #求平均利润    if M2>M1: M1=M2; u=u+1;    else:        print('最佳购进量:',u-1);        break

三、概率论相关知识补充

超几何分布

设有N件产品，其中有M件次品，现从中任取n件，用X表示其中的次品数，求其分布律。
$$P(X=k)=\frac{C_M^k{C}_{N-M}^{n-k}}{C_N^n},k=0,1,2,\cdots ,l.$$

$$l=min(M,n)$$

泊松分布

若
$$P(X=k)=e^{-\lambda }\frac{{\lambda }^k}{k!},k=0,1,2,\cdots$$
其中 $\lambda >0$ 是常数，则称 $X$ 服从参数为 $\lambda$ 的泊松 (Poisson) 分布. 记作 $X\sim \pi (\lambda )$ 或 $P(\lambda )$

应用场合 在某个时段内：

• 市级医院急诊病人数;
• 某地区拨错号的电话呼唤次数;
• 某地区发生的交通事故的次数.
• 一本书一页中的印刷错误数.

正态分布
若随机变量 $X$ 的概率密度函数为
$$f(x)=\frac{1}{\sqrt{2\pi }\sigma }{e}^{-\frac{(x-\mu {)}^2}{2{\sigma }^2}},-\infty <x<+\infty$$
其中 $\mu ,\sigma$ 为常数，且 $\sigma >0$, 则称随机变量 $X$ 服从参数 为 $\mu ,{\sigma }^2$ 的正态分布.服从正态分布的随机变量为正态变量. 记作 $X\sim N\left(\mu ,{\sigma }^2\right)$

二维正态分布
若 $r.v.(X,Y)$ 的联合 $d.f.$ 为
$$f(x,y)=\frac{1}{2\pi {\sigma }_1{\sigma }_2\sqrt{1-{\rho }^2}}\times e^{-\frac{1}{2\left(1-{\rho }^2\right)}\left[\frac{{\left(x-{\mu }_1\right)}^2}{{\sigma }_1^2}-2\rho \frac{\left(x-{\mu }_1\right)\left(y-{\mu }_2\right)}{{\sigma }_1{\sigma }_2}+\frac{{\left(y-{\mu }_2\right)}^2}{{\sigma }_2^2}\right]}$$
则称 $(X,Y)$ 服从参数为 ${\mu }_1,{\sigma }_1^2,{\mu }_2,{\sigma }_2^2,\rho$ 的正态分布,记作 $(X,Y)\sim N\left({\mu }_1,{\sigma }_1^2;{\mu }_2,{\sigma }_2^2;\rho \right).$ 其中 $,{\sigma }_1,{\sigma }_2>0,-1<\rho <1.$

• 一本书一页中的印刷错误数.

正态分布
若随机变量 $X$ 的概率密度函数为
$$f(x)=\frac{1}{\sqrt{2\pi }\sigma }{e}^{-\frac{(x-\mu {)}^2}{2{\sigma }^2}},-\infty <x<+\infty$$
其中 $\mu ,\sigma$ 为常数，且 $\sigma >0$, 则称随机变量 $X$ 服从参数 为 $\mu ,{\sigma }^2$ 的正态分布.服从正态分布的随机变量为正态变量. 记作 $X\sim N\left(\mu ,{\sigma }^2\right)$

二维正态分布
若 $r.v.(X,Y)$ 的联合 $d.f.$ 为
$$f(x,y)=\frac{1}{2\pi {\sigma }_1{\sigma }_2\sqrt{1-{\rho }^2}}\times e^{-\frac{1}{2\left(1-{\rho }^2\right)}\left[\frac{{\left(x-{\mu }_1\right)}^2}{{\sigma }_1^2}-2\rho \frac{\left(x-{\mu }_1\right)\left(y-{\mu }_2\right)}{{\sigma }_1{\sigma }_2}+\frac{{\left(y-{\mu }_2\right)}^2}{{\sigma }_2^2}\right]}$$
则称 $(X,Y)$ 服从参数为 ${\mu }_1,{\sigma }_1^2,{\mu }_2,{\sigma }_2^2,\rho$ 的正态分布,记作 $(X,Y)\sim N\left({\mu }_1,{\sigma }_1^2;{\mu }_2,{\sigma }_2^2;\rho \right).$ 其中 $,{\sigma }_1,{\sigma }_2>0,-1<\rho <1.$