采样-MCMC

材料

https://www.bilibili.com/video/BV17D4y1o7J2

1.随机抽样问题

$假设已知随机变量X的概率密度函数f(x),求X^2的期望?$
$可以通过公式E(g(X))=\int_{-\infty}{+\infty}f(x)g(x)dx,令g(X)=X^2即可$
$但是这个积分怎么求解?$
$不妨抽取样本X,然后计算g(X),进而算术平均$
$E(g(X))=\frac{1}{N}\sum_{i=1}^{N}g(x_i)$
$\color{red}{问题:如何抽取服从特定分布的样本呢?}$

2.马尔科夫链MC采样

$意义不细说了,举个例子$
$股票系统有三种状态牛市,熊市,横盘$
$状态转义矩阵$
$\begin{pmatrix} & 牛市 & 熊市 & 横盘\\ 牛市 & 0.9 & 0.075 & 0.025\\ 熊市 & 0.15 & 0.8 & 0.05\\ 横盘 & 0.25 & 0.25 & 0.5\\ \end{pmatrix}$
$那么此时当系统的初始状态经过N次变换后终究会变成稳态$
$相关马尔科夫链的验证代码如下$

点击查看代码

import numpy as np
import random


# 返回一个数组中某个元素数量
def count_element(source_array, element):
    count = 0
    for x in source_array:
        if x == element:
            count += 1
    return count

# 矩阵的左标识原状态,矩阵的上标识转移后的状态,矩阵的值是原状态转移到现状态的概率
transfer_matrix = np.array([[0.9, 0.075, 0.025], [0.15, 0.8, 0.05], [0.25, 0.25, 0.5]])
# transfer_m = np.array([[0.7, 0.2, 0.1], [0.2, 0.4, 0.4], [0.1, 0.4, 0.5]]) # 这个数据算出来的是1/3 ,1/3,1/3 这种均数
print(transfer_matrix)
# for i in range(100):
#     print(np.dot(transfer_m, transfer_m))
# first sample
# 随便选个初始值
start_m = np.array([[1, 0, 0]])
for i in range(20):
    print(np.dot(start_m, transfer_matrix))
    start_m = np.dot(start_m, transfer_matrix)

# 初始值
s = []
for i in range(100):
    s.append('C')

for i in range(100):
    # print(i, "迭代")
    print("样本", i, "中A,B,C的数量", count_element(s, 'A'), count_element(s, 'B'), count_element(s, 'C'))
    prob = None
    for i in range(100):
        if s[i] == 'A':
            prob = transfer_matrix[0]
        elif s[i] == 'B':
            prob = transfer_matrix[1]
        elif s[i] == 'C':
            prob = transfer_matrix[2]
        rnd = random.randint(0, 100) / 100
        if rnd < prob[0]:
            s[i] = 'A'
        elif rnd < prob[0] + prob[1]:
            s[i] = 'B'
        else:
            s[i] = 'C'

[0.25 0.25 0.5 ]]
[[0.9 0.075 0.025]]
[[0.8275 0.13375 0.03875]]
[[0.7745 0.17875 0.04675]]
[[0.73555 0.212775 0.051675]]
[[0.70683 0.238305 0.054865]]
[[0.685609 0.2573725 0.0570185]]
[[0.6699086 0.2715733 0.0585181]]
[[0.65828326 0.28213131 0.05958543]]
[[0.64967099 0.28997265 0.06035636]]
[[0.64328888 0.29579253 0.06091859]]
[[0.63855852 0.30011034 0.06133114]]
[[0.635052 0.30331295 0.06163505]]
[[0.63245251 0.30568802 0.06185947]]
[[0.63052533 0.30744922 0.06202545]]
[[0.62909654 0.30875514 0.06214832]]
[[0.62803724 0.30972343 0.06223933]]
[[0.62725186 0.31044137 0.06230677]]
[[0.62666957 0.31097368 0.06235675]]
[[0.62623785 0.31136835 0.0623938 ]]
[[0.62591777 0.31166097 0.06242126]]
样本 0 中A,B,C的数量 0 0 100
样本 1 中A,B,C的数量 23 20 57
样本 2 中A,B,C的数量 34 32 34
样本 3 中A,B,C的数量 41 34 25
样本 4 中A,B,C的数量 46 40 14
样本 5 中A,B,C的数量 54 36 10
样本 6 中A,B,C的数量 51 40 9
样本 7 中A,B,C的数量 48 47 5
样本 8 中A,B,C的数量 53 39 8
样本 9 中A,B,C的数量 57 36 7
样本 10 中A,B,C的数量 59 35 6
样本 11 中A,B,C的数量 57 34 9
样本 12 中A,B,C的数量 61 35 4
样本 13 中A,B,C的数量 62 34 4
样本 14 中A,B,C的数量 66 29 5
样本 15 中A,B,C的数量 64 27 9
样本 16 中A,B,C的数量 71 23 6
样本 17 中A,B,C的数量 62 29 9
样本 18 中A,B,C的数量 66 29 5
样本 19 中A,B,C的数量 70 28 2
样本 20 中A,B,C的数量 66 31 3
样本 21 中A,B,C的数量 65 31 4
样本 22 中A,B,C的数量 66 32 2
样本 23 中A,B,C的数量 67 28 5
样本 24 中A,B,C的数量 64 29 7
样本 25 中A,B,C的数量 62 31 7
样本 26 中A,B,C的数量 66 32 2
样本 27 中A,B,C的数量 63 34 3
样本 28 中A,B,C的数量 66 31 3
样本 29 中A,B,C的数量 65 28 7
样本 30 中A,B,C的数量 64 33 3

$经过运行可知,数据分布逐渐收敛,并且采集的样本基本也是满足最终收敛的分布的$
$\color{red}{这里有一个很有意思的问题}$
transfer_matrix = np.array([[0.9, 0.075, 0.025], [0.15, 0.8, 0.05], [0.25, 0.25, 0.5]])
transfer_m = np.array([[0.7, 0.2, 0.1], [0.2, 0.4, 0.4], [0.1, 0.4, 0.5]])
当我用下面的那个转移矩阵算出来的收敛分布是[1/3 ,1/3,1/3],有点出乎我的意料
这里留一个有意思的课题,根据状态转移矩阵怎么算出最终的收敛分布