概率论整理(三)

大数定律背后的理论支撑

一组独立同分布的随机变量：(X_1,X_2,X_3,...,X_n)，期望μ，方差(σ^2)，则这组随机变量的均值为

(M_n={X_1+X_2+X_3+...+X_n\over n})

• 用样本均值估计总体均值，为什么可行？

(M_n)的期望E(M_n)=E(X_1+X_2+X_3+...+X_n\over n)

根据数学期望的齐次性和可加性，有

E(M_n)=E(X_1+X_2+X_3+...+X_n\over n)=({1\over n}(EX_1+EX_2+...+EX_n))

由于是独立同分布，则上式等于

(nμ\over n)=μ=E(X_i)

结果我们可以知道，一组独立同分布样本均值的期望就等于随机变量的期望。

• 样本数越多，估计越准

样本均值的方差：var(M_n)

var(M_n)=var(X_1+X_2+X_3+...+X_n\over n)

根据方差的性质，以及是独立分布，有

var(M_n)=var(X_1+X_2+X_3+...+X_n\over n)=({1\over n^2}(varX_1+varX_2+...varX_n))=(nσ^2\over n^2)=(σ^2\over n)

通过推导，我们知道n个随机变量的方差是单一随机变量的(1\over n)，样本均值的方差变小了。随着样本量的增大，样本均值的分布越接近于总体均值。当n趋近于无穷大的时候，也就是说当样本量非常非常大的时候，通过抽样得到的抽样样本去计算所得到的平均值就等于总体均值。

验证大数定律

• 样本均值与随机变量的期望

以15000个服从参数为(10,0.4)的二项分布随机变量为总体，观察随着样本数量增多，样本均值越来越接近分布期望。

这里需要补充一下二项分布的数学期望EX=np，方差VX=np(1-p)，推导方式可以参考概率论整理(二) 中的切比雪夫不等式。

import numpy as np
from scipy.stats import binom
import matplotlib.pyplot as plt

if __name__ == '__main__':

    # 一组实验是进行10次
    n = 10
    # 单次命中概率为0.4
    p = 0.4
    # 总共进行15000次
    sample_size = 15000
    # 期望值
    expected_value = n * p
    N_samples = range(1, sample_size, 10)
    # 进行3组实验
    for k in range(3):
        # 创建一个二项分布对象
        binom_rv = binom(n=n, p=p)
        # 每组实验进行15000次采样,生成15000个二项分布的随机变量
        X = binom_rv.rvs(size=sample_size)
        # 得到前1、11、21...14991个随机变量的均值
        sample_average = [X[:i].mean() for i in N_samples]
        # 画出每组实验的横纵坐标，横坐标为二项分布间隔10次的采样次数，纵坐标为二项分布累计满10次的均值
        plt.plot(N_samples, sample_average, label='average of sample {}'.format(k))
    # 画出数学期望的横纵坐标，横坐标为二项分布间隔10次的采样次数，纵坐标为数学期望
    plt.plot(N_samples, expected_value * np.ones_like(sample_average), ls='--',
             label='true expected value:n*p={}'.format(n * p), c='k')
    plt.legend()
    plt.grid(ls='--')
    plt.show()

从上图可以看到，3次数据采样中，无论是哪一次，它们的累计均值都随着采样次数的增大而逼近二项分布的数学期望。

• 样本均值的方差与分布

100万个服从均值为0，标准差为20的正态分布随机变量数据，每次从正态分布总体中随机抽取5个样本，计算样本均值，重复1万次观察样本均值的分布；再每次从正态分布总体中随机抽取50个样本，计算样本均值，重复1万次观察样本均值的分布。

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm

if __name__ == '__main__':

    # 创建一个正态分布的随机变量，并对该随机变量采集1000000个数据
    norm_rvs = norm(loc=0, scale=20).rvs(size=1000000)
    # 对这1000000个正态分布的数据画出直方图
    plt.hist(norm_rvs, density=True, alpha=0.3, color='b', bins=100, label='original')

    mean_array = []
    # 进行10000次抽取
    for i in range(10000):
        # 每次从该1000000个正态分布的数据中随机抽取不重复的5个数据
        sample = np.random.choice(norm_rvs, size=5, replace=False)
        # 获取这5个随机数据的均值
        mean_array.append(np.mean(sample))
    # 对这10000次抽取的均值画出直方图
    plt.hist(mean_array, density=True, alpha=0.3, color='r', bins=100, label='sample size=5')

    mean_array = []
    # 进行10000次抽取
    for i in range(10000):
        # 每次从该1000000个正态分布的数据中随机抽取不重复的50个数据
        sample = np.random.choice(norm_rvs, size=50, replace=False)
        # 获取这50个随机数据的均值
        mean_array.append(np.mean(sample))
    # 对这10000次抽取的均值画出直方图
    plt.hist(mean_array, density=True, alpha=0.3, color='g', bins=100, label='sample size=50')

    plt.gca().axes.set_xlim(-60, 60)
    plt.legend(loc='best')
    plt.grid(ls='--')
    plt.show()

上图中蓝色的部分是原始的正态分布数据；红色的部分是每次从原始数据中抽取5个数据，连续抽取10000次得到的均值数据分布；绿色的部分是每次从原始数据中抽取50个数据，连续抽取10000次得到的均值数据分布。从上图我们可以发现，随着每次选取的样本数据增多，样本均值分布的图像越来越向数学期望集中，同时方差不断的减小(分布越来越窄，越来越高)。

蒙特卡罗方法

蒙特卡罗方法的核心思想就是通过模拟出来的大量样本集或者随机过程去近似我们想要研究的问题对象。

典型应用

1. 近似计算不规则面积/体积/积分。
2. 模拟随机过程，预测随机过程可能性结果的区间范围。
3. 利用马尔科夫链-蒙特卡罗方法(MCMC)进行未知参数的统计推断。
4. 用蒙特卡罗方法近似计算圆面积

计算方式为(圆形面积\over 正方形面积)≈(圆内点的个数\over 点的总个数)

这里我们要计算π值，大数定律告诉我们，当点的数量越来越多的时候，π值就会越来越接近于真实值。

(πr^2\over (2r)^2)≈(圆内点的个数\over 点的总个数)可得到π≈4*(圆内点的个数\over 点的总个数)

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Circle
from scipy.stats import uniform

if __name__ == '__main__':

    n = 100000  # 点的总数
    r = 1.0  # 半径
    ox, oy = (0., 0.)  # 圆心
    # 生成点的横纵坐标(均匀分布)
    uniform_x = uniform(ox - r, 2 * r).rvs(n)
    uniform_y = uniform(oy - r, 2 * r).rvs(n)
    # 生成的点到原点的距离，如果小于等于半径则在圆内赋值为1，大于半径则在圆外赋值为0
    d_array = np.sqrt((uniform_x - ox)**2 + (uniform_y - oy)**2)
    res = sum(np.where(d_array <= r, 1, 0))
    pi = (res / n) * 4
    print(pi)

    fig, ax = plt.subplots(1, 1)
    ax.plot(uniform_x, uniform_y, 'ro', alpha=0.2, markersize=0.3)
    plt.axis('equal')
    circle = Circle(xy=(ox, oy), radius=r, alpha=0.5)
    ax.add_patch(circle)
    plt.grid(ls='--')
    plt.show()

运行结果

3.14004

中心极限定理

之前我们说过一维随机变量的中心极限定理

满足

，现在我们来看一下多个独立同分布的随机变量的中心极限定理。

n个满足独立同分布的随机变量(X_1、X_2、X_3...X_n)，期望为μ，方差为(σ^2)，构造新的随机变量(Z_n={X_1+X_2+X_3+...+X_n-nμ\over \sqrt{n}σ})

1. 该随机变量的期望E(Z_n)=E({X_1+X_2+X_3+...+X_n-nμ\over \sqrt{n}σ})=(EX_1+X_2+X_3+...+X_n-nμ\over \sqrt{n}σ)=0
2. 该随机变量的方差var(Z_n)=var({X_1+X_2+X_3+...+X_n-nμ\over \sqrt{n}σ})=(varX_1+varX_2+varX_3+...+varX_n+0\over nσ^2)=1，方差这里需要注意两点：1、常数抽取出来需要变平方；2、常数的方差为0.
3. 随着样本个数n的增大，随机变量的分布逐渐趋向于一个标准正态分布，当n->∞时，随机变量的分布收敛于一个标准正态分布。
4. 这里需要注意的是，中心极限定理对(X_i)随机变量的原始分布没有任何要求，非常具有一般性，(X_i)可以是任意分布。
5. 随机变量(Z_n={X_1+X_2+X_3+...+X_n-nμ\over \sqrt{n}σ})服从标准正态分布。
6. (S_n=X_1+X_2+X_3+...+Xn)趋近于一个均值为nμ，方差为(nσ^2)的正态分布，大量样本的独立随机因素的叠加是趋近于一个正态分布的，不需要知道随机变量的概率分布律(PMF)或者概率密度函数(PDF)。
7. 中心极限定理的模拟与验证

从一个服从参数p=0.3的几何分布中进行采样，分三组实验，分别每次采样2个、5个、50个样本，每组实验各重复100000次，观察(Z_n={X_1+X_2+X_3+...+X_n-nμ\over \sqrt{n}σ})的分布图像。

import numpy as np
from scipy.stats import geom
import matplotlib.pyplot as plt

if __name__ == '__main__':

    fig, ax = plt.subplots(2, 2)
    # 创建一个几何分布对象
    geom_rv = geom(p=0.3)
    # 对该几何分布进行1000000次采样
    geom_rvs = geom_rv.rvs(size=1000000)
    # 获取该几何分布的均值和方差
    mean, var, _, _ = geom_rv.stats(moments='mvsk')
    # 画出原始几何分布的直方图
    ax[0, 0].hist(geom_rvs, bins=100, density=True)
    ax[0, 0].set_title('geom distribution:p=0.3')
    ax[0, 0].grid(ls='--')

    n_array = [0, 2, 5, 50]
    for i in range(1, 4):
        Z_array = []
        n = n_array[i]
        # 分别从原始数据中随机抽取2，5，50个数据，重复100000次
        for j in range(100000):
            sample = np.random.choice(geom_rvs, n)
            # 使用抽取的数据构建新的随机变量Zn
            Z_array.append((sum(sample) - n * mean) / np.sqrt(n * var))
        # 画出每次抽取并构建的新随机变量Zn的直方图
        ax[i // 2, i % 2].hist(Z_array, bins=100, density=True)
        ax[i // 2, i % 2].set_title('n={}'.format(n))
        ax[i // 2, i % 2].set_xlim(-3, 3)
        ax[i // 2, i % 2].grid(ls='--')
    plt.show()

上图中左上角的是原始的几何分布的数据直方图，右上角是每次随机从原始数据中采样2个数据，并重复10万次的(Z_n)；左下角的是每次随机从原始数据中采样5个数据，并重复10万次的(Z_n)；右下角是每次随机从原始数据中采样50个数据，并重复10万次的(Z_n) 。这里每次采样的数据其实就是不同的独立同分布的随机变量(X_1、X_2、X_3...X_n)，随着采集的数据个数增大，样本总量也在增多，在上图中也可以看到随着样本总量的增加，新的随机变量(Z_n)越来越趋近于一个标准正态分布。

随机过程

随机过程是一串有限或无限的随机变量序列。

蒙特卡罗模拟赌博场景

赌徒和庄家赌抛硬币，如果为正面，本轮赌徒胜，庄家付给赌徒1元。如果为反面，本轮庄家胜，赌徒付给庄家1元。赌徒有初始赌本10元，手上的钱一旦输光则退出赌局。

赌博结果依托于每次抛硬币的结果，每一轮赌博就是一个伯努利试验，赢的概率是p=0.5。赌博的过程就是一串伯努利试验构成的随机过程，每轮赌局中赢则赌本增加1元，输则赌本减少1元。

采用蒙特卡罗方法，样本数为100000个赌徒，每个赌徒最多赌10000轮，如果在此过程中输光了赌本，则提前退出，如果到了10000轮还有赌本，赌局也停止。

import pandas as pd
import random

if __name__ == '__main__':

    sample_list = []
    person_num = 100000  # 赌徒人数
    round_num = 10000  # 总轮数

    for person in range(1, person_num + 1):
        money = 10  # 初始赌本
        for round in range(1, round_num + 1):
            result = random.randint(0, 1)  # 随机生成0或1
            if result == 1:
                money += 1
            else:
                money -= 1
            if money == 0:
                break
        sample_list.append([person, round, money])

    sample_df = pd.DataFrame(sample_list, columns=['person', 'round', 'money'])
    sample_df.set_index('person', inplace=True)
    print('总轮数:{}，总人数：{}'.format(round_num, person_num))
    print('输光赌本提前出局的人数:{}'.format(person_num - len(sample_df[sample_df['round'] == round_num])))
    print('赌满全场且盈利的人数：{}'.format(len(sample_df[sample_df['money'] > 10])))
    print('堵满全场且亏损的人数：{}'.format(len(sample_df[sample_df['money'] < 10][sample_df['money'] > 0])))

运行结果(运行时间较长)

总轮数:10000，总人数：100000
输光赌本提前出局的人数:92086
赌满全场且盈利的人数：7873
堵满全场且亏损的人数：24

从以上的结果可以看出，最终能挣钱的不到10%，亏钱人数大于90%。但如果把总轮数调低，比如是100，实验结果大概40%的人可以盈利，60%的人亏损；如果调到1000，实验结果大概20%的人可以盈利，80%的人亏损。试验结果：和庄家1:1的对赌，随着轮数的增加，基本上都破产了。即便双方都不出千，轮数越多，赌徒均亏损，因为庄家的资金量是无穷的。

股价的变化

利用目前的股价(S_t)去预测时间Δt之后的股价(S_{t+1})

金融中有一个公式

(S_{t+1}=S_t+μS_tΔt+σS_tɛ\sqrt{Δt})

其中μ为股票收益率的期望值，设为15%；σ表示股票的波动率，设为0.2；Δt=(T\over n)，T表示整数年份，n表示具体步数，比如T为一年，如果n取244，那么Δt的粒度就是每个交易日(一年有244个交易日)；随机变量ɛ服从标准正态分布，决定了每日股价(S_i)是一个随机变量，而由股价构成的序列是一个随机过程。

假设目前股价为(S_0)=10块钱，用蒙特卡罗方法估计1年后股价概率分布。

import scipy
import matplotlib.pyplot as plt
from math import sqrt

if __name__ == '__main__':

    s0 = 10.0  # 初始股价
    T = 1.0  # 时间1年
    n = 244 * T  # 股票交易日
    mu = 0.15  # 股票收益率期望值
    sigma = 0.2  # 股票波动率
    n_simulation = 10000  # 样本数
    dt = T / n
    s_array = []  # 价格序列
    # 遍历每一个样本
    for i in range(n_simulation):
        s = s0
        # 遍历每一个交易日
        for j in range(int(n)):
            # 标准正态分布的概率值
            e = scipy.random.normal()
            # 获取新的股价
            s = s + mu * s * dt + sigma * s * e * sqrt(dt)
        s_array.append(s)

    plt.hist(s_array, alpha=0.6, bins=30, density=True, edgecolor='k')
    plt.grid(ls='--')
    plt.show()

通过上图我们可以看出相同初始值的10000个样本对应的1年后的股价的分布特征。变成20快和7.5以下的概率都很小，变成12.5的概率是最大的。

• 股票价格变化过程的展现

监测从(T_0)时刻起到1年后的时间段内，每隔Δt时间间隔点由实时价格随机变量构成的序列。由244个时间点的股价数据构成序列，设定100个样本，模拟出100条价格曲线。

import scipy
import matplotlib.pyplot as plt
from math import sqrt
import numpy as np

if __name__ == '__main__':

    s0 = 10.0  # 初始股价
    T = 1.0  # 时间1年
    n = 244 * T  # 股票交易日
    mu = 0.15  # 股票收益率期望值
    sigma = 0.2  # 股票波动率
    n_simulation = 100  # 样本数
    dt = T / n
    random_series = np.zeros(int(n), dtype=float)
    x = range(0, int(n))
    # 遍历每一个样本
    for i in range(n_simulation):
        random_series[0] = s0
        # 遍历每一个交易日
        for j in range(1, int(n)):
            # 标准概率分布随机变量
            e = scipy.random.normal()
            # 记录每一个交易日的新价格
            random_series[j] = random_series[j - 1] + mu * random_series[j - 1] * dt + sigma * random_series[j - 1] * e * sqrt(dt)
        # 画出每一个样本的数据序列
        plt.plot(x, random_series)

    plt.grid(ls='--')
    plt.show()

从上图中，我们可以看出虽然每一支股票的价格变化都不同，但是可以看出完整的变化区间以及变化密度(大量的股票都在中心区段变化，少量的股票在两端发生变化)。这个变化过程体现了蒙特卡罗方法的精髓，让我们从宏观上把握了整个随机过程。

两类重要的随机过程

• 到达过程

关注某种“到达”的事件是否发生，比如在一个服务窗口前，顾客的到达时刻；十字路口车辆依次通过的时刻。相邻间隔时间相互独立，且随机过程无记忆性。这里要细化为两种过程，在到达时间为离散的情况下，随机过程为伯努利过程，相邻间隔时间服从几何分布；另一类是到达时间为连续的情况下，随机过程为泊松过程，相邻间隔时间服从指数分布。

比如窗口上一时刻是否有顾客，都不对下一时刻是否来顾客提供信息；以及之前提到的博彩的信息一样，这一轮赌徒是赢还是输，不对下一轮赌局的输赢带来任何影响。

• 马尔科夫过程

未来的数据与历史数据有关联、有联系。最为重要的是马尔科夫模型，未来的数据只依赖于当前数据，而与过去的历史数据无关。

马尔科夫链的三要素

马尔科夫过程是随着时间而发生状态变换的过程。它可分为离散时间的马尔科夫链和连续时间的马尔科夫链两类，我们首先考虑离散时间的马尔科夫链。

离散时间的马尔科夫链在确定的离散时间点上发生变化。它一共有三要素：

1. 离散时间
2. 状态空间
3. 转移概率
4. 离散时间马尔科夫链的状态空间

我们通常用n表示时刻，(X_n)表示马尔科夫链此时的状态。所有的状态构成一个集合S={1,2,3,...,m}，该集合我们称为离散时间马尔科夫链的状态空间。

• 离散时间马尔科夫链的转移概率

转移概率(p_{i,j})：当前的状态是i，下一个状态等于j的概率。它本质上是条件概率：

(p_{i,j})=P((X_{n+1}=j|X_n=i)),i,j∈S

• 转移概率图举例：

在上图中，张三拥有三种状态——宅家、运动、吃美食，那么

状态空间={宅家、运动、吃美食}

如果今天宅在家中，明天宅家的概率为0.2；而如果今天宅在家中，明天出去吃美食的概率为0.6；如果今天宅在家中，明天出去运动的概率是0.2。

马尔科夫链的基本性质

只要n时刻的马尔科夫链状态为i，无论过去发生什么，无论马尔科夫链如何到达状态i，下一时刻转移到状态j的概率一定是转移概率(p_{i,j})。比如张三今天吃美食，不管他昨天是运动还是宅在家中，他明天出去运动的概率都是0.3。下一个状态只与当前状态有关，而与更早的历史状态无关。

对任意时刻n，状态空间中的任意状态i,j∈S，以及时刻n前任意可能的状态序列(i_0,i_1,i_2,...,i_{n-1})，均有：

P((X_{n+1}=j|X_n=i,X_{n-1}=i_{n-1},...,X_0=i_0))=P((X_{n+1}=j|X_n=i))=(p_{i,j})

• 转移概率的性质
• 首先(p_{i,j})满足非负性
• 其次(\sum_{j=1}^mp_{ij}=1)对所有的i都成立
• 当i=j时，(p_{i,j})状态发生了一次特殊的转移，即自身转移
• 状态转移矩阵

状态空间S中，任意两个状态i和j之间都有一个转移概率(p_{i,j})，并且满足(p_{i,j})≥0；所有状态间转移概率按照顺序组织成一个二维矩阵，其中第i行第j列的元素就是(p_{i,j})

上图就是状态转移矩阵，它刻画了对应的马尔科夫链的本质特征。

我们依然使用张三的例子来解读状态转移矩阵

我们令状态1为宅在家中，状态2为运动，状态3为吃美食。那么这个马尔科夫链就是一个3*3的矩阵，这里我们可以看到(p_{33})=0.6，它表示今天吃美食，明天也吃美食；(p_{32})=0.3，它表示今天吃美食，明天去运动。

多步转移概率的计算

之前我们知道，从状态i一步到达状态j的转移概率为：(p_{i,j})=P((X_{n+1}=j|X_n=i))，而m步的状态转移概率为：

(p^m(i,j)=P(X_{n+m}=j|X_n=i))

例：社会阶层流动的马尔科夫链，状态1是贫困，状态2是中产，状态3是财富自由，该链的状态转移矩阵为

假设爷爷处于贫困(状态1)，那么父亲处于中产(状态2)，而孙子处于财富自由(状态3)的概率有多大

该问题用条件概率表示就为P((X_2=3,X_1=2|X_0=1))

使用联合概率公式           =(P(X_2=3,X_1=2,X_0=1)\over P(X_0=1))           具体可以参考概率论整理(二) 中的二维离散型随机变量条件分布

                                         =(P(X_2=3,X_1=2,X_0=1)\over P(X_1=2,X_0=1))⋅(P(X_1=2,X_0=1)\over P(X_0=1))      此处为分子分母同乘(P(X_1=2,X_0=1))，保持不变

使用联合概率公式          =(P(X_2=3|X_1=2,X_0=1)⋅P(X_1=2|X_0=1))

根据马尔科夫链性质      =(P(X_2=3|X_1=2)⋅P(X_1=2|X_0=1))

                                        =(p_{23}⋅p_{12})=0.2*0.2=0.04

例：假设爷爷是贫困水平(状态1)，孙子处于财富自由(状态3)的概率有多大？

根据描述，该问题用条件概率表示就为                                             P((X_2=3|X_0=1))

父亲这一代可以处于任意状态(贫困、中产、财富自由)                 =(P(X_2=3,X_1=1|X_0=1)+P(X_2=3,X_1=2|X_0=1)+P(X_2=3,X_1=3|X_0=1))

根据上面的例子推导可知                                                                 =(p_{13}p_{11}+p_{23}p_{12}+p_{33}p_{13})

                                                                                                           =(\sum_{k=1}^3p_{k3}p_{1k})=0.1*0.7+0.2*0.2+0.4*0.1=0.15

• 多步转移与矩阵乘法

转移概率(P(X_2=3|X_0=1)=p_{13}p_{11}+p_{23}p_{12}+p_{33}p_{13})是转移矩阵

第一行和第三列点乘的结果。正好对应了从状态1到状态3两步状态转移的概率值。

概率转移矩阵的二次幂

，得到的结果矩阵就包含所有状态间通过两步到达的概率值。

import numpy as np

if __name__ == '__main__':

    A = np.array([[0.7, 0.2, 0.1], [0.3, 0.5, 0.2], [0.2, 0.4, 0.4]])
    print(np.dot(A, A))

运行结果

[[0.57 0.28 0.15]
 [0.4  0.39 0.21]
 [0.34 0.4  0.26]]

从结果的第一行，第三列，我们可以看到它刚好等于之前第二例的计算结果为0.15。

这里是两步的状态转移概率，如果是n步的话，其实就是n次幂。

import numpy as np

if __name__ == '__main__':

    A = np.array([[0.7, 0.2, 0.1], [0.3, 0.5, 0.2], [0.2, 0.4, 0.4]])

    def get_matrix_pow(n):
        ret = A
        for i in range(n):
            ret = np.dot(ret, A)
        print(ret)

    get_matrix_pow(3)
    get_matrix_pow(5)
    get_matrix_pow(10)
    get_matrix_pow(20)
    get_matrix_pow(100)
    get_matrix_pow(1000)

运行结果

[[0.4879 0.3288 0.1833]
 [0.4554 0.3481 0.1965]
 [0.4422 0.3552 0.2026]]
[[0.471945 0.338164 0.189891]
 [0.465628 0.341871 0.192501]
 [0.463018 0.343384 0.193598]]
[[0.46814979 0.34038764 0.19146257]
 [0.46804396 0.34044963 0.1915064 ]
 [0.46800013 0.34047531 0.19152456]]
[[0.46808512 0.34042552 0.19148935]
 [0.46808509 0.34042554 0.19148937]
 [0.46808508 0.34042555 0.19148937]]
[[0.46808511 0.34042553 0.19148936]
 [0.46808511 0.34042553 0.19148936]
 [0.46808511 0.34042553 0.19148936]]
[[0.46808511 0.34042553 0.19148936]
 [0.46808511 0.34042553 0.19148936]
 [0.46808511 0.34042553 0.19148936]]

通过以上的结果，我们会发现，随着n值的增大，它的结果会逐渐收敛于一个固定的概率矩阵。

状态转移矩阵的n次幂：n步状态转移概率收敛于

[[0.46808511 0.34042553 0.19148936]
 [0.46808511 0.34042553 0.19148936]
 [0.46808511 0.34042553 0.19148936]]

它的意义在于，无论你当前处于哪个阶层，n代之后(相对较大值)，变成各个阶层的概率是一样的，变成贫穷的概率是0.46，变成中产的概率是0.34，变成财富自由的概率是0.19。

路径概率问题

任何一个给定状态序列的概率：(P(X_0=i_0,X_1=i_1,...,X_n=i_n)=P(X_0=i_0)p_{i_0i_1}p_{i_1i_2}...p_{i_{n-1}i_n})

证明：(P(X_0=i_0,X_1=i_1,...,X_n=i_n))是一个联合概率，我们把(X_0=i_0,X_1=i_1,...,X_{n-1}=i_{n-1})视为条件A，(X_n=i_n)视为结果B，根据条件概率乘法公式(P (AB)=P (A) P (B|A))，有

(P(X_0=i_0,X_1=i_1,...,X_n=i_n))=(P(X_0=i_0,X_1=i_1,...,X_{n-1}=i_{n-1}))(P(X_n=i_n|X_0=i_0,X_1=i_1,...,X_{n-1}=i_{n-1}))

由于这是一个马尔科夫问题，则(P(X_n=i_n|X_0=i_0,X_1=i_1,...,X_{n-1}=i_{n-1}))=(P(X_n=i_n|X_{n-1}=i_{n-1}))=(p_{i_{n-1}}p_{i_n})，这就是一个转移概率，从状态(i_{n-1})转移到状态(i_n)，则上式可以得出

(P(X_0=i_0,X_1=i_1,...,X_n=i_n))=(P(X_0=i_0,X_1=i_1,...,X_{n-1}=i_{n-1}))(p_{i_{n-1}}p_{i_n})

进一步得到                                         =(P(X_0=i_0,...,X_{n-2}=i_{n-2})p_{i_{n-2}i_{n-1}}p_{i_{n-1}}i_n)

最终可得                                            =(P(X_0=i_0)p_{i_0i_1}p_{i_1i_2}...p_{i_{n-1}i_n})

例：从太爷爷开始，太爷爷是贫穷，爷爷是贫穷，爸爸是中产，儿子是财富自由，孙子中产，重孙子贫穷。

社会阶层流动的马尔科夫链：状态1是贫困，状态2是中产，状态3是财富自由，状态转移矩阵依然为

路径概率是：P((X_0=1))(p_{11}p_{12}p_{23}p_{32}p_{21})=1*0.7*0.2*0.2*0.4*0.3=0.00336，P((X_0=1))表示太爷爷是贫穷的概率，定为1.

马尔科夫过程的两种典型极限状态

• 极限与初始状态无关的情况

社会流动概率转移矩阵

随着n增大，收敛于

。当n->∞时，矩阵中的每一个值都会收敛于一个极限值，这个极限值不依赖于初始状态。这也就是说你的老祖宗无论处于哪一个阶层，过了很多代后，到你这一代处于任何一个阶层的概率都是一定的，跟你的老祖宗没有关系了。

当时间n比较小的时候，n步转移概率矩阵中的值还会比较依赖于他的初始状态i，而当n不断增大时，这种依赖性将会逐渐消失。

• 极限依赖于初始状态的情况

随着n->∞时，矩阵中的每一个值会收敛于一个极限值，但这些极限值会依赖于最初位于哪一个初始状态。具体表现在矩阵中不是每一行都相等。

上图表示一只绵羊在四个区域中移动，每次只能移动一个位置，如果在区域2和3之间移动则平安无事，如果到区域1或4，则立即被老虎吃掉。它的状态转移概率矩阵如下

import numpy as np

if __name__ == '__main__':

    # 状态转移概率矩阵
    A = np.array([[1, 0, 0, 0], [0.2, 0.4, 0.4, 0], [0, 0.4, 0.4, 0.2], [0, 0, 0, 1]])

    def get_matrix_pow(n):
        ret = A
        for i in range(n):
            ret = np.dot(ret, A)
        print(ret)

    get_matrix_pow(5)
    get_matrix_pow(10)
    get_matrix_pow(50)
    get_matrix_pow(100)
    get_matrix_pow(1000)

运行结果

[[1.       0.       0.       0.      ]
 [0.468928 0.131072 0.131072 0.268928]
 [0.268928 0.131072 0.131072 0.468928]
 [0.       0.       0.       1.      ]]
[[1.         0.         0.         0.        ]
 [0.55705033 0.04294967 0.04294967 0.35705033]
 [0.35705033 0.04294967 0.04294967 0.55705033]
 [0.         0.         0.         1.        ]]
[[1.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00]
 [5.99994291e-01 5.70899077e-06 5.70899077e-06 3.99994291e-01]
 [3.99994291e-01 5.70899077e-06 5.70899077e-06 5.99994291e-01]
 [0.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00]]
[[1.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00]
 [6.00000000e-01 8.14814391e-11 8.14814391e-11 4.00000000e-01]
 [4.00000000e-01 8.14814391e-11 8.14814391e-11 6.00000000e-01]
 [0.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00]]
[[1.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00]
 [6.00000000e-01 4.92092769e-98 4.92092769e-98 4.00000000e-01]
 [4.00000000e-01 4.92092769e-98 4.92092769e-98 6.00000000e-01]
 [0.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00]]

根据结果我们可以看出，随着n的不断增大，转移概率收敛于

这里状态1和4被称为吸收态，一旦到达吸收态，就没有机会转移到其他状态；

n步转移概率矩阵收敛的极限值依赖于初始状态：最开始羊位于区域2，它死于区域1的概率是0.6；最开始位于区域3，死于区域1的概率是0.4；n步转移概率矩阵中，状态2和3两列都收敛于0：无论羊初始位于哪个状态，最终都死于区域1或4，存活的概率为0。

• 两个马尔科夫过程的比较

通过观察，我们可以看到在阶层流动性中

任取两个状态,i和j，从状态i出发都可以到达状态j。

从羊入虎口图中，从状态2出发可以到达任意状态，从状态1出发，到达不了除自身外的其他状态。

马尔科夫链中的常返类和周期性

• 可达与常返

对于状态i，如果对于每个从状态i出发可达的状态j，相应的从状态j出发反过来也可到达状态i，那么状态i就是常返的。

数学语言描述就是，令状态i的可达状态集合为A(i)，对于集合中的每一个状态j，如果i∈A(j)，那么状态i就是常返的。

在阶层流动的状态变化中，所有的状态都是常返的；羊入虎口中状态2、3是非常返，状态1、4是常返。

• 常返类

在上图中，我们可以看到，状态1和状态2是常返状态，状态3是非常返状态。状态1的可达状态集是{状态1，状态2}，因此状态1和状态2构成了一个常返类。

于是该图就可以分为1个常返类(包含两个状态)和1个非常返状态(类)。上图中虚线部分就是常返类，状态3则是非常返类。

• 常返类的重要特性
• 常返类只进不出。
• 不管开局如何，终将进入常返类。
• 有多个常返类的马尔科夫链，一定不会收敛于一个唯一的稳态分布。
  1. ​

1.  在羊入虎口的马尔科夫链中，就有两个常返类。单一常返类如果是周期的，也不收敛。    ​

1. 在上图中就是一个周期性的，n->∞，如果n为奇数，则从状态1转移到状态2；如果n为偶数，则从状态2转移到状态1。它是不收敛的。
2. 对于常返类R：只要存在n≥1和R中特定的状态i，使得经过n步之后可以到达R中的任意状态，常返态就是非周期的。在上图中并不存在这样一个n值它既可以到达状态1，也可以到达状态2.
3. ​

1. 在上图中，n=3时，从状态1可以到达常返类中所有状态，非周期。比如1-2-2-2，1-2-3-3，1-2-3-4，1-2-2-1。

马尔科夫链的稳态分析和求法

稳态的概念源于：当n->∞时，n步转移概率矩阵中每个数值的收敛情况。

在稳态下：对于每一个状态j，n步转移概率会趋近于一个独立于初始状态i的极限值，记作(π_j)。在稳态下处于状态j的概率也同样趋近于极限值(π_j)，即(P(X_n=j)->π_j)(当n->∞)。随着n->∞，马尔科夫链要收敛于一个稳态分布，必须是非周期的；如果要求稳态分布是唯一的，马尔科夫链必须只含有一个常返类。

• 稳态的求法

还是以之前的阶层流动的马尔可夫链为例，令贫穷阶层为状态1，中产阶层为状态2，财务自由为状态3，马尔科夫链到达稳态后，三个概率趋近于极限值：(π_1,π_2,π_3)，由分布概率的归一性原则，有：(π_1+π_2+π_3=1)，稳态的本质：到达稳态后，经过下一个时间的状态转移，马尔科夫链的概率分布保持不变。

1. 针对状态1：(0.7π_1+0.3π_2+0.2π_3=π_1)
2. 针对状态2：(0.2π_1+0.5π_2+0.4π_3=π_2)
3. 针对状态3：(0.1π_1+0.2π_2+0.4π_3=π_3)

当然可以使用矩阵的形式

包含归一性(π_1+π_2+π_3=1)，可得

1. (π_1=0.46808511)
2. (π_2=0.34042553)
3. (π_3=0.19148936)

方程组的解和n步转移概率矩阵各行的值

完全一致。

隐马尔科夫模型

隐马尔科夫模型(HMM)是一种统计模型，广泛应用于语音识别，词性自动标注等自然语言处理领域。

首先由一个隐藏的马尔科夫链生成一个状态随机序列，再由状态随机序列中的每一个状态对应生成各自的观测，由这些观测组成一个观测随机序列。因此隐马尔科夫模型伴随着两条线：一个是观测随机序列这条“明线”，另一个是隐藏着的状态随机序列这条“暗线”。

典型实例1：盒子摸球实验

有3个盒子，编号分别是1、2、3号，每个盒子里都装着个数不等的黑球和白球。

1. 1号盒子：黑球2个，白球8个
2. 2号盒子：黑球6个，白球4个
3. 3号盒子：黑球4个，白球6个

每次先随机出现一个盒子，然后我们从随机出现的盒子中随机摸出一个球，并且记录下球的颜色，最后把球放回盒子，然后循环往复。

• 实验中的两个序列

每次只能看到被摸出的球的颜色，无法知道每次随机出现的盒子编号，这是一个大的前提。

1. 试验中会依次出现不同编号的盒子，由于无法观测到盒子的编号，因此盒子的序列被称为隐含状态序列，是暗线。
2. 能够观察到的是球的颜色，因此球的颜色序列就是观测序列，是明线。

例如，试验重复7次，其中一种可能的观测序列为：O={黑，黑，白，白，白，黑，黑}；假定每次盒子随机出现的过程是一个马尔科夫过程，状态集合为Q={盒子1，盒子2，盒子3}，令第一次各个盒子出现的概率为：1号盒子0.3，2号盒子0.5，3号盒子0.2，即π=((0.3,0.5,0.2)^T)，状态转移图如下所示

则各个盒子之间相互转移的概率矩阵为：

• 试验中观测序列的生成
• 1号盒子：黑球2个，白球8个，放回式摸球，这个是一个古典概型，摸出黑球的概率为0.2，摸出白球的概率为0.8：是从特定隐含状态中生成指定规则的观测概率。
• 同样把2号盒子和3号盒子的观测概率集中在一个矩阵中，得到观测概率矩阵

，

观测集合V={黑球，白球}

重复7次上述过程，得到两个序列：

1. 长度为7的隐藏状态序列：I={2号盒，2号盒，1号盒，3号盒，1号盒，2号盒，3号盒}(这是无法观测到的)
2. 长度为7的观测序列：O={黑球，黑球，白球，黑球，白球，白球，黑球}

状态序列和观测序列

状态转移和观测输出概率

贝叶斯统计推断：最大后验

贝叶斯定理

P((A_i|B))=(P(B|A_i)P(A_i)\over P(B))