初看MCMC采样的时候，一脑子的雾水。不尽会问，马尔科夫链稳定后怎么采样啊，都稳定了每次采样出来的东西不是一样的吗？这时候特别希望有个例子能看看，怎么采样的，采出来的是什么东西。很不幸，网上的例子全是连续函数的版本的，没有离散版本的。那么今天，我就把我的代码分享给大家。
'''
Created on 2018年5月16日
p:输入的概率分布，离散情况采用元素为概率值的数组表示
N:认为迭代N次马尔可夫链收敛
Nlmax:马尔可夫链收敛后又取的服从p分布的样本数
isMH:是否采用MH算法，默认为True
1）输入我们任意选定的马尔科夫链状态转移矩阵Q，平稳分布π(x)，设定状态转移次数阈值n1，需要的样本个数n2
2）从任意简单概率分布采样得到初始状态值x0
3）for t=0 to n1+n2−1:
a) 从条件概率分布Q(x|xt)中采样得到样本x∗
b) 从均匀分布采样u∼uniform[0,1]
c) 如果u<α(xt,x∗)=π(x∗)Q(x∗,xt), 则接受转移xt→x∗，即xt+1=x∗
d) 否则不接受转移，即xt+1=xt
样本集(xn1,xn1+1,...,xn1+n2−1)即为我们需要的平稳分布对应的样本集。
'''
from __future__ import division
import matplotlib.pyplot as plt
import numpy as np
from array import array
def mcmc(Pi ,Q,N=1000,Nlmax=100000,isMH=False):
X0 = np.random.randint(len(Pi))# 第一步：从均匀分布（随便什么分布都可以）采样得到初始状态值x0
T = N+Nlmax-1
result = [0 for i in range(T)]
t = 0
while t < T-1:
t = t + 1
# 从条件概率分布Q(x|xt)中采样得到样本x∗
# 该步骤是模拟采样，根据多项分布，模拟走到了下一个状态
#（也可以将该步转换成一个按多项分布比例的均匀分布来采样）
x_cur = np.argmax(np.random.multinomial(1,Q[result[t-1]])) # 第二步：取下一个状态 ，采样候选样本
if isMH:
'''
细致平稳条件公式：πi Pij=πj Pji,∀i,j
'''
a = (Pi[x_cur] * Q[x_cur][result[t-1]]) /(Pi[result[t-1]] * Q[result[t-1]][x_cur]) # 第三步：计算接受率
acc = min(a ,1)
else: #mcmc
acc = Pi[x_cur] * Q[x_cur][result[t-1]]
u = np.random.uniform(0 ,1) # 第四步：生成阈值
if u< acc: # 第五步：是否接受样本
result[t]=x_cur
else:
result[t]= result[t-1]
return result
def count(q, n):
L = array("d")
l1 = array("d")
l2 = array("d")
for e in q:
L.append(e)
for e in range(n):
l1.append(L.count(e))
for e in l1:
l2.append(e / sum(l1))
return l1, l2
if __name__ == '__main__':
Pi = np.array([0.5, 0.2, 0.3]) # 目标的概率分布
#状态转移矩阵，但是不满足在 平衡状态时和 Pi相符
#我们的目标是按照某种条件改造Q ，使其在平衡状态时和Pi相符
#改造方法就是，构造矩阵 P，且 P(i,j)=Q(i,j)α(i,j)
# α(i, j) = π(j)Q(j, i)
# α(j, i) = π(i)Q(i, j)
Q = np.array([[0.9, 0.075, 0.025],
[0.15, 0.8, 0.05],
[0.25, 0.25, 0.5]])
a = mcmc(Pi,Q)
l1 = ['state%d' % x for x in range(len(Pi))]
plt.pie(count(a, len(Pi))[0], labels=l1, labeldistance=0.3, autopct='%1.2f%%')
plt.title("markov:" +str(Pi))
plt.show()
可以看见。我们的转移概率矩阵Q的平衡状态是[0.61.0.31.0.8],和我们的目标分布 Pi [0.5, 0.2, 0.3]是不相同的，所以我们的目标是构造一个转移概率矩阵 P， P(i,j)=Q(i,j)α(i,j)，使得P的平衡状态满足PI。代码的结果为：
由于mcmc采样效率低（接受率低），需要多采一点样本才看起来更像。