原文：/articles/66987#2
1什么是MCMC，什么时候使用它
MCMC只是一种从分布中抽样的算法。这个术语代表“马尔科夫蒙特卡罗”，因为它是一种使用了 “马尔科夫链”的“蒙特卡罗”（即随机）方法。MCMC只是一种蒙特卡罗方法。
（文中有些使用R语言实现的没有改动）
2为什么我要从分布抽样呢？
从分布抽烟是解决一些问题的最简单的方法
也许在贝叶斯推断中最常见的方式是使用MCMC从某些模型的后验概率分布中抽取样本。有了这些示例，你就可以问这样的问题：什么是参数的平均值和可信区间？”。
例如，加入你有合适的参数模型的后验概率密度是某个函数f(x,y)。然后，计算参数x的平均值，你可以这样计算
你可以简单的读作“x乘参数（x,y）的概率，积分在x和y可能采用的所有可能值上。”
另一种计算这个值得方法是模拟观察值K，
从f(x,y)计算样本平均值为
其中x(j）为第j个样本的x的值。如果这些样本是服从独立分布的样本，那么随着K趋于无穷大，x的估计平均值会收敛于真平均数。
假设我们的目标分布是一个正态分布，平均数为m,标准偏差为s。显然，这种分布的平均值是m，让我们试着通过从分布中抽取样本来显示。
例1估计一个平均值为m，标准差为s的正态分布的平均值。
import numpy as np
import statsmodels.formula.api as sm
import matplotlib.pyplot as plt
if __name__ == "__main__":
mu = 0
sigma = 1
samplenum = 10000
#设置随机数生成器种子，每次运行的随机数都是一样的。
np.random.seed(0)
samples = np.random.normal(mu,sigma,samplenum)
mu_sample = np.mean(samples)
#计算累加平均数
x = list(samples.cumsum(0))
n = [i for i in range(1,len(x)+1)]
cummean = [a/b for a,b in zip(x,n)]
马尔可夫链蒙特卡罗
定义
设
表示在t时间一些随机变量的值。马尔可夫链起始于
产生一系列样点
，接着一系列随机步骤。
马尔可夫链满足马尔可夫属性。马尔可夫属性是已知现在状态的条件下，将来所处的状态与过去状态无关，即“遗忘性”；基本上不管你如何达到某个状态x，x的转移概率不变：
从一个步骤到下一个步骤的转换是由转移核描述，它可以由状态i到状态j的转变的概率（或连续变量的概率密度）描述为
表示状态j在t时间（步骤）的链的概率，并且定义
为可能的状态的概率的向量。然后，给出
，我们利用Chapman-Kolmogorov方程计算
。
这个概率是我们在状态k乘以从k到i的转变的概率，在所有可能的状态k上求和。设P是概率转移矩阵——矩阵的第i,j元素是P(i→j)，并且将上述等式重写为
我们可以轻松地迭代这个方程：
平稳分布
如果有一些矢量
满足
那么$\vec\pi^$是这个马尔可夫链的平稳分布。直观地，认为这个系统将设置的状态的最终特征集；运行足够长的时间，该系统会“遗忘”它的初始状态，则这个向量的第i个元素是这个系统将在状态i的概率。
如果这个过程是不可约的和非周期性的，马尔可夫链将会有一个平稳分布。
在数学上，
是特征值为1的左特征向量。
这里有一个快速的定义，使事情更具体（但请注意，这无关MCMC本身–这只是对马尔可夫链的观点！）假设我们有一个三态马尔可夫过程。设P为链的概率转移矩阵：
p = np.array([[0.5,0.25,0.25],[0.2, 0.1, 0.7],[0.25, 0.25, 0.5]])
print(p[1,1])
请注意，P行的共计为一：
入口P[i,j]给出了从状态i到状态j的概率（所以这就是上述的P(i→j)）。
注意，与行不同，列不一定共计为1：
这个函数取一个状态向量x（其中x[i]是处于状态i的概率），并通过乘以转移矩阵P来迭代它，使系统前进n个步骤。
def iterate(x,p,n):
res = np.zeros((n+1,x.size))
res[0,] = x
for i in range(n):
res[i+1,] = x = np.dot(x,p)
return res从状态1中的系统开始（所以x是向量[ 1,0,0 ]，指示存在处于状态1的100％概率，并且没有处于任何其他状态的机会），并且迭代10个步骤 ：
n = 10
x = np.array([1,0,0])
y1 = iterate(np.array([0,1,0]),p,n)同样，对于其他两个可能的起始状态：
y2 = iterate(np.array([1,0,0]),p,n)
y3 = iterate(np.array([0,0,1]),p,n)这表明在平稳分布上的收敛性。
X=np.linspace(0,n+1,11)
plt.plot(X,y1,"-",X, y2, ":",X, y3, "-.")
plt.legend(loc="best",ncol=3,labels=["line11","line12","line13","line21","line22","line23","line31","line32","line33"])
plt.show()
这意味着不管起始分布如何，不管它在哪里开始，在大约10次或更多次迭代之后，链处于状态1的概率为32％。
我们可以用R的特征函数来提取系统的主导特征向量（这里的t()转置矩阵得到左特征向量）。
然后添加点到之前显示我们接近收敛的图中：
根据特征向量的定义，乘以特征向量的转移矩阵返回特征向量本身：
（严格地说，这应该是由V相乘的特征值，但这些矩阵的主特征值总是1）。
在这里运行的函数用一个状态（这一次，只是一个整数代表其中系统的状态1,2,3），如上相同的转移矩阵，和一些运行的步骤。每个步骤，查看可能转换到的地方，并选择1（这使用R的样本函数）。
这是100步左右的链：
绘制我们在每个状态随着时间的推移的时间分数。
运行这步再长一点（5000步）
存在稳定分布的充分（但不是必要的）条件是细致平衡，其为：
这意味着链是可逆的。这种情况意味着一个平稳分布存在的原因是
将状态j的细致平衡方程的两侧求和
左边的项等于
的第k个元素，右边的项可以作为因子
然后，因为
因为P是一个转移概率函数，由总概率定律决定，概率为1）所以右边是
，所以我们有
其适用于所有k，所以
马尔可夫链具有平稳分布，如果我们运行它们足够长的时间，我们可以看看链在哪里花费时间，并得到该平稳分布的合理估计。
