【统计分析】(task4) 应用随机过程（更新ing）

学习总结

• task4学习datawhale萌弟的随机过程。概率论中我们会对随机变量本身进行研究（如随机变量的分布函数、密度函数等等）。不像多维随机变量、大数定律中的无穷多个随机变量是相互独立的，随机过程是也是研究无穷多个随机变量，但它们很多情况下不是相互独立的，注意随机变量还加上了时间维度（不是随机变量，而是普通变量）。
• 白噪声随机过程、滑动平均随机过程都是宽平稳过程，而宽平稳过程的一个显然特性就是均值函数与方差函数随时间不变。宽平稳过程在实际中容易操作，因此在时间序列中运用十分广泛。
• 描述随机过程的武器——概率空间、随机变量等；能否研究一个随机过程的关键就是减少问题的维度-这也是物理的核心思想。而达到这个目的的神器有：
  • 马尔科夫过程(Markov Processes)—随机过程的精华，即随机过程的每一步的结果最多只与上一步有关，而与其它无关。
  • 马尔可夫链(Markov chain)—即用数学语言表述makov过程；如果一个过程是markov过程，这个过程就得到了神简化，你只需要知道第n步是如何与第n-1步相关的，一般由一组条件概率表述，就可以求得整个过程。
• 随机过程是在随机变量的基础上添加一个时间维度，因此随机过程的数字特征也由随机变量的数字变成了以时间$t$为变量的函数。随机过程 { X ( t ) , t ∈ T } \{X(t), t \in T\} {X(t),t∈T}的常见的数字特征有：均值函数、方差函数、协方差函数、自相关函数等等。
• 宽平稳和独立增量过程：
  • 宽平稳：白噪声随机过程、滑动平均随机过程都是宽平稳过程，而宽平稳过程的一个显然特性就是均值函数与方差函数随时间不变。宽平稳过程在实际中容易操作，因此在时间序列中运用十分广泛。
  • 随机过程任意两个时间的随机变量$X(t)$与$X(s)$之间一般是不独立的。但是，大量的事实证明许多随机过程的增量是相互独立的，即独立增量过程。

零、导言

• 概率论：仅仅研究的是在某一个时点，某个随机变量（随机向量）的分布情况，如：2022年07月03日这天的某个股票价格在开盘前是一个随机变量，因为在开盘前这个价格可以是$(0,+\infty )$的任意值，但是每个价格的概率是有差别的。
• 随机过程：描述的是一段时间内随机变量的分布的变化，如：连续的两天的某只股票在开盘前的价格，这两个随机变量的分布按照时间变化的性质可以使用随机过程来描述。
• 总结：随机过程就是在随机变量的基础上加入了时间维度（时间维度不是随机变量，只是普通变量），随机过程与概率论有很强的关联性。

一、随机过程的基本概念

在概率论中，随机变量的个数主要有一维随机变量、多维随机变量、大数定律中也涉及无穷个随机变量，但它们之间是相互独立的。我们不仅需要对随机变量本身进行研究（如随机变量的分布函数、密度函数等等），也要开始研究随机现象的变化过程。

在这个前提下，我们必须考虑无限多个随机变量的一次观测，与大数定律不同的是：这些无限多个随机变量可能并不是相互独立的。因此，我们将需要研究的无穷多个随机变量（很多情况不是相互独立的）称为随机过程。

1.1 随机过程的几个基本概念

（1）概率空间$(\Omega ,\mathcal{F},P)$

样本点or基本事件$\omega$：随机试验的可能结果称为样本点或基本事件。
样本空间$\Omega$：样本点的全体称为样本空间。
样本空间 $\Omega$ 称为必然事件， 空集 $\varnothing$ 称为不可能事件。
事件：$\Omega$ 的子集 $A$ 由基本事件组成， 称为事件。

在实际研究问题中，往往并不是对样本空间的所有子集（所有事件）都关心，而是关心某些重要事件以及重要事件的可能性大小，因此需要一个概念来刻画我们感兴趣的事件：

概率空间$(\Omega ,\mathcal{F},P)$：设 $\Omega$ 是一个样本空间（或任意一个集合)， $\mathcal{F}$ 是 $\Omega$ 的某些子集组成的集合族（可以理解为集合的集合）。如果满足

• (1) $\Omega \in \mathcal{F}$；
• (2) 若 $A\in \mathcal{F}$， 则 $A^{\mathrm{c}}=\Omega \backslash A\in \mathcal{F}$;
• (3) 若 $A_n\in \mathcal{F}(n=1,2,\cdots )$， 则 ${\bigcup }_{n=1}{A}_n\in \mathcal{F}$。

则称 $\mathcal{F}$ 为 $\Omega$ 上的一个 $\sigma$ 代数， $(\Omega ,\mathcal{F})$ 称为可测空间，$(\Omega ,\mathcal{F},P)$为概率空间。

（2）随机过程：

随机过程是概率空间$(\Omega ,\mathcal{F},P)$上的一族随机变量 { X ( t ) , t ∈ T } \{X(t), t \in T \} {X(t),t∈T}，其中$T$被称为指标集或参数集。随机过程的概念，可以理解为一组按时间不断排列的随机变量，这组随机变量中的随机变量间可能有关系，可能没有关系。

案例一：随机漫步的傻瓜

一个醉汉傻瓜在路上行走，这个醉汉每次只会向前或者向后走一步，他会以概率$p$前进一步，以概率1-p后退一步。现在，我们需要模拟醉汉的行走轨迹。（假设$X(t)$表示醉汉在$t$时刻的位置，则$X(t)$就是一个随机过程，这个随机过程叫做随机游动。）

模拟醉汉傻瓜的前进轨迹（一条轨迹）：

# 模拟醉汉傻瓜的前进轨迹（一条轨迹）
from scipy.stats import bernoulli
def go_random(p, n_steps):
    random_step = [] #记录前进的位置
    random_step.append(0)  # 初始位置在0点
    r = np.array(bernoulli.rvs(p, size=n_steps))
    r[r==0] = -1
    random_step += np.cumsum(r).tolist()
    print(random_step) # 打印当前位置
    return random_step

random_steps = go_random(0.6, 10)
plt.figure(figsize=(8, 6))
plt.scatter(x=np.arange(len(random_steps)), y=random_steps)
plt.plot(np.arange(len(random_steps)), random_steps)
plt.xlabel("t")
plt.ylabel("X(t)")
plt.title("随机游动1次观测")
plt.xticks(np.arange(11))
plt.yticks(np.arange(-10, 11))
plt.show()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21

模拟醉汉傻瓜的前进轨迹（10条轨迹）：

# 模拟醉汉傻瓜的前进轨迹（10条轨迹）
import matplotlib.pyplot as plt
plt.rcParams["font.family"]="Songti SC"

from scipy.stats import bernoulli
def go_random(p, n_steps):
    random_step = [] #记录前进的位置
    random_step.append(0)  # 初始位置在0点
    r = np.array(bernoulli.rvs(p, size=n_steps))
    r[r==0] = -1
    random_step += np.cumsum(r).tolist()
    return random_step

plt.figure(figsize=(8, 6))
for i in range(10):
    random_steps = go_random(0.5, 10)
    plt.scatter(x=np.arange(len(random_steps)), y=random_steps)
    plt.plot(np.arange(len(random_steps)), random_steps, label='第'+str(i+1)+'次模拟')
plt.xlabel("t")
plt.ylabel("X(t)")
plt.title("随机游动10次观测")
plt.xticks(np.arange(11))
plt.yticks(np.arange(-10, 11))
plt.legend()
plt.show()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25

上图：傻瓜在每个时刻的位置都是一个随机变量，如在$t=2$处，有三个观测点，这三个观测点分别都是来自于$t=2$时的随机变量$X(2)$的一个采样（$X(2)$是一个随机变量而不是一个数值或者一个普通变量）。

（3）有限维分布族：

对于一个或者有限多个随机变量来说，了解随机变量的分布函数就能掌握随机变量的全部面貌。但是对于随机过程 { X ( t ) , t ∈ T } \{X(t), t \in T\} {X(t),t∈T}来说，由于随机变量的个数是无限的，因此想要了解随机过程的全貌，必须知道：

• 每个时刻$t$的随机变量$t\in T,X(t)$的分布函数 F ( t , x ) = P { X ( t ) ⩽ x } F(t, x)=P\{X(t) \leqslant x\} F(t,x)=P{X(t)⩽x}
• 任意两个时刻的$t_1,t_2\in T$的二维随机变量$(X(t_1),X(t_2))$的二维联合分布函数 F t 1 , t 2 ( x 1 , x 2 ) = P { X ( t 1 ) ⩽ x 1 , X ( t 2 ) ⩽ x 2 } F_{t_1, t_2}(x_{1}, x_{2})=P\left\{X\left(t_{1}\right) \leqslant x_{1}, X\left(t_{2}\right) \leqslant x_{2}\right\} Ft1​,t2​​(x1​,x2​)=P{X(t1​)⩽x1​,X(t2​)⩽x2​}
• 任意$n$个时刻的随机变量的联合分布函数 F t 1 , t 2 , ⋯   , t n ( x 1 , x 2 , ⋯   , x n ) = P { X ( t 1 ) ⩽ x 1 , X ( t 2 ) ⩽ x 2 , ⋯   , X ( t n ) ⩽ x n } F_{t_{1}, t_{2}, \cdots, t_{n}}\left(x_{1}, x_{2}, \cdots, x_{n}\right)=P\left\{X\left(t_{1}\right) \leqslant x_{1}, X\left(t_{2}\right) \leqslant x_{2}, \cdots, X\left(t_{n}\right) \leqslant x_{n}\right\} Ft1​,t2​,⋯,tn​​(x1​,x2​,⋯,xn​)=P{X(t1​)⩽x1​,X(t2​)⩽x2​,⋯,X(tn​)⩽xn​}

因此，描述完整的随机过程的全貌需要一族分布函数，这族分布函数叫做随机过程 { X ( t ) , t ∈ T } \{X(t), t \in T\} {X(t),t∈T}的有限维分布族:
{ F t 1 , t 2 , ⋯   , t n ( x 1 , x 2 , ⋯   , x n ) , t 1 , t 2 , ⋯   , t n ∈ T , n ⩾ 1 } \left\{F_{t_{1}, t_{2}, \cdots, t_{n}}\left(x_{1}, x_{2}, \cdots, x_{n}\right), t_{1}, t_{2}, \cdots, t_{n} \in T\right. ,n \geqslant 1\} {Ft1​,t2​,⋯,tn​​(x1​,x2​,⋯,xn​),t1​,t2​,⋯,tn​∈T,n⩾1}

（4）随机过程的数字特征：

• 在概率论中，有时候不需要知道随机变量的方方面面，因此有时候我们不需要存储一整个随机变量的分布函数，因为这是件十分困难的事情。我们可能需要知道的是这个随机变量的某个特征，如：平均水平、中间位置水平、离散程度等等，这些数字对应的随机变量数字特征：数学期望、方差/标准差、分位数等等。
• 在随机过程中，随机过程 { X ( t ) , t ∈ T } \{X(t), t \in T\} {X(t),t∈T}的有限维分布族能反映随机过程的方方面面，但是有时候求随机过程的有限维分布族是一件十分困难的事情，因此我们也想像随机变量一样使用“数字特征”的概念。
  • 随机过程 { X ( t ) , t ∈ T } \{X(t), t \in T\} {X(t),t∈T}的常见的数字特征有：均值函数、方差函数、协方差函数、自相关函数等等。
  • 随机变量的数字特征是一个数字，而随机过程的数字特征是一个函数。因为随机过程是在随机变量的基础上添加一个时间维度，因此随机过程的数字特征也由随机变量的数字变成了以时间$t$为变量的函数。

（4.1）均值函数：

设 { X ( t ) , t ∈ T } \{X(t), t \in T\} {X(t),t∈T} 是一随机过程，称 $X(t)$ 的期望 ${\mu }_X(t)=E[X(t)]$ 为过程的均值函数（如果存在的话)。

（1）模拟醉汉傻瓜的前进轨迹（1000条轨迹），并绘制均值函数

# 模拟醉汉傻瓜的前进轨迹（1000条轨迹），并绘制均值函数（p=0.5）
from scipy.stats import bernoulli
def go_random(p, n_steps):
    random_step = [] #记录前进的位置
    random_step.append(0)  # 初始位置在0点
    r = np.array(bernoulli.rvs(p, size=n_steps))
    r[r==0] = -1
    random_step += np.cumsum(r).tolist()
    return random_step

mean_step = []
plt.figure(figsize=(8, 6))
for i in range(1000):
    random_steps = go_random(0.5, 10)
    mean_step.append(random_steps) #存储每次实验的结果
    plt.scatter(x=np.arange(len(random_steps)), y=random_steps)
    plt.plot(np.arange(len(random_steps)), random_steps)
mean_var = np.array(mean_step).mean(axis=0)
mean_var = np.mean(np.array(mean_step), axis=0)
plt.plot(np.arange(len(random_steps)), mean_var, lw=3, c='red', ls='--', label=r'均值函数$p=0.5$')
plt.xlabel("t")
plt.ylabel("X(t)")
plt.title("随机游动1000次观测")
plt.xticks(np.arange(11))
plt.yticks(np.arange(-10, 11))
plt.legend()
plt.show()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27

（2）模拟醉汉傻瓜的前进轨迹（1000条轨迹），并绘制均值函数（p=0.2）

# 模拟醉汉傻瓜的前进轨迹（1000条轨迹），并绘制均值函数（p=0.2）
from scipy.stats import bernoulli
def go_random(p, n_steps):
    random_step = [] #记录前进的位置
    random_step.append(0)  # 初始位置在0点
    r = np.array(bernoulli.rvs(p, size=n_steps))
    r[r==0] = -1
    random_step += np.cumsum(r).tolist()
    return random_step

mean_step = []
plt.figure(figsize=(8, 6))
for i in range(1000):
    random_steps = go_random(0.2, 10)
    mean_step.append(random_steps) #存储每次实验的结果
    plt.scatter(x=np.arange(len(random_steps)), y=random_steps)
    plt.plot(np.arange(len(random_steps)), random_steps)
mean_var = np.array(mean_step).mean(axis=0)
mean_var = np.mean(np.array(mean_step), axis=0)
plt.plot(np.arange(len(random_steps)), mean_var, lw=3, c='red', ls='--', label=r'均值函数$p=0.2$')
plt.xlabel("t")
plt.ylabel("X(t)")
plt.title("随机游动1000次观测")
plt.xticks(np.arange(11))
plt.yticks(np.arange(-10, 11))
plt.legend()
plt.show()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27

（4.2）方差函数

与均值函数对应的，方差函数$\mathrm{Var}[X(t)]=\gamma (t,t)$通俗来说就是每个时刻的随机变量的方差形成的关于时间$t$的函数。

案例：$X(t)=X_0+tV(a⩽t⩽b)$， 其中 $X_0$ 和 $V$ 是相互独立且均服从 $N(0,1)$ 分布的随机变量，求随机过程 { X ( t ) , t ∈ T } \{X(t), t \in T \} {X(t),t∈T}的方差函数。

均值函数：
μ X ( t ) = E [ X ( t ) ] = E ( X 0 + t V ) = E ( X 0 ) + t E ( V ) = 0 μX(t)=E[X(t)]=E(X0+tV)=E(X0)+tE(V)=0

μX​(t)​=E[X(t)]=E(X0​+tV)=E(X0​)+tE(V)=0​
方差函数：
$$\begin{array}{rl}\mathrm{Var}[X(t)]& =\gamma (t,t)=E[X^2(t)]-E[X(t){]}^2\\ & =E[X_0^2+t^2{V}^2+2t{X}_{0}V]\\ & =1+t^2\end{array}$$
【栗子】模拟随机过程{X(t) t in T} 1000次的观测轨迹，并绘制均值函数和方差函数

# 模拟随机过程{X(t) t in T} 1000次的观测轨迹，并绘制均值函数和方差函数
def x_t_random(a, b, n_step):
    sample_list = []
    t_T = np.linspace(a, b, n_step)
    X0 = np.random.randn()
    for t in t_T:
        sample_list.append(X0 + t * np.random.randn())
    return sample_list
a, b, n_step, n_times = 0, 1, 100, 1000  # n_times代表模拟1000次
n_times_sample_list = [x_t_random(a, b, n_step) for _ in range(n_times)]  #所有实验结果
plt.figure(figsize=(8, 6))
for i in range(n_times):
    plt.plot(np.linspace(a, b, n_step), n_times_sample_list[i])
mean_var = np.array(n_times_sample_list).mean(axis=0)
var_var = np.array(n_times_sample_list).var(axis=0)
plt.plot(np.linspace(a, b, n_step), mean_var, lw=4, c='red', ls='-', label=r'均值函数$\mu_{X}(t)=0$')
plt.plot(np.linspace(a, b, n_step), var_var, lw=4, c='blue', ls='-', label=r'方差函数$\operatorname{Var}[X(t)]=1+t^2$')
plt.xlabel("t")
plt.ylabel("X(t)")
plt.title("随机过程1000次观测")
plt.xticks(np.linspace(a, b, 10))
plt.xlim(0,1)
plt.legend()
plt.show()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24

（4.3）协方差函数：

γ ( t 1 , t 2 ) = E { [ X ( t 1 ) − μ X ( t 1 ) ] [ X ( t 2 ) − μ X ( t 2 ) ] } ( t 1 , t 2 ∈ T ) \gamma\left(t_{1}, t_{2}\right)=E\left\{\left[X\left(t_{1}\right)-\mu_{X}\left(t_{1}\right)\right]\left[X\left(t_{2}\right)-\mu_{X}\left(t_{2}\right)\right]\right\}\left(t_{1}, t_{2} \in T\right) γ(t1​,t2​)=E{[X(t1​)−μX​(t1​)][X(t2​)−μX​(t2​)]}(t1​,t2​∈T)

在概率论中，讨论两个随机比变量的相关程度可以使用协方差，而$n$维随机变量则可以使用协方差矩阵。因此，在随机过程中，随机变量增加了一个时间维度后，协方差函数就变成了二元函数。

案例：$X(t)=X_0+tV(a⩽t⩽b)$， 其中 $X_0$ 和 $V$ 是相互独立且均服从 $N(0,1)$ 分布的随机变量，求随机过程 { X ( t ) , t ∈ T } \{X(t), t \in T \} {X(t),t∈T}的协方差函数。
γ ( t 1 , t 2 ) = E [ ( X ( t 1 ) − 0 ) ( X ( t 2 ) − 0 ) ] = E [ ( X 0 + t 1 V ) ( X 0 + t 2 V ) ] = E ( X 0 2 ) + t 1 t 2 E ( V 2 ) = 1 + t 1 t 2 γ(t1,t2)=E[(X(t1)−0)(X(t2)−0)]=E[(X0+t1V)(X0+t2V)]=E(X20)+t1t2E(V2)=1+t1t2

γ(t1​,t2​)​=E[(X(t1​)−0)(X(t2​)−0)]=E[(X0​+t1​V)(X0​+t2​V)]=E(X02​)+t1​t2​E(V2)=1+t1​t2​​

（4.4）自相关函数：

$$R_X(s,t)=E[X(s)X(t)]$$ $$(s,t\in T)$$

# 模拟随机过程{X(t) t in T} 1000次的观测轨迹，并绘制协方差函数
def x_t_random(a, b, n_step):
    sample_list = []
    t_T = np.linspace(a, b, n_step)
    X0 = np.random.randn()
    for t in t_T:
        sample_list.append(X0 + t * np.random.randn())
    return sample_list
a, b, n_step, n_times = 0, 1, 30, 1000  # n_times代表模拟1000次
n_times_sample_arr = np.array([x_t_random(a, b, n_step) for _ in range(n_times)]) #所有实验结果
fig = plt.figure(figsize=(8,6))
ax = fig.add_subplot(projection='3d')
t1 = np.linspace(a, b, n_step)
t2 = np.linspace(a, b, n_step)
t1, t2 = np.meshgrid(t1, t2)
gamma_t = []
for i in range(n_step):
    tmp = []
    for j in range(n_step):
        tmp.append(np.mean((n_times_sample_arr[:,i]-np.mean(n_times_sample_arr[:,i]))*(n_times_sample_arr[:, j]-np.mean(n_times_sample_arr[:,j]))))
    gamma_t.append(tmp)
gamma_t = np.array(gamma_t)
print(gamma_t.shape)
ax.plot_surface(t1, t2, gamma_t, cmap='plasma')
ax.set_xlabel("t1")
ax.set_ylabel("t2")
ax.set_zlabel(r"$\gamma(t_1, t_2)$")
ax.set_title(r"$X(t)=X_{0}+t V(a \leqslant t \leqslant b)$协方差函数")
plt.show()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29

# 模拟随机过程{X(t) t in T} 1000次的观测轨迹，并绘制自相关函数
def x_t_random(a, b, n_step):
    sample_list = []
    t_T = np.linspace(a, b, n_step)
    X0 = np.random.randn()
    for t in t_T:
        sample_list.append(X0 + t * np.random.randn())
    return sample_list
a, b, n_step, n_times = 0, 1, 30, 1000  # n_times代表模拟1000次
n_times_sample_arr = np.array([x_t_random(a, b, n_step) for _ in range(n_times)]) #所有实验结果
fig = plt.figure(figsize=(8,6))
ax = fig.add_subplot(projection='3d')
t1 = np.linspace(a, b, n_step)
t2 = np.linspace(a, b, n_step)
t1, t2 = np.meshgrid(t1, t2)
R_t = []
for i in range(n_step):
    tmp = []
    for j in range(n_step):
        tmp.append(np.mean((n_times_sample_arr[:,i])*(n_times_sample_arr[:, j])))
    R_t.append(tmp)
R_t = np.array(R_t)
print(R_t.shape)
ax.plot_surface(t1, t2, R_t, cmap='plasma')
ax.set_xlabel("t1")
ax.set_ylabel("t2")
ax.set_zlabel(r"$\gamma(t_1, t_2)$")
ax.set_title(r"$X(t)=X_{0}+t V(a \leqslant t \leqslant b)$自相关函数")
plt.show()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29

1.2 随机过程的基本类型

对于两个随机变量来说，两个随机变量的关系可以分为：一个随机变量与另一个随机变量之间是否同分布、一个随机变量与另一个随机变量之间是否独立。

对于随机过程来说，由于在任意$n$个时刻可以形成一个$n$维随机变量，这个$n$维随机变量由它的联合分布确定，而取不同的时间所得到的$n$维随机变量也是不一样的，因此可以根据随机过程内的$n$维随机变量之间的关系，即：是否同分布和是否独立，可以将随机过程分为平稳过程与独立增量过程等等。

（1）严平稳过程：

• 平稳：有一种情况是理论上保证“理论上的不会变化太大”，另一种情况是现实中看起来不会变化太大。
• 严平稳：描述的是随机变量在经过同一个时间段后的随机变量之间的分布保持不变，由于随机变量的分布可以完全确定随机变量，而随机变量的分布保持一致就说明了随机变量在“理论上不会变化太大”，具体来说：

【定义】如果随机过程 { X ( t ) , t ∈ T } \{X(t), t \in T\} {X(t),t∈T} 对任意的 $t_1,t_2,\cdots ,t_n\in T$ 和任意的 $h$ (使得 $t_i+h\in T,i=1,2,\cdots ,n)$ 有 $$\left(X\left(t_1+h\right),X\left(t_2+h\right),\cdots ,X\left(t_n+h\right)\right)$$ 与 $\left(X\left(t_1\right),X\left(t_2\right),\cdots ,X\left(t_n\right)\right)$ 具有相同的联合分布， 记为：
$$\left(X\left(t_1+h\right),X\left(t_2+h\right),\cdots ,X\left(t_n+h\right)\right)\stackrel{\mathrm{d}}{=}\left(X\left(t_1\right),X\left(t_2\right),\cdots ,X\left(t_n\right)\right)$$
则称 { X ( t ) , t ∈ T } \{X(t), t \in T\} {X(t),t∈T} 为严平稳的。

对于严平稳过程而言, 有限维分布关于时间是平移不变的，这个要求在绝大多数情况都过于严格，而且不容易被验证，毕竟计算一个有限维分布本身就是一件十分困难的事。通常情况下，我们需要实操性更强的平稳性，那就是所谓的“看起来平稳”的弱平稳。

（2）宽平稳过程：

宽平稳过程通俗来说就是看起来平稳。在现实世界中，我们主要通过研究随机变量的一般水平（均值）和分散程度来大致刻画随机变量，并不需要完整知道随机变量的分布。因此，所谓的看起来平稳指的就是：随机过程的均值函数随时间平移以及随机过程的方差函数也是随着时间是无关的。

【定义】如果随机过程 $X(t)$ 的所有二阶矩都存在， 并且 $E[X(t)]=\mu$, 协 方差函数 $\gamma (t,s)$ 只与时间差 $t-s$ 有关， 则称 { X ( t ) , t ∈ T } \{X(t), t \in T\} {X(t),t∈T} 为宽平稳过程或二阶平稳过程。

案例一：白噪声序列（随机过程）

设 { X n , n = 0 , 1 , ⋯   } \left\{X_{n}, n=0,1, \cdots\right\} {Xn​,n=0,1,⋯} 为一列两两互不相关的随机变量序列， 满足 $E\left(X_n\right)=0(n=0,1,\cdots )$， 且
E ( X m X n ) = { 0 , m ≠ n σ 2 , m = n E\left(X_{m} X_{n}\right)= {0,m≠nσ2,m=n

E(Xm​Xn​)={0,σ2,​m=nm=n​
则随机过程 { X n , n = 0 , 1 , ⋯   } \left\{X_{n}, n=0,1, \cdots\right\} {Xn​,n=0,1,⋯}为白噪声序列。

为什么宽平稳过程需要“协方差函数 $\gamma (t,s)$ 只与时间差 $t-s$ 有关”这个条件。如果一个随机过程的协方差函数与所在时间有关，那么意味着你这个过程的性质随时都在发生变化；与之对应的，如果一个过程的协方差函数只与时间间隔有关而与所在时刻无关，那么这意味着这个过程在所有时间起点上的性质都是相同的，这就看起来很平稳。

白噪声序列的平稳性是如何体现的：

• 白噪声过程的一个样本轨迹是什么样子的
• 白噪声过程的均值函数和方差函数

白噪声过程：

# 白噪声过程
def get_white_noise(n_steps):
    white_noise_ls = []
    for step in range(n_steps):
        noise_sample = np.random.randn()
        white_noise_ls.append(noise_sample)
    return white_noise_ls
n_steps = 1000
t = np.arange(n_steps)
plt.figure(figsize=(8,6))
plt.plot(t, np.array(get_white_noise(n_steps=n_steps)), lw=3, color='blue', alpha=0.75)
plt.xlabel("t")
plt.ylabel("X(t)")
plt.title("白噪声过程")
plt.show()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

白噪声过程的均值和方差函数：

# 白噪声过程的均值和方差函数
def get_white_noise(n_steps):
    white_noise_ls = []
    for step in range(n_steps):
        noise_sample = np.random.randn()
        white_noise_ls.append(noise_sample)
    return white_noise_ls

n_steps, n_times = 10, 1000  # n_steps随机过程的时间t, n_times模拟这个过程的次数
simulate_result = np.array([get_white_noise(n_steps=n_steps) for _ in range(n_times)])
print(simulate_result.shape)
t_ls = np.arange(n_steps)
plt.figure(figsize=(8,6))
for i in range(n_times):
    plt.plot(t_ls, simulate_result[i])  # 绘制每次的轨迹
plt.plot(t_ls, np.mean(simulate_result, axis=0), lw=4, color='red', label='mean_func')
plt.plot(t_ls, np.var(simulate_result, axis=0), lw=4, color='blue', label='var_func')
plt.xlabel("t")
plt.ylabel("X(t)")
plt.title("白噪声过程及其均值、方差函数")
plt.xticks(np.arange(n_steps))
plt.legend()
plt.show()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23

案例二：滑动平均过程

设 { ε n , n = 0 , ± 1 , ⋯   } \left\{\varepsilon_{n}, n=0, \pm 1, \cdots\right\} {εn​,n=0,±1,⋯} 为一列两两互不相关的有相同均值 $\mu$ 和相同 方差 ${\sigma }^2$ 的随机变量序列， $a_1,a_2,\cdots ,a_k$ 为任意 $k$ 个实数。考虑下列定义的序列：
$$X_n=a_1{\epsilon }_n+a_2{\epsilon }_{n-1}+\cdots +a_k{\epsilon }_{n-k+1},n=0,\pm 1,\pm 2,\cdots$$
那么随机过程 { X ( t ) } \{X(t) \} {X(t)}被称为滑动平均序列。

滑动平均序列（随机过程）的平稳性：

• 一条滑动平均序列的样本轨迹
• 滑动平均序列的均值和方差函数

# 获得一条滑动平均序列
def get_esp_ls(n_steps):
    esp_ls = []
    for i in range(n_steps):
        esp_ls.append(np.random.randn())
    return esp_ls

def get_moving_avg(esp_ls, ak_ls):
    xn = []
    for i in range(len(esp_ls)-len(ak_ls)+1):
        tmp = 0
        for j in range(len(ak_ls)):
            tmp += esp_ls[i+j]*ak_ls[j]
        xn.append(tmp)
    return xn

n_steps = 100
esp_ls = get_esp_ls(n_steps=n_steps)
xn = get_moving_avg(esp_ls=esp_ls, ak_ls=[0.2, 0.6, 0.2])
t = np.arange(n_steps)
plt.figure(figsize=(8,6))
plt.plot(t, esp_ls, lw=2, color='blue', alpha=0.75, label=r'$\{ \epsilon_t \}$')
plt.plot(t[2:], xn, color='red', alpha=0.75, label=r'$\{ X(t) \}$')
plt.legend()
plt.show()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25

# 滑动平均序列(随机过程)的均值函数和方差函数
def get_esp_ls(n_steps):
    esp_ls = []
    for i in range(n_steps):
        esp_ls.append(np.random.randn())
    return esp_ls

def get_moving_avg(esp_ls, ak_ls):
    xn = []
    for i in range(len(esp_ls)-len(ak_ls)+1):
        tmp = 0
        for j in range(len(ak_ls)):
            tmp += esp_ls[i+j]*ak_ls[j]
        xn.append(tmp)
    return xn

n_steps, n_times = 100, 1000
simulate_result = np.array([get_moving_avg(get_esp_ls(n_steps=n_steps), [0.2, 0.6, 0.2]) for _ in range(n_times)])
t = np.arange(n_steps)
plt.figure(figsize=(8,6))
for i in range(n_times):
    plt.plot(t[2:], simulate_result[i])
plt.plot(t[2:], np.mean(simulate_result, axis=0), lw=4, color='red', label='均值函数')
plt.plot(t[2:], np.var(simulate_result, axis=0), lw=4, color='blue', label='方差函数')
plt.legend()
plt.show()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26

总结上面两个栗子：白噪声随机过程、滑动平均随机过程都是宽平稳过程，而宽平稳过程的一个显然特性就是均值函数与方差函数随时间不变。宽平稳过程在实际中容易操作，因此在时间序列中运用十分广泛。

（3）独立增量过程：

对于两个随机变量来说，两个随机变量的关系可以分为：

• 一个随机变量与另一个随机变量之间是否同分布（同分布，即平稳性）
• 一个随机变量与另一个随机变量之间是否独立。

现在来看另一个关系：一个随机变量与另一个随机变量之间是否独立。一般来说，随机过程任意两个时间的随机变量$X(t)$与$X(s)$之间是不独立的，因为在很多时候，前一个时刻的数据会深刻影响下一个时刻的数据，如：昨天股票跌停了，今天的股票价格必定会受到昨日信息的影响。

因此，随机过程任意两个时间的随机变量$X(t)$与$X(s)$之间一般是不独立的。但是，大量的事实证明许多随机过程的增量是相互独立的，这种随机过程被称为独立增量过程：

【定义】如果对任意 $t_1,t_2,\cdots ,t_n\in T,t_1<t_2<\cdots <t_n$， 随机变量 $X\left(t_2\right)-$ $X\left(t_1\right),X\left(t_3\right)-X\left(t_2\right),\cdots ,X\left(t_n\right)-X\left(t_{n-1}\right)$ 是相互独立的， 则称 { X ( t ) , t ∈ T } \{X(t), t \in T\} {X(t),t∈T} 是 独立增量过程。

如果对任意 $t_1,t_2$, 有 $$X\left(t_1+h\right)-X\left(t_1\right)\stackrel{\mathrm{d}}{=}X\left(t_2+h\right)-X\left(t_2\right)$$, 则称 { X ( t ) , t ∈ T } \{X(t), t \in T\} {X(t),t∈T} 是 平稳增量过程。

兼有独立增量和平稳增量的过程称为平稳独立增量过程。

【布朗运动】我们来看一个平稳独立增量过程的随机过程案例：布朗运动，先来通俗定义布朗运动。考虑在一直线上的简单的、 对称的随机游动。 设质点每经过 $\Delta t$ 时间，随机地以概率 $p=1/2$ 向右移 $\Delta x>0$，以概率 $q=1/2$ 向左移动一个 $\Delta x$ ，且每次移动相互独立。 若 $X_t$ 表示 $t$ 时刻 质点的位置，且有 $\Delta t\to 0$ 时 $\Delta x=c\sqrt{\Delta t}$ ，则
$$X_t\sim N\left(0,c^{2}t\right)$$

这个随机过程就是著名的“布朗运动”，将布朗运动与股票价格行为联系在一起，进而建立起维纳过程的数学模型是本世纪的一项具有重要意义的金融创新，在现代金融数学中占有重要地位。迄今，普遍的观点仍认为，股票市场是随机波动的，随机波动是股票市场最根本的特性，是股票市场的常态。布朗运动假设是现代资本市场理论的核心假设。（《随机模拟》章节对布朗运动进行详细的阐述）

我们来看看平稳独立增量的随机过程有什么有用的性质，我们先用python模拟一段随机漫步，观察随机漫步的均值函数和增量的均值函数：

# 模拟随机游走并绘制均值函数
def get_brown(p, n_steps):
    brown_ls = []
    x_0 = 0
    brown_ls.append(x_0)
    for i in range(n_steps-1):
        p_i = np.random.rand()
        if p_i <= p:
            brown_ls.append(x_0+1)
            x_0 += 1
        else:
            brown_ls.append(x_0-1)
            x_0 -= 1
    return brown_ls
simulate_brown_ls = []
p, n_steps, n_times = 0.8, 20, 10000
for i in range(n_times):
    simulate_brown_ls.append(get_brown(p = p, n_steps=n_steps))
simulate_brown_arr = np.array(simulate_brown_ls)
mean_func = np.mean(simulate_brown_arr, axis=0)
t = np.arange(n_steps)
plt.figure(figsize=(8,6))
for i in range(n_times):
    plt.plot(t, simulate_brown_arr[i, :])
plt.plot(t, mean_func, lw=4, color='red', alpha=0.75, label='均值函数')
plt.xlabel("t")
plt.ylabel("X(t)")
plt.title("随机游走及其均值函数")
plt.xticks(np.arange(0, n_steps, 2))
plt.yticks(np.arange(-20, 20, 2))
plt.legend()
plt.show()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32

# 模拟随机游走的增量过程并绘制均值函数
def get_del_brown(p, n_steps):
    brown_ls = []
    x_0 = 0
    brown_ls.append(x_0)
    for i in range(n_steps-1):
        p_i = np.random.rand()
        if p_i <= p:
            brown_ls.append(x_0+1)
            x_0 += 1
        else:
            brown_ls.append(x_0-1)
            x_0 -= 1
    delta_brown_ls = []
    for i, xi in enumerate(brown_ls):
        if i != 0:
            delta_brown_ls.append(xi - brown_ls[i-1])
    return delta_brown_ls
simulate_del_brown_ls = []
p, n_steps, n_times = 0.8, 20, 10000
for i in range(n_times):
    simulate_del_brown_ls.append(get_del_brown(p = p, n_steps=n_steps))
simulate_del_brown_ls = np.array(simulate_del_brown_ls)
mean_func = np.mean(simulate_del_brown_ls, axis=0)
t = np.arange(1, n_steps)
plt.figure(figsize=(8,6))
for i in range(n_times):
    plt.plot(t, simulate_del_brown_ls[i, :])
plt.plot(t, mean_func, lw=4, color='red', alpha=0.75, label='均值函数')
plt.xlabel("t")
plt.ylabel("X(t)")
plt.title("随机游走的增量过程及其均值函数")
plt.xticks(np.arange(0, n_steps, 2))
plt.yticks(np.arange(-2, 2, 1))
plt.legend()
plt.show()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36

总结上面2个栗子：如果一个随机过程是平稳独立增量过程，即：增量过程平稳且独立（均值函数是一个常数），那么这个随机过程的均值函数是一个线性函数。

以上的平稳过程和独立增量过程等等的讨论都是基于随机过程的一些很宏观的性质进行讨论的，在实际的运用中，我们需要对实际应用场景进行讨论和建模，因此会衍生出类似于：泊松过程、马尔可夫过程等等的随机过程，这些随机过程往往是基于应用场景进行讨论和定义的。

二、泊松过程

2.1 什么是计数过程

泊松过程是一种十分常见的计数过程，只要我们想要研究所观察的事件出现的次数, 都可以使用计数过程来描述，如：某个时间段内正在银行排队等待服务的顾客书、某条道路上的车辆数量、某个时间段内某地区的出生人数或死亡人数等等。那怎么样准确定义计数过程呢？

随机过程 { N ( t ) , t ⩾ 0 } \{N(t), t \geqslant 0\} {N(t),t⩾0} 称为计数过程， 如果 $N(t)$ 表示从 0 到 $t$ 时 刻某一特定事件 $A$ 发生的次数， 它具备以下两个特点:

• （1） $N(t)⩾0$ 且取值为整数;
• （2）当 $s<t$ 时, $N(s)⩽N(t)$ 且 $N(t)-N(s)$ 表示 $(s,t]$ 时间内事件 $A$ 发生的 次数。

2.2 泊松过程的几个定义

由于计数过程的概念太大，如：某个事件出现的次数服从什么分布、某两个事件出现需要等待的时间满足什么分布等等，直接使用计数过程描述事物出现的次数会十分困难。泊松过程是一个定义更精确的计数过程，具体来说：

（1）泊松分布的第一个定义

【定义一】计数过程 { N ( t ) , t ⩾ 0 } \{N(t), t \geqslant 0\} {N(t),t⩾0} 称为参数为 $\lambda (\lambda >0)$ 的 Poisson 过程， 如果

• (1) $N(0)=0$；
• (2) 过程有独立增量；
• (3) 在任一长度为 $t$ 的时间区间中事件发生的次数服从均值为 $\lambda t$ 的 Poisson 分布, 即对一切 $s⩾0,t>0$, 有
  P { N ( t + s ) − N ( s ) = n } = e − λ t ( λ t ) n n ! , n = 0 , 1 , 2 , ⋯ P\{N(t+s)-N(s)=n\}=\mathrm{e}^{-\lambda t} \frac{(\lambda t)^{n}}{n !}, \quad n=0,1,2, \cdots P{N(t+s)−N(s)=n}=e−λtn!(λt)n​,n=0,1,2,⋯

从泊松过程的定义可以看到，$N(t+s)-N(s)$ 的分布为 P { N ( t + s ) − N ( s ) = n } = e − λ t ( λ t ) n n ! , n = 0 , 1 , 2 , ⋯ P\{N(t+s)-N(s)=n\}=\mathrm{e}^{-\lambda t} \dfrac{(\lambda t)^{n}}{n !}, \quad n=0,1,2, \cdots P{N(t+s)−N(s)=n}=e−λtn!(λt)n​,n=0,1,2,⋯与所在时间点t无关，只与两个时刻的经历的时间段有关，因此蕴含了泊松过程是平稳增量过程（具体验证和证明后面再做）。还是由 P { N ( t + s ) − N ( s ) = n } = e − λ t ( λ t ) n n ! , n = 0 , 1 , 2 , ⋯ P\{N(t+s)-N(s)=n\}=\mathrm{e}^{-\lambda t} \frac{(\lambda t)^{n}}{n !}, \quad n=0,1,2, \cdots P{N(t+s)−N(s)=n}=e−λtn!(λt)n​,n=0,1,2,⋯可以看到，一段时间段$t$内，事件发生的次数服从均值为 $\lambda t$ 的 Poisson 分布，即$$E[N(t)]=\lambda t$$，因此$\lambda$代表了某个时间段内时间的平均发生次数，简称发生率。

【排队论应用】设某银行从早上8:30营业，到中午12:30休息，顾客以25人/h的速度进入银行办业务，则在11:30～12:30内没有客户办理业务的概率是多少？

解：假设客户到银行办理业务的过程为泊松过程，设早上8:30开始为$t=0$时刻，那么11:30为$t=3$时刻，12:30为$t=4$时刻，顾客以25人/h的速度进入银行办业务代表了$\lambda =25$，因此在11:30～12:30内没有客户办理业务的概率可以表达为：
P { N ( 4 ) − N ( 3 ) = 0 } = e − 25 ( 25 ) 0 0 ! = e − 25 P\{N(4)-N(3)=0\}=\mathrm{e}^{-25} \frac{(25)^{0}}{0 !}=\mathrm{e}^{-25} P{N(4)−N(3)=0}=e−250!(25)0​=e−25

为什么很多描述事物出现的次数的事件能用泊松过程来描述呢？
小概率事件原理。一般来说，事物发生的次数应该使用二项分布来描述，即假设事物发生的概率为$p$，那么$N$次实验中事物发生$x$次应该服从二项分布，那为什么这里可以使用泊松分布来描述呢？这是因为二项分布在某些情况下会逼近泊松分布，为了描述二项分布逼近泊松分布，我们需要来看看什么情况下二项分布会逼近泊松分布。

# 二项分布逼近泊松分布
import seaborn as sns 
def binormial_vs_poisson(n, p):
    samples_times = 10000
    binormial_samples = np.random.binomial(n=[n]*samples_times, p=p)
    poisson_samples = np.random.poisson(lam=n*p, size=samples_times)
    sns.histplot(binormial_samples, label='binormail', kde=True,  alpha=0.75, color='orange', stat='probability')
    sns.histplot(poisson_samples, label='poisson', kde=True,  alpha=0.75, stat='probability')
    plt.xlabel("X")
    plt.ylabel("p(x)")
    plt.legend()
    plt.plot()
binormial_vs_poisson(100, 0.5)  
binormial_vs_poisson(1000, 0.3) 
1
2
3
4
5
6
7
8
9
10
11
12
13
14

从上述实验可以看到，当二项分布中的每个伯努利实验的事件发生的概率很小，而实验次数很多的话，二项分布会更加接近于泊松分布，如：道路上发生交通事故的概率很小等等。这个近似想法推广到随机过程，那就是虽然在很短时间内单次事件发生的概率很小，但是如果很多个很短的时间连接在一起，小概率事件的发生将会以一个稳定的速率进行，这个描述很接近平稳增量过程的描述，因此可以引出泊松过程的第二个定义：

（2）泊松过程的第二个定义

（定义二）满足以下条件 的计数过程 { N ( t ) , t ⩾ 0 } \{N(t), t \geqslant 0\} {N(t),t⩾0} 是 Poisson 过程，即：

• (1)$N(0)=0$；
• (2)过程有平稳独立增量;
• (3)存在 $\lambda >0$, 当 $h\to 0$ 时, 有
  P { N ( t + h ) − N ( t ) = 1 } = λ h + o ( h ) P\{N(t+h)-N(t)=1\}=\lambda h+o(h) P{N(t+h)−N(t)=1}=λh+o(h)
• (4) 当 $h\to 0$ 时, 有
  P { N ( t + h ) − N ( t ) ⩾ 2 } = o ( h ) P\{N(t+h)-N(t) \geqslant 2\}=o(h) P{N(t+h)−N(t)⩾2}=o(h)

如何理解以上的几个条件呢？

• 条件（1）说明在过程的一开始，计数为0；
• 条件（2）要求过程有平稳独立增量，目的就是像前文所说的通过连接很多个很短的时间，使得小概率事件以稳定的速率进行，从而达到二项分布近似泊松分布的目的。
• 条件（3）说明在很短时间范围内，事件发生一次的概率很小（h很小），从而满足小概率事件的目的。
• 条件（4）说明了小时间区间范围内事件发生两次及以上的概率几乎为0（几乎不可能）。

除了使用一个时间段内事件发生的次数服从泊松分布这个角度去定义泊松过程以外，我们还能从两次事件发生的时间的分布来定义泊松过程，而这个角度是最适合拿来做模拟实验的角度。如果$T_n$ 表示第 $n$ 次 $(n=1,2,\cdots )$ 事件发生的时刻, 规定 $T_0=0$ 。 $X_n$ 表示第 $n$ 次与 第 $n-1$ 次事件发生的时间间隔，那么$X_n(n=1,2,\cdots )$ 服从参数为 $\lambda$ 的指数分布, 且相互独立。

（3）泊松分布的第三个定义

（定义三）如果每次事件发生的时间间隔$X_1,X_2,...$相互独立，且服从同一个参数为$\lambda$的指数分布，那么计数过程 { N ( t ) , t ⩾ 0 } \{N(t), t \geqslant 0\} {N(t),t⩾0} 是参数为 $\lambda$ 的 Poisson 过程。

# 使用定义三模拟泊松过程的一条样本轨迹
def simulate_poisson(lmd, n_steps):
    exp_samples = np.random.exponential(scale=1./lmd, size=n_steps)
    t_samples = np.cumsum(exp_samples)
    return t_samples
lmd=10
n_steps=10
t_samples = simulate_poisson(lmd=lmd, n_steps=n_steps)
count_x = 0
plt.figure(figsize=(8,6))
for step in range(n_steps-1):
    count_x += 1
    plt.plot(t_samples[step:step+2], [count_x]*2, color='red', alpha=0.75)
    plt.scatter(t_samples[step], [count_x], color='red', s=20)
    plt.scatter(t_samples[step+1], [count_x], color='white', s=20, linewidths=2, edgecolors='gray')
plt.yticks(np.arange(10))
plt.xlabel("t")
plt.ylabel("N(t)")
plt.title("模拟泊松过程的样本轨迹")
plt.plot()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

2.3 泊松过程的合成与分解

有的时候，在现实生活中不只有一个泊松过程，如：假设有一个车站，车站内有两种公交车，分别为：3号和15号，如果车站内为了乘坐3号线的乘客到来服从泊松过程 { N 1 ( t ) } \{N_1(t)\} {N1​(t)}，车站内为了乘坐15号线的乘客到来服从泊松过程 { N 2 ( t ) } \{N_2(t)\} {N2​(t)}，现在问为了坐车而来的乘客是一个什么样的过程，即 { X ( t ) = N 1 ( t ) + N 2 ( t ) } \{X(t) = N_1(t)+N_2(t)\} {X(t)=N1​(t)+N2​(t)}是一个什么过程？这就涉及到将两个乃至于多个泊松过程的合成过程，因此：

（1）泊松过程的合成：
设 { N 1 ( t ) } \left\{N_{1}(t)\right\} {N1​(t)} 和 { N 2 ( t ) } \left\{N_{2}(t)\right\} {N2​(t)} 是强度为 ${\lambda }_1$ 和 ${\lambda }_2$ 的泊松过程， 且相互独立, 则 { N 1 ( t ) + N 2 ( t ) } \left\{N_{1}(t)+N_{2}(t)\right\} {N1​(t)+N2​(t)} 是 ${\lambda }_1+{\lambda }_2$ 的泊松过程。

反过来，车站内为了乘车的乘客到来服从强度为$\lambda$的泊松过程，而每位乘客选择3号线的概率为$p$，选择15号线的概率为$1-p$。因此，乘坐三号线的乘客到来服从强度为$\lambda p$的泊松过程，乘坐15号线的乘客到来服从强度为$(1-p)\lambda$的泊松过程。具体来说：

（2）泊松过程的分解：设 { N ( t ) } \{N(t)\} {N(t)} 是强度为 $\lambda$ 的泊松过程， 若每个事件独立地以概率 $p$ 为类型 1 ， 以 $1-p$ 为类型 2 ， 令 { N 1 ( t ) } \left\{N_{1}(t)\right\} {N1​(t)} 和 { N 2 ( t ) } \left\{N_{2}(t)\right\} {N2​(t)} 分别表示到 $t$ 为止类型1和 类型2发生的个数， 则 { N 1 ( t ) } \left\{N_{1}(t)\right\} {N1​(t)} 和 { N 2 ( t ) } \left\{N_{2}(t)\right\} {N2​(t)} 分别是强度为 $\lambda p$ 和 $\lambda (1-p)$ 的泊松过程， 且相互独立。

【案例】设有一个交通十字路口，车辆分别从东南西北驶进十字路口服从强度为5，6，7，8的泊松过程，每辆车在经过十字路口时直行、左转、右转和掉头的概率分别为：0.5，0.2，0.2和0.1，请问驶出十字路口的各个方向的泊松过程的强度是多少？

解决办法如下：

2.4 复合泊松过程

除了多个泊松过程的合成和分解，有时候我们需要研究多个步骤的“复合”，如：顾客先进入商店，再选择买东西。

复合泊松过程的一个场景：假设顾客进入商店的过程是一个泊松过程，这个客户进入商店后买东西的金额是服从某个分布的随机变量，现在想要研究某个时间段内商店的营业额。这个问题是典型的复合泊松过程问题，具体来说：

{ N ( t ) , t ⩾ 0 } \{N(t), t \geqslant 0\} {N(t),t⩾0} 是强度为 $\lambda$ 的泊松过程, { Y k , k = 1 , 2 , ⋯   } \left\{Y_{k}, k=1,2, \cdots\right\} {Yk​,k=1,2,⋯} 是一 列独立同分布随机变量， 且与 { N ( t ) , t ⩾ 0 } \{N(t), t \geqslant 0\} {N(t),t⩾0} 独立， 令
$$X(t)=\sum _{k=1}^{N(t)}{Y}_k,t⩾0,$$
则称 { X ( t ) , t ⩾ 0 } \{X(t), t \geqslant 0\} {X(t),t⩾0} 为复合泊松过程。

在复合泊松过程中，我们可能更加关注的是$X(t)$的均值和方差，即：若 $$\mathrm{E}\left({\mathrm{Y}}_{\mathrm{i}}^2\right)<+\infty$$, 则 $$\mathrm{E}[\mathrm{X}(\mathrm{t})]=\lambda \mathrm{tE}\left({\mathrm{Y}}_1\right),\mathrm{Var}[\mathrm{X}(\mathrm{t})]=\lambda \mathrm{t}E\left({\mathrm{Y}}_1^2\right)$$

2.5 泊松过程的案例

（1）道路口的车辆数统计

目标：掌握如何使用简单的simpy进行仿真的流程。

案例说明：为了改善道路的路面情况（道路经常维修，坑坑洼洼），因此想统计一天中有多少车辆经过，因为每天的车辆数都是随机的，一般来说有两种技术解决这个问题：

• （1）在道路附近安装一个计数器或安排一个技术人员，在一段长时间的天数（如365天）每天24h统计通过道路的车辆数。
• （2）使用仿真技术大致模拟下道路口的场景，得出一个近似可用的仿真统计指标。

由于方案（1）需要花费大量的人力物力以及需要花费大量的调研时间，虽然能得出准确的结果，但是有时候在工程应用中并不允许。因此，我们选择方案（2），我们通过一周的简单调查，得到每天的每个小时平均车辆数：[30, 20, 10, 6, 8, 20, 40, 100, 250, 200, 100, 65, 100, 120, 100, 120, 200, 220, 240, 180, 150, 100, 50, 40]，通过利用平均车辆数进行仿真：

# 模拟仿真研究该道路口一天平均有多少车经过
import simpy

class Road_Crossing:
    def __init__(self, env):
        self.road_crossing_container = simpy.Container(env, capacity = 1e8, init = 0)
    
def come_across(env, road_crossing, lmd):
    while True:
        body_time = np.random.exponential(1.0/(lmd/60))  # 经过指数分布的时间后，泊松过程记录数+1
        yield env.timeout(body_time)  # 经过body_time个时间
        yield road_crossing.road_crossing_container.put(1)

hours = 24  # 一天24h
minutes = 60  # 一个小时60min
days = 3   # 模拟3天
lmd_ls = [30, 20, 10, 6, 8, 20, 40, 100, 250, 200, 100, 65, 100, 120, 100, 120, 200, 220, 240, 180, 150, 100, 50, 40]   # 每隔小时平均通过车辆数
car_sum = []  # 存储每一天的通过路口的车辆数之和
print('仿真开始：')
for day in range(days):
    day_car_sum = 0   # 记录每天的通过车辆数之和
    for hour, lmd in enumerate(lmd_ls):
        env = simpy.Environment()
        road_crossing = Road_Crossing(env)
        come_across_process = env.process(come_across(env, road_crossing, lmd))
        env.run(until = 60)  # 每次仿真60min
        if hour % 4 == 0:
            print("第"+str(day+1)+"天，第"+str(hour+1)+"时的车辆数：", road_crossing.road_crossing_container.level)
        day_car_sum += road_crossing.road_crossing_container.level
    car_sum.append(day_car_sum)
print("每天通过交通路口的的车辆数之和为：", car_sum)
'''
仿真开始：
第1天，第1时的车辆数： 23
第1天，第5时的车辆数： 15
第1天，第9时的车辆数： 275
第1天，第13时的车辆数： 105
第1天，第17时的车辆数： 223
第1天，第21时的车辆数： 160
第2天，第1时的车辆数： 31
第2天，第5时的车辆数： 5
第2天，第9时的车辆数： 246
第2天，第13时的车辆数： 102
第2天，第17时的车辆数： 214
第2天，第21时的车辆数： 150
第3天，第1时的车辆数： 21
第3天，第5时的车辆数： 5
第3天，第9时的车辆数： 234
第3天，第13时的车辆数： 94
第3天，第17时的车辆数： 203
第3天，第21时的车辆数： 154
每天通过交通路口的的车辆数之和为： [2571, 2405, 2408]
'''
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53

（2）复合泊松过程案例

仿真“每天的商店营业额”这个复合泊松过程：

首先，我们假设每个小时进入商店的平均人数为：[10, 5, 3, 6, 8, 10, 20, 40, 100, 80, 40, 50, 100, 120, 30, 30, 60, 80, 100, 150, 70, 20, 20, 10]，每位顾客的平均花费为：10元（大约一份早餐吧），请问每天商店的营业额是多少？

# 模拟仿真研究该商店一天的营业额
import simpy

class Store_Money:
    def __init__(self, env):
        self.store_money_container = simpy.Container(env, capacity = 1e8, init = 0)
    
def buy(env, store_money, lmd, avg_money):
    while True:
        body_time = np.random.exponential(1.0/(lmd/60))  # 经过指数分布的时间后，泊松过程记录数+1
        yield env.timeout(body_time) 
        money = np.random.poisson(lam=avg_money)
        yield store_money.store_money_container.put(money)

hours = 24  # 一天24h
minutes = 60  # 一个小时60min
days = 3   # 模拟3天
avg_money = 10
lmd_ls = [10, 5, 3, 6, 8, 10, 20, 40, 100, 80, 40, 50, 100, 120, 30, 30, 60, 80, 100, 150, 70, 20, 20, 10]   # 每个小时平均进入商店的人数
money_sum = []  # 存储每一天的商店营业额总和
print('仿真开始：')
for day in range(days):
    day_money_sum = 0   # 记录每天的营业额之和
    for hour, lmd in enumerate(lmd_ls):
        env = simpy.Environment()
        store_money = Store_Money(env)
        store_money_process = env.process(buy(env, store_money, lmd, avg_money))
        env.run(until = 60)  # 每次仿真60min
        if hour % 4 == 0:
            print("第"+str(day+1)+"天，第"+str(hour+1)+"时的营业额：", store_money.store_money_container.level)
        day_money_sum += store_money.store_money_container.level
    money_sum.append(day_money_sum)
print("每天商店的的营业额之和为：", money_sum)
'''
仿真开始：
第1天，第1时的营业额： 83
第1天，第5时的营业额： 42
第1天，第9时的营业额： 1028
第1天，第13时的营业额： 917
第1天，第17时的营业额： 598
第1天，第21时的营业额： 797
第2天，第1时的营业额： 139
第2天，第5时的营业额： 150
第2天，第9时的营业额： 1054
第2天，第13时的营业额： 1046
第2天，第17时的营业额： 473
第2天，第21时的营业额： 823
第3天，第1时的营业额： 138
第3天，第5时的营业额： 22
第3天，第9时的营业额： 996
第3天，第13时的营业额： 1076
第3天，第17时的营业额： 580
第3天，第21时的营业额： 666
每天商店的的营业额之和为： [11691, 11823, 12284]
'''
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55

三、马尔可夫过程

现在不谈时间序列分析，想看看基于什么假设下，我们比较容易“预测”未来，这个假设叫做——马尔可夫假设。马尔可夫假设说的是假设未来只与现在有关，与过去是没有关系的，这个马尔可夫假设在某些情况下是对现实的合理简化。

栗子：要预测某高考生被啥级别大学录取，不需要知道初中小学的努力和成绩，知道高考成绩即可。即马尔可夫假设的一个合理解释。人生就像是一个马尔可夫链，你的未来取决于你当下正在做的事，当下做的决定，而与过去做完的事无关，骄兵必败。

3.1 马尔可夫性

马尔可夫性质是因俄国数学家安德烈·马尔可夫而得名，其含义是：当一个随机过程在给定现在状态及所有过去状态情况下，其未来状态的条件概率分布仅依赖于当前状态；

即在给定现在状态时，它与过去状态（即该过程的历史路径）是条件独立的，那么此随机过程即具有马尔可夫性质。马尔可夫假设则是用来假设一个模型具有马尔可夫性质，比如隐马尔可夫模型。具体来说，马尔可夫性是指：

设 { X ( t ) , t ∈ T } \{X(t), t \in T\} {X(t),t∈T} 为一随机过程， $E$ 为其状态空间，若对任意的 $t_1<t_2<\cdots <t_n<t$ ，任意的 $x_1,x_2,\cdots ,x_n,x,\in E$ ，随机变量 $\mathrm{X}(t)$ 在已知变量 $X\left(t_1\right)=x_1,\cdots ,X\left(t_n\right)=x_n$ 之下的条件分布函数只与 $X\left(t_n\right)=x_n$ 有关，而与 $X\left(t_1\right)=x_1,\cdots ,X\left(t_{n-1}\right)=x_{n-1}$ 无关，即条件分布函数满足等式
$$F\left(x,t\mid x_n,x_{n-1},\cdots ,x_2,x_1,t_n,t_{n-1},\cdots ,t_2,t_1\right)=F\left(x,t\mid x_n,t_n\right)$$
即
P { X ( t ) ⩽ x ∣ X ( t n ) = x n , ⋯   , X ( t 1 ) = x 1 } = P { X ( t ) ⩽ x ∣ X ( t n ) = x n } P\left\{X(t) \leqslant x \mid X\left(t_{n}\right)=x_{n}, \cdots, X\left(t_{1}\right)=x_{1}\right\}=P\left\{X(t) \leqslant x \mid X\left(t_{n}\right)=x_{n}\right\} P{X(t)⩽x∣X(tn​)=xn​,⋯,X(t1​)=x1​}=P{X(t)⩽x∣X(tn​)=xn​}
此性质称为马尔可夫性，亦称无后效性或无记忆性。

3.2 初始分布、转移概率与转移概率矩阵

栗子：社会学家通过长时间的观察与研究发现：决定一个人的收入阶层最重要的因素就是其父母的收入阶层。现在，社会学家将社会的人口分成三个层级：上层人、中层人和底层人，分别用1、2、3来表示。通过研究父代和子代的社会阶层关系，社会学家发现：如果一个人的收入属于底层人，则它的孩子属于底层收入的概率为0.65，属于中层收入的概率为0.28，属于上层收入的概率为0.07等等，可以将规律汇总：

假设现在的社会阶层的分布为：$[0.110,0.680,0.210]$，请问N年后社会阶层的分布是如何？

在上面的案例中，社会在还没开始衍化时的阶层分布为初始分布，即：$[0.110,0.680,0.210]$;父代状态为1（上层），子代状态为1（上层）的概率为0.52，这个概率就是转移概率，记作：$p_{11}$，代表由状态1一步转移到状态1的概率。按照这个描述，应该还有$p_{12},p_{13},p_{21},...,p_{33}$，这些都是转移概率。将所有的转移概率形成一个矩阵，即：

这个将所有的转移概率囊括其中的概率矩阵称为转移概率矩阵，转移概率矩阵的每一行的概率和为1，这是显然的，因为由某个状态转化到所有已知状态的概率和肯定是1。

问：“假设现在的社会阶层的分布为：$[0.110,0.680,0.210]$，请问N年后社会阶层的分布是如何？
可以使用python模拟这个过程：

# 模拟N年后的社会阶层的分布情况
def get_years_dist(p0, P, N):
    P1 = P
    for i in range(N):
        P1 = np.matmul(P1, P)
    return np.matmul(p0, P1)

p0 = np.array([0.110, 0.680, 0.210])
P = np.array([
    [0.52, 0.36, 0.12],
    [0.18, 0.67, 0.15],
    [0.07, 0.28, 0.65]
])
print("初始阶层分布为",p0)
N = 20
for year in range(N):
    print("第"+str(year+1)+"年后的阶层分布为",np.round(get_years_dist(p0, P, year+1), 6)
'''
初始阶层分布为 [0.11 0.68 0.21]
第1年后的阶层分布为 [0.218375 0.511604 0.270021]
第2年后的阶层分布为 [0.224545 0.496996 0.278459]
第3年后的阶层分布为 [0.225715 0.491792 0.282493]
第4年后的阶层分布为 [0.225669 0.489856 0.284475]
第5年后的阶层分布为 [0.225435 0.489097 0.285468]
第6年后的阶层分布为 [0.225247 0.488783 0.285971]
第7年后的阶层分布为 [0.225127 0.488645 0.286228]
第8年后的阶层分布为 [0.225058 0.488582 0.28636 ]
第9年后的阶层分布为 [0.22502  0.488552 0.286428]
第10年后的阶层分布为 [0.225    0.488537 0.286464]
第11年后的阶层分布为 [0.224989 0.488529 0.286482]
第12年后的阶层分布为 [0.224983 0.488526 0.286491]
第13年后的阶层分布为 [0.22498  0.488524 0.286496]
第14年后的阶层分布为 [0.224979 0.488523 0.286499]
第15年后的阶层分布为 [0.224978 0.488522 0.2865  ]
第16年后的阶层分布为 [0.224977 0.488522 0.286501]
第17年后的阶层分布为 [0.224977 0.488522 0.286501]
第18年后的阶层分布为 [0.224977 0.488522 0.286501]
第19年后的阶层分布为 [0.224977 0.488522 0.286501]
第20年后的阶层分布为 [0.224977 0.488522 0.286501]
'''
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40

从上面看出马尔科夫过程的转化规律：马尔可夫过程只需要知道初始分布以及转移概率矩阵就可以知道往后任意一时间点的状态分布，计算规律非常简单，只需要使用初始分布与转移概率矩阵的N次幂的矩阵乘法，N代表了经过了N个时间点。

3.3 平稳分布

观察N年的阶层分布的转化规律，是不是发现在转化了16年后，阶层分布开始稳定在一个确定的分布：$[\mathrm{0.2249770.4885220.286501}]$，这是偶然吗？我们来用一个实验来验证下这个观点。在这个案例，我们将使用三个不同的初始分布：$[0.1,0.8,0.1]$，$[0.8,0.1,0.1]$和$[0.1,0.1,0.8]$，观察在三个初始分布不同情况下，最终是否会有稳定的分布？

# 使用3个不同的初始分布，模拟N年后的社会阶层的分布情况
def get_years_dist(p0, P, N):
    P1 = P
    for i in range(N):
        P1 = np.matmul(P1, P)
    return np.matmul(p0, P1)

p1 = np.array([0.1, 0.8, 0.1])
p2 = np.array([0.8, 0.1, 0.1])
p3 = np.array([0.1, 0.1, 0.8])
p_list = [p1, p2, p3]
P = np.array([
    [0.52, 0.36, 0.12],
    [0.18, 0.67, 0.15],
    [0.07, 0.28, 0.65]
])

for i, p_i in enumerate(p_list):
    print("+++++++++++++++p"+str(i+1)+"的实验+++++++++++++++++")
    print("初始阶层分布为",p_i)
    N = 30
    for year in range(N):
        print("第"+str(year+1)+"年后的阶层分布为",np.round(get_years_dist(p_i, P, year+1), 6))
'''
+++++++++++++++p1的实验+++++++++++++++++
初始阶层分布为 [0.1 0.8 0.1]
第1年后的阶层分布为 [0.22735 0.53024 0.24241]
第2年后的阶层分布为 [0.230634 0.504982 0.264385]
第3年后的阶层分布为 [0.229333 0.495394 0.275273]
第4年后的阶层分布为 [0.227693 0.49155  0.280757]
第5年后的阶层分布为 [0.226532 0.48992  0.283548]
第6年后的阶层分布为 [0.225831 0.489191 0.284978]
第7年后的阶层分布为 [0.225435 0.488851 0.285714]
第8年后的阶层分布为 [0.225219 0.488687 0.286094]
第9年后的阶层分布为 [0.225104 0.488605 0.28629 ]
第10年后的阶层分布为 [0.225043 0.488564 0.286392]
第11年后的阶层分布为 [0.225012 0.488544 0.286445]
第12年后的阶层分布为 [0.224995 0.488533 0.286472]
第13年后的阶层分布为 [0.224986 0.488527 0.286486]
第14年后的阶层分布为 [0.224982 0.488525 0.286493]
第15年后的阶层分布为 [0.22498  0.488523 0.286497]
第16年后的阶层分布为 [0.224978 0.488522 0.286499]
第17年后的阶层分布为 [0.224978 0.488522 0.2865  ]
第18年后的阶层分布为 [0.224977 0.488522 0.286501]
第19年后的阶层分布为 [0.224977 0.488522 0.286501]
第20年后的阶层分布为 [0.224977 0.488522 0.286501]
第21年后的阶层分布为 [0.224977 0.488522 0.286501]
第22年后的阶层分布为 [0.224977 0.488522 0.286501]
第23年后的阶层分布为 [0.224977 0.488522 0.286501]
第24年后的阶层分布为 [0.224977 0.488522 0.286501]
第25年后的阶层分布为 [0.224977 0.488522 0.286501]
第26年后的阶层分布为 [0.224977 0.488522 0.286501]
第27年后的阶层分布为 [0.224977 0.488522 0.286501]
第28年后的阶层分布为 [0.224977 0.488522 0.286501]
第29年后的阶层分布为 [0.224977 0.488522 0.286501]
第30年后的阶层分布为 [0.224977 0.488522 0.286501]
+++++++++++++++p2的实验+++++++++++++++++
初始阶层分布为 [0.8 0.1 0.1]
第1年后的阶层分布为 [0.31058 0.46465 0.22477]
第2年后的阶层分布为 [0.260872 0.48606  0.253068]
第3年后的阶层分布为 [0.240859 0.490433 0.268708]
第4年后的阶层分布为 [0.232334 0.490538 0.277128]
第5年后的阶层分布为 [0.22851  0.489896 0.281594]
第6年后的阶层分布为 [0.226718 0.48934  0.283942]
第7年后的阶层分布为 [0.225851 0.48898  0.285169]
第8年后的阶层分布为 [0.225421 0.48877  0.285809]
第9年后的阶层分布为 [0.225204 0.488654 0.286142]
第10年后的阶层分布为 [0.225094 0.488591 0.286315]
第11年后的阶层分布为 [0.225037 0.488558 0.286405]
第12年后的阶层分布为 [0.225008 0.488541 0.286451]
第13年后的阶层分布为 [0.224993 0.488532 0.286475]
第14年后的阶层分布为 [0.224985 0.488527 0.286488]
第15年后的阶层分布为 [0.224981 0.488524 0.286494]
第16年后的阶层分布为 [0.224979 0.488523 0.286498]
第17年后的阶层分布为 [0.224978 0.488522 0.286499]
第18年后的阶层分布为 [0.224978 0.488522 0.2865  ]
第19年后的阶层分布为 [0.224977 0.488522 0.286501]
第20年后的阶层分布为 [0.224977 0.488522 0.286501]
第21年后的阶层分布为 [0.224977 0.488522 0.286501]
第22年后的阶层分布为 [0.224977 0.488522 0.286501]
第23年后的阶层分布为 [0.224977 0.488522 0.286501]
第24年后的阶层分布为 [0.224977 0.488522 0.286501]
第25年后的阶层分布为 [0.224977 0.488522 0.286501]
第26年后的阶层分布为 [0.224977 0.488522 0.286501]
第27年后的阶层分布为 [0.224977 0.488522 0.286501]
第28年后的阶层分布为 [0.224977 0.488522 0.286501]
第29年后的阶层分布为 [0.224977 0.488522 0.286501]
第30年后的阶层分布为 [0.224977 0.488522 0.286501]
+++++++++++++++p3的实验+++++++++++++++++
初始阶层分布为 [0.1 0.1 0.8]
第1年后的阶层分布为 [0.16267 0.41761 0.41972]
第2年后的阶层分布为 [0.189139 0.455882 0.35498 ]
第3年后的阶层分布为 [0.205259 0.472925 0.321816]
第4年后的阶层分布为 [0.214388 0.480861 0.30475 ]
第5年后的阶层分布为 [0.21937  0.484687 0.295943]
第6年后的阶层分布为 [0.222032 0.486578 0.291391]
第7年后的阶层分布为 [0.223438 0.487528 0.289034]
第8年后的阶层分布为 [0.224175 0.488011 0.287814]
第9年后的阶层分布为 [0.22456  0.488258 0.287182]
第10年后的阶层分布为 [0.22476  0.488386 0.286854]
第11年后的阶层分布为 [0.224865 0.488451 0.286684]
第12年后的阶层分布为 [0.224919 0.488485 0.286596]
第13年后的阶层分布为 [0.224947 0.488503 0.286551]
第14年后的阶层分布为 [0.224961 0.488512 0.286527]
第15年后的阶层分布为 [0.224969 0.488517 0.286515]
第16年后的阶层分布为 [0.224973 0.488519 0.286508]
第17年后的阶层分布为 [0.224975 0.48852  0.286505]
第18年后的阶层分布为 [0.224976 0.488521 0.286503]
第19年后的阶层分布为 [0.224976 0.488521 0.286502]
第20年后的阶层分布为 [0.224977 0.488521 0.286502]
第21年后的阶层分布为 [0.224977 0.488521 0.286502]
第22年后的阶层分布为 [0.224977 0.488522 0.286502]
第23年后的阶层分布为 [0.224977 0.488522 0.286501]
第24年后的阶层分布为 [0.224977 0.488522 0.286501]
第25年后的阶层分布为 [0.224977 0.488522 0.286501]
第26年后的阶层分布为 [0.224977 0.488522 0.286501]
第27年后的阶层分布为 [0.224977 0.488522 0.286501]
第28年后的阶层分布为 [0.224977 0.488522 0.286501]
第29年后的阶层分布为 [0.224977 0.488522 0.286501]
第30年后的阶层分布为 [0.224977 0.488522 0.286501]
'''
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121

无论初始分布如何，只要转移概率矩阵不变，最终整个社会的阶层分布都会是$[\mathrm{0.2249770.4885220.286501}]$。因此，$[\mathrm{0.2249770.4885220.286501}]$称为该转移概率矩阵的平稳分布，平稳分布与初始分布是没有关系的，只与转移概率矩阵有关。

如果存在 状态空间$E$ 上的概率分布 $\pi =\left({\pi }_0,{\pi }_{1,\dots }\right)$, 满足矩阵方程
$$\pi =\pi P$$
则称 $\pi$ 为马尔可夫链的平稳分布。

平稳分布：某个分布如果经过转移概率矩阵的转化后还是这个分布不变，那该分布就是对应转移概率矩阵的平稳分布。

最后，我们从现实的意义来评价下这个社会阶层衍化的案例：社会阶层的转移概率矩阵相当于各阶层的人跃迁的通道，在现实的角度类似于国民教育，而平稳分布于初始分布是无关的，只与转移概率矩阵是有关系的。如果zf不加以干预，让某种不好的转化延续下去，将导致社会阶层的固化。

3.4 马尔可夫过程的案例

艾滋病发展过程分为四个阶段（状态），急性感染期（状态 1）、无症状期（状态 2）， 艾滋病前期（状态 3）, 典型艾滋病期（状态 4）。艾滋病发展过程基本上是一个不可逆的过程,即：状态1 -> 状态2 -> 状态3 -> 状态4。现在收集某地600例艾滋病防控数据，得到以下表格：

若一个人此时是无症状期（状态2）在10次转移之后，这个人的各状态的概率是多少？

# 研究无症状期病人在10期转移后的状态分布
def get_years_dist(p0, P, N):
    P1 = P
    for i in range(N):
        P1 = np.matmul(P1, P)
    return np.matmul(p0, P1)

p0 = np.array([0, 1, 0, 0])
P = np.array([
    [10.0/80, 62.0/80, 5.0/80, 3.0/80],
    [0, 140.0/290, 93.0/290, 57.0/290],
    [0, 0, 180.0/220, 40.0/220],
    [0, 0, 0, 1]
])
N = 10
print(str(N)+"期转移后，状态分布为：", np.round(get_years_dist(p0, P, N), 4))
# 10期转移后，状态分布为： [0.000e+00 3.000e-04 1.048e-01 8.948e-01]
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17

四、实战建模

4.1 泊松过程实战：十字路口交通压力仿真（泊松过程的合成和分解）

复现之前在“泊松过程的合成与分解”章节的十字路口的例子，先来回顾下：

设有一个交通十字路口，车辆分别从东南西北驶进十字路口服从强度为5，6，7，8的泊松过程，每辆车在经过十字路口时直行、左转、右转和掉头的概率分别为：0.5，0.2，0.2和0.1，请问驶出十字路口的各个方向的泊松过程的强度是多少？

前面的计算是使用结论的方式计算的，这次用仿真。

class CrossRoad:
    def __init__(self, env):
        self.left = simpy.Container(env, capacity=1e8, init=0)  # 记录左转的次数
        self.right = simpy.Container(env, capacity=1e8, init=0)  # 记录右转的次数
        self.straight = simpy.Container(env, capacity=1e8, init=0)   # 记录直行的次数
        self.turn = simpy.Container(env, capacity=1e8, init=0)   # 记录掉头的次数


def car_cross(env, cross_road, lmd):    # 模拟从某个方向来的车辆在左转、右转、直走和掉头的次数
    while True:
        yield env.timeout(np.random.exponential(scale=1. / lmd))
        east_choice = np.random.choice(['left', 'right', 'straight', 'turn'], p=[0.2, 0.2, 0.5, 0.1])
        if east_choice == 'left':
            cross_road.left.put(1)
        if east_choice == 'right':
            cross_road.right.put(1)
        if east_choice == 'straight':
            cross_road.straight.put(1)
        if east_choice == 'turn':
            cross_road.turn.put(1)


if __name__ == "__main__":
    car_num = np.zeros([4, 4])  # 记录从🧍每个方向到各个方向的车流量，行代表从某个方向出去，列代表到达哪个方向
    total_times = 60 * 24 * 100  # 60min * 24hour * 100day
    for i, (direction, lmd) in enumerate(zip(["east", "south", "west", "north"], [5.0, 6.0, 7.0, 8.0])):
        env = simpy.Environment()
        cross_road = CrossRoad(env)
        direction_cross = env.process(car_cross(env, cross_road, lmd))
        env.run(until=total_times)
        if direction == "east":
            car_num[i, 0] = cross_road.turn.level  # east -> east
            car_num[i, 1] = cross_road.left.level  # east -> south
            car_num[i, 2] = cross_road.straight.level   # east -> west
            car_num[i, 3] = cross_road.right.level    # east -> north
        elif direction == "south":
            car_num[i, 0] = cross_road.right.level  # south -> east
            car_num[i, 1] = cross_road.turn.level
            car_num[i, 2] = cross_road.left.level
            car_num[i, 3] = cross_road.straight.level
        elif direction == "west":
            car_num[i, 0] = cross_road.straight.level  # west -> east
            car_num[i, 1] = cross_road.right.level
            car_num[i, 2] = cross_road.turn.level
            car_num[i, 3] = cross_road.left.level
        else:
            car_num[i, 0] = cross_road.left.level  # north -> east
            car_num[i, 1] = cross_road.straight.level
            car_num[i, 2] = cross_road.right.level
            car_num[i, 3] = cross_road.turn.level
    directions_lmd = np.round(np.sum(car_num, axis=0) / total_times, 1)
    print("东路口的lmd=", directions_lmd[0])
    print("南路口的lmd=", directions_lmd[1])
    print("西路口的lmd=", directions_lmd[2])
    print("北路口的lmd=", directions_lmd[3])
'''
东路口的lmd= 6.8
南路口的lmd= 7.0
西路口的lmd= 6.0
北路口的lmd= 6.2
'''
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61

4.2 马尔可夫过程实战：预测上证指数的涨跌问题

用马尔可夫链预测上证指数的点位，这个案例十分有价值。上海证券综合指数简称“上证指数”或“上证综指”，上证指数的样本股是在上海证券交易所全部上市股票，包括A股和B股，因此上证指数反映了上海证券交易所上市股票价格的变动情况，自1991年7月15日起正式发布。下面对上证指数的计算方法做一个简介：

• 在众多股票中抽取多只具有代表性的股票；
• 按单价或总值加权平均，或不加权平均；
• 计算算术平均数、几何平均数，或兼顾价格与总值；

首先，我们需要获取上证指数的历史数据。

import baostock as bs
from pylab import mpl
plt.style.use('ggplot')
# 登陆系统
lg = bs.login()
# 显示登陆返回信息
print('login respond error_code:'+lg.error_code)
print('login respond  error_msg:'+lg.error_msg)

stock_code = 'sh.000001'
stock_name = '上证综指'
rs = bs.query_history_k_data_plus(stock_code,
    "date,code,open,high,low,close,preclose,volume,amount,pctChg",
    start_date='2020-01-02', end_date='2022-06-30', frequency="d")
print('query_history_k_data_plus respond error_code:'+rs.error_code)
print('query_history_k_data_plus respond  error_msg:'+rs.error_msg)
# 打印结果集
data_list = []
while (rs.error_code == '0') & rs.next():
    # 获取一条记录，将记录合并在一起
    data_list.append(rs.get_row_data())
result = pd.DataFrame(data_list, columns=rs.fields)
result['date'] = pd.to_datetime(result['date'])
result['close'] = result['close'].astype("float")

# 绘图
plt.figure(figsize=(8, 6))
plt.plot(result['date'],result['close'],c="#e74c3c")
plt.legend((stock_name,),loc = 'best',fontsize=15)
plt.xlabel('交易日')
plt.xticks()
plt.show()
'''
login success!
login respond error_code:0
login respond  error_msg:success
query_history_k_data_plus respond error_code:0
query_history_k_data_plus respond  error_msg:success
'''
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39

从上证指数的收盘价信息可以知道：自疫情到现在，点数维持在2600点到3800点之间，但是指数在日价格变化是不敏感的，我们可以研究月的变化。因此，我们需要提取每个月末的上证指数的收盘价：

month_close = pd.DataFrame(result.groupby([result['date'].dt.year, result['date'].dt.month])['date', 'close'].last().values, columns=['date', 'close'])
# month_close

# 上证指数自疫情以来月末收盘价的范围：
print("月末最低收盘价：", month_close['close'].min())
print("月末最高收盘价：", month_close['close'].max())
1
2
3
4
5
6

上证指数自疫情以来的月末收盘价为2700～3700之间。按照马尔可夫链的习惯，我们需要将这个范围分成几个离散的状态，因此：

• (2700-2900]之间为状态1
• (2900-3100]之间为状态2
• (3100-3300]之间为状态3
• (3300-3500]之间为状态4
• (3500-3700]之间为状态5

# 获取上证指数的状态
month_close['state'] = 0
month_close.loc[(month_close['close']>2700)&(month_close['close']<=2900), 'state'] = 1
month_close.loc[(month_close['close']>2900)&(month_close['close']<=3100), 'state'] = 2
month_close.loc[(month_close['close']>3100)&(month_close['close']<=3300), 'state'] = 3
month_close.loc[(month_close['close']>3300)&(month_close['close']<=3500), 'state'] = 4
month_close.loc[(month_close['close']>3500)&(month_close['close']<=3700), 'state'] = 5
# month_close

# 需要获取转移概率矩阵，这里我们留个心眼，把后半年的数据作为验证数据，看看马尔可夫链在上证指数的预测上表现如何？
# 分割分析数据集与验证数据集
pred_data = month_close.loc[month_close['date']>='2022-01-01', :]  # 验证数据集
analysis_data = month_close.loc[month_close['date']< '2022-01-01', :]
state_list = [1, 2, 3, 4, 5]
P = np.zeros([len(state_list), len(state_list)])
for state1 in state_list:
    index1 = np.where(analysis_data['state']==state1)[0]
    index1 = np.array([i for i in index1 if i+1 < len(analysis_data)])
    for state2 in state_list:
        P[state1-1, state2-1] = (analysis_data.iloc[(index1+1), -1]==state2).sum()
P_trans = np.zeros_like(P)
for i in range(P.shape[0]):
    for j in range(P.shape[1]):
        P_trans[i, j] = P[i, j] / np.sum(P[i, :], axis=0)
# P_trans
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25

将拿到的转移概率矩阵与初始分布做六个月的转移，即

• 初始分布为：2021-12-31的状态5
• 转移概率矩阵为：P

p0 = np.array([0, 0, 0, 0, 1])
for i in range(6):
    p0 = np.matmul(p0, P_trans)
    print("第"+str(i+1)+"月预测的指数状态分布为：", p0)
    print("第"+str(i+1)+"月真实的指数状态为：",pred_data.iloc[i, :]['state'])
'''
第1月预测的指数状态分布为： [0.         0.         0.         0.28571429 0.71428571]
第1月真实的指数状态为： 4
第2月预测的指数状态分布为： [0.         0.         0.03571429 0.34693878 0.61734694]
第2月真实的指数状态为： 4
第3月预测的指数状态分布为： [0.         0.         0.06122449 0.36771137 0.57106414]
第3月真实的指数状态为： 3
第4月预测的指数状态分布为： [0.         0.         0.07657617 0.37762911 0.54579472]
第4月真实的指数状态为： 2
第5月预测的指数状态分布为： [0.         0.         0.08549172 0.38304399 0.53146429]
第5月真实的指数状态为： 3
第6月预测的指数状态分布为： [0.         0.         0.09062636 0.3861148  0.52325885]
第6月真实的指数状态为： 4
'''
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19

上面只是使用简单的一阶马尔科夫模型，可以使用其他时间序列模型。

附：时间安排

Reference

[1] https://github.com/Git-Model/Modeling-Universe/tree/main/Data-Story
[2] http://www.stat.yale.edu/~pollard/Courses/251.spring2013/
[3] 关于随机过程的其他理解参考https://www.zhihu.com/question/26694486/answer/349872296。