1. 求离散马尔科夫链的平稳分布的两种解法
假设离散马尔科夫链的转移矩阵为P PP,平稳分布为π \piπ,则平稳分布满足:
P π = π P \pi = \pi Pπ=π
1.1 迭代法
求平稳分布的一种简单方法是迭代法,即随机初始化初始分布π 0 \pi_0π
0
,利用上式不断迭代求解下一时刻的状态分布直到状态分布收敛,则求得平稳分布。
"""通过markov_marix迭代得到`平稳分布`. """
pi0 = zeros(num)
pi0[0] = 1
pi = list([pi0])
for i in range(200):
print('epoch %02d' % (i+1), pi[-1])
pi.append(markov_marix.dot(pi[-1]))
if sum(abs(pi[-1]-pi[i])) <= 1e-7:
print('break...', pi[-1])
break
steady_vec = pi[-1] # 平稳分布
print('平稳分布:', steady_vec)
1.2 特征分解法
另一种解法是利用特征分解,由于平稳分布满足P π = π P\pi=\piPπ=π,与特征方程A x = λ x Ax=λxAx=λx联系可知平稳分布π ππ就是转移矩阵P PP的特征值为1对应的特征向量(归一化),因此可以直接对转移矩阵进行特征分解来求平稳分布。
"""特征分解求平稳分布. """
values, vecs = eig(markov_marix,)
for i in range(len(values)):
if abs(values[i]-1) <= 1e-9:
print(values[i])
print(vecs[:, i]/sum(vecs[:, i]))
模拟结果如下,两种方法求得平稳分布一致
2. MCMC方法: 按照一定概率分布生成想要的样本
按照一定的概率分布生成我们想要的样本,实在实际应用中非常重要的方法, 相关的方法也有许多了, 对于给定的一维函数 f(x) , 有许多常用的采样方法是的样本符合 f(x) 的概率分布, 但是想要模拟一个其分量为相关随机变量的随机向量 X ,确实一件困难的事. 下面我们介绍一个强有力的方法来生成分布近似为 X 的分布的随机向量. 这种方法被称为 MCMC 方法.
2.1 马氏链
考虑一组随机变量: X0,X1,… . 这里把 Xn 解释为"在时间点n的系统的状态, 并假设 Xn 可能取值的集合,即系统所能达到的状态的集合为 1,…,N . 如果存在一组数
P
i
j
,
i
,
j
=
1
,
.
.
.
,
N
P_{ij},i,j=1,...,N
Pij,i,j=1,...,N , 使得过程无论什么时候处于状态 i , 不用管前面的状态如何, 其下一个状态是 j 的概率为
P
i
j
P_{ij}
Pij , 我们便称 {Xn,n≥0} , 构成一个转移矩阵为
P
i
j
,
i
,
j
=
1
,
.
.
.
,
N
P_{ij},i,j=1,...,N
Pij,i,j=1,...,N 的马氏链. 显然,从上一个时间点上的状态 i 转移到下一个时间点上的所有的状态的概率之和应该是1:
现在的主要问题是: 我们想要生成一个符合我们要求的稳定概率的马氏链, 那么我们通过运行这个马氏链来获得样本,这个样本就应该符合我们想要的概率.
2.2 Hastings-Metropolos算法
在介绍这个算法之前, 我们先来给出一个结论: 定理:[细致平稳条件] 如果非周期马氏链的转移矩阵 P 和分布 π(x) 满足:
π ( i ) P i j = π ( j ) P i j , f o r a l l i , j ( 1 ) π(i)P_{ij}=π(j)P_{ij},\ for\ all\ i,j\ (1) π(i)Pij=π(j)Pij, for all i,j (1)
则 π ( x ) \pi(x) π(x) 是马氏链的平稳分布,上式被称为细致平稳条件(detailed balance condition)。
我们可以通过这个条件, 来构建任意我们想要的稳定概率的转移矩阵.
原文链接:https://blog.youkuaiyun.com/weixin_43486780/article/details/104407875
http://logicgogh.github.io/articles/MCMC