Python：利用蒙特卡洛方法模拟验证概率分布

利用 MonteCarlo 模拟验证概率分布

有这样一道题目：

已知两个独立随机变量 $x,y$，随机变量 $x$ 服从几何分布 $\mathrm{G}\mathrm{e}\mathrm{o}\mathrm{m}(p)$，$y$ 服从区间 $[0,1]$ 上的均匀分布 $\mathrm{U}(0,1)$，求新的随机变量 $z=xy$ 的概率分布。

这个题目可以使用数学方法，将其答案显式地写出来，但是验证解出来的答案是否正确，就可以使用蒙特卡洛方法了。我们可以先写出自己的答案，然后编程看看使用蒙特卡洛方法模拟出来的结果与我们自己计算出来的结果是否一致。

1/ 使用数学方法解题

第一步我们先用高数的知识解题，这一步如果看不懂，可以跳过，直接看第二步的编程模拟部分，我会把结果写出来，重要的是学会蒙特卡洛方法的思路，而不是学会如何解这道题。

首先，由题设知：
$$\begin{gathered} F_Y(y)=\{\begin{array}{ll}0,& y<0\\ y,& 0<y<1\\ 1,& y>1\end{array} \\ P(x=k)=(1-p{)}^{k-1}p \end{gathered}$$
故：

F Z ( z ) = P { Z ≤ z } = P { X Y ≤ z } = P { Y ≤ z } ⋅ P { X = 1 } + P { 2 Y ≤ z } ⋅ P { X = 2 } + P { 3 Y ≤ z } ⋅ P { X = 3 } + ⋯ + P { k Y ≤ z } ⋅ P { X = k } = P { Y ≤ z } ⋅ p + P { Y ≤ z 2 } ⋅ ( 1 − p ) p + ⋯ + P { Y ≤ z k } ⋅ ( 1 − p ) k − 1 p \begin{aligned} F_Z(z) & = P\{Z\le z\}=P\{XY\le z\} \\ & = P\{Y\le z \}\cdot P\{X=1 \}+P\{2Y\le z \}\cdot P\{X=2 \}+P\{3Y\le z \}\cdot P\{X=3 \}+\cdots+P\{kY\le z \}\cdot P\{X=k \} \\ & = P\{Y\le z \}\cdot p+P\{Y\le \frac{z}{2} \}\cdot (1-p)p+\cdots+P\{Y\le \frac{z}{k} \}\cdot (1-p)^{k-1}p \end{aligned} FZ​(z)​=P{Z≤z}=P{XY≤z}=P{Y≤z}⋅P{X=1}+P{2Y≤z}⋅P{X=2}+P{3Y≤z}⋅P{X=3}+⋯+P{kY≤z}⋅P{X=k}=P{Y≤z}⋅p+P{Y≤2z​}⋅(1−p)p+⋯+P{Y≤kz​}⋅(1−p)k−1p​

当 $z<0$ 时，$F_Z(z)=0$

当 $0<z\le 1$ 时，$F_Z(z)=zp+\frac{1}{2}z(1-p)p+\frac{1}{3}z(1-p{)}^{2}p+\cdots +\frac{1}{k}z(1-p{)}^{k-1}p$

当 $1<z\le 2$ 时，$F_Z(z)=p+\frac{1}{2}z(1-p)p+\frac{1}{3}z(1-p{)}^{2}p+\cdots +\frac{1}{k}z(1-p{)}^{k-1}p$

$⋮$

综上所述：
$$F_Z(z)=\{\begin{array}{ll}0,& z<0\\ zp+\frac{1}{2}z(1-p)p+\frac{1}{3}z(1-p{)}^{2}p+\cdots +\frac{1}{k}z(1-p{)}^{k-1}p,& 0<z\le 1\\ p+\frac{1}{2}z(1-p)p+\frac{1}{3}z(1-p{)}^{2}p+\cdots +\frac{1}{k}z(1-p{)}^{k-1}p,& 1<z\le 2\\ ⋮& ⋮\\ p+(1-p)p+(1-p{)}^{2}p+\cdots +\frac{1}{k}z(1-p{)}^{k-1}p,& k-1<z\le k\end{array}$$

2/ 使用蒙特卡洛方法验证

算出来的答案还不知道是否正确，我们可以使用蒙特卡洛方法来验证。其基本思想就是通过生成大量的数据，模拟分布的情况，在数据量足够大的情况，可以较好的把问题模拟出来。

代码在文章的末尾会附上。

首先，根据算出来的答案，可以整理成为：
$$F_z(z)=p\sum _{m=1}^j(1-p{)}^{m-1}+zp\sum _{k=j+1}^{\infty }\frac{1}{k}(1-p{)}^{k-1},j<z\le j+1,j=0,1,2,3,...$$

在代码实现上，不能将 $k$ 一直计算到无穷大，由于当 $k$ 大于一定的数时，对于整个函数的贡献很小，故设定了一个最大的 $k$ 值 $k_{max}=200$。

根据蒙特卡洛方法，我们利用 Python 的 NumPy 库，产生几何分布和在 $[0,1]$ 上的均匀分布随机数，即生成大量的 $X$ 和 $Y$，然后让 $Z=XY$，通过统计，计算在不同的区间上所包含的数据点，画出直方图：

图 1. 利用 MonteCarlo 模拟出的直方图，其中几何分布的 $p=0.1$，所选取的数据点数为 $20000$ 个，每个区间的宽度为 $1$。由于当 $z>40$ 时出现的概率很低，故将区间的最大值设为 $40$。

概率分布函数对 $z$ 求导，得到：
$$f_z(z)=p\sum _{k=j+1}^{\infty }\frac{1}{k}(1-p{)}^{k-1},j<z\le j+1,j=0,1,2,3,...$$
可以知道概率密度函数在不同的区间上有不同的取值，同一区间范围取值相同，即在概率分布函数上看，应该是会随着区间的不同，有不一样的斜率，并且曲线斜率在递减。

图 2. MonteCarlo 直方图（蓝色）和概率密度函数（橙色）对比图。从左到右依次为选择 200、2000、20000 个数据点做出的曲线。

由题目计算出来的函数为概率分布函数，将直方图每个区间的进行 np.cumsum() 函数累加，就可以算出蒙特卡洛模拟的概率分布：

图 3. MonteCarlo （蓝色）和计算出的函数（橙色）对比图。从左到右依次为选择 200、2000、20000 个数据点做出的曲线，蓝色曲线随着数据点选取数量的增加，越接近橙色曲线。

附上模拟的代码：

import numpy as np
import matplotlib.pyplot as plt

np.random.seed(222)

# 把计算得到的函数写成一个函数
def distribution_z(z, p, max_k = 200):
    import math
    j = int(math.floor(z))
    A = 0
    for m in range(1, j + 1):
        A += (1 - p) ** (m - 1)
    A *= p
    
    B = 0
    for k in range(j + 1, max_k + 1):
        a = (1 - p) ** (k - 1)
        a /= k
        B += a
    B *= z * p

    return A + B

def pdf_z(z, p, max_k = 200):
    import math
    j = int(math.floor(z))
    B = 0
    for k in range(j + 1, max_k + 1):
        a = (1 - p) ** (k - 1)
        a /= k
        B += a
    return B * p

p = 0.1
# 选取数据点，点越多越精确
dataPoints = 20000

Unit = np.random.rand(dataPoints)
Geom = np.random.geometric(p, dataPoints)
distri_of_Monte = Geom * Unit

# 概率密度函数 PDF
plt.hist(distri_of_Monte, bins = 40, range = (0, 40))
points_of_z = np.arange(0, 41, 0.01)
pdf_of_z = np.array([pdf_z(zi, p) for zi in points_of_z]) * dataPoints
plt.plot(points_of_z, pdf_of_z)
# print(pdf_of_z)
plt.show()

hist, bin_edges = np.histogram(distri_of_Monte, bins = 40, range = (0, 40))

# 概率分布函数 CDF
hist_list = np.cumsum(hist) / dataPoints

plt.plot(bin_edges[1:], hist_list)

points_of_z = np.arange(1, 41, 0.1)
distri_of_z = [distribution_z(zi, p) for zi in points_of_z]

plt.plot(points_of_z, distri_of_z)

plt.show()