拓端tecdat|python贝叶斯随机过程：马尔可夫链Markov-Chain，MC和Metropolis-Hastings，MH采样算法可视化

介绍

本文，我们说明了贝叶斯学习和 计算统计一些结果。

1.   
2.  from math import pi
3.  from pylab import *

马尔可夫链的不变测度

考虑一个高斯 AR(1) 过程， 

， 其中 

 是标准高斯随机变量的独立同分布序列，独立于 。假使 

.。然后，具有均值的高斯分布 

 和方差 

 是马尔可夫链的平稳分布。我们用马尔可夫链的单个轨迹所取值的直方图来检查这个属性。

1.   
2.  f=lambda x,m,sq: np.exp(-(x-m)**2/(2*sq))/np.sqrt(2*pi*sq)
3.   
4.   
5.  plt.hist

第二个例子

我们在这里考虑一个马尔可夫链的例子，它的状态空间 

 是开单位区间。如果链条在 ，它​等概率 

 ​选择两个区间之一 

 或者 

 ，然后移动到一个点，  它均匀分布在选定的区间内。马尔可夫链的不变分布有 cdf， 

。 通过微分，我们可以得到相关的密度： 

 。对所有 

， 我们现在用马尔可夫链取值的直方图检查这个属性。

1.   
2.  x=arange(1,m)/m
3.   
4.  for i in range(p-1):
5.  [a,b]=rand(2)
6.   
7.  plt.hist

我们还可以说明直方图如何收敛到平稳分布的密度。这可以通过使用 matplotlib 中的“动画”模块的动态动画来完成。下面是python代码。

1.   
2.   
3.  anm = animation.FuncAnimation

以这个例子结束，这是一个动画。

1.   
2.   
3.  data = []
4.  for i in range(p-1):
5.  [a,b]=npr.rand(2
6.   
7.  if ((i+1)%100==0):
8.  data.append
9.   
10.  anim = animation.Func

我们现在用一个例子来说明大数定律。如 

。 那么，我们期望 

,

1.   
2.  x=np.arange/(p)
3.  for i in range(p-1):
4.  [a,b]=npr.rand
5.  m=np.cumsum(g(m))/np.arange(1,p+1)
6.  plot

对称随机游走 Metropolis Hasting 算法

我们现在考虑一个目标分布，它是两个高斯分布的混合，一个集中在  ，另一个集中在 

。

  是中心标准正态分布的密度。

为了针对此分布，我们根据对称随机游走 Metropolis Hasting 算法进行采样。当链条处于状态时 ，我们提出一个候选 

， 根据 

 ，其中 

 。然后我们接受 

 ，有概率 

， 其中 

. 否则， 

.

1.  from IPython.display import HTML
2.   
3.   
4.  rc('animation', html='jshtml')
5.  ani

独立Metropolis Hasting 算法

我们再次考虑一个目标分布，它是两个高斯分布的混合，一个集中在  ，另一个集中在 ，

，其中  是中心标准正态分布的密度。为了针对这种分布，我们根据具有独立提议的 Metropolis Hasting 算法进行采样。当链条处于状态时 ，我们提出一个候选 

 ，根据 

 ，其中 

 。然后我们接受 

 有概率 

， 其中 

 和  是密度 

.。否则， 

.。

1.   
2.  mc=npr.randn*np.one
3.  data=[]
4.  for i in range:
5.  v=sig*npr+sft
6.  alpha
7.  if (npr.rand()<alpha):
8.  mc[i+1] = v
9.  if ((i+1)%r==0):
10.  data.append
11.   
12.  x=np.linspac
13.   
14.   
15.   
16.  anim = animation.FuncAn