python贝叶斯随机过程:马尔可夫链Markov-Chain,MC和Metropolis-Hastings,MH采样算法可视化

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

由Kaizong Ye,Liao Bao撰写

马尔可夫链的不变测度

from math import pi
from pylab import *

考虑一个高斯 AR(1) 过程, , 其中  是标准高斯随机变量的独立同分布序列,独立于 。假使 .。然后,具有均值的高斯分布  和方差  是马尔可夫链的平稳分布。我们用马尔可夫链的单个轨迹所取值的直方图来检查这个属性。

f=lambda x,m,sq: np.exp(-(x-m)\*\*2/(2\*sq))/np.sqrt(2\*pi\*sq)


plt.hist

第二个例子

我们在这里考虑一个马尔可夫链的例子,它的状态空间  是开单位区间。如果链条在 ,它等概率  选择两个区间之一  或者  ,然后移动到一个点,  它均匀分布在选定的区间内。马尔可夫链的不变分布有 cdf, 。 通过微分,我们可以得到相关的密度:  。对所有 , 我们现在用马尔可夫链取值的直方图检查这个属性。


视频

贝叶斯推断线性回归与R语言预测工人工资数据

探索见解

去bilibili观看

探索更多视频

x=arange(1,m)/m

for i in range(p-1):
    \[a,b\]=rand(2)
    
plt.hist

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

anm = animation.FuncAnimation

R语言中实现马尔可夫链蒙特卡罗MCMC模型

阅读文章


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

data = \[\]
for i in range(p-1):
    \[a,b\]=npr.rand(2
    
    if ((i+1)%100==0):
        data.append
    
anim = animation.Func


随时关注您喜欢的主题


我们现在用一个例子来说明大数定律。如 。 那么,我们期望 ,

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

对称随机游走 Metropolis Hasting 算法

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

为了针对此分布,我们根据对称随机游走 Metropolis Hasting 算法进行采样。

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

当链条处于状态时 ,我们提出一个候选 , 根据  ,其中  。然后我们接受  ,有概率 , 其中 . 否则, .

from IPython.display import HTML


rc('animation', html='jshtml')
ani

独立Metropolis Hasting 算法

我们再次考虑一个目标分布,它是两个高斯分布的混合,一个集中在  ,另一个集中在 ,其中  是中心标准正态分布的密度。

为了针对这种分布,我们根据具有独立提议的 Metropolis Hasting 算法进行采样。当链条处于状态时 ,我们提出一个候选  ,根据  ,其中  。然后我们接受  有概率 , 其中  和  是密度 .。否则, .。

mc=npr.randn*np.one
data=\[\]
for i in range:
    v=sig*npr+sft
    alpha
    if (npr.rand()<alpha): 
        mc\[i+1\] = v        
    if ((i+1)%r==0):
        data.append

x=np.linspac


    
anim = animation.FuncAn


可下载资源

关于作者

Kaizong Ye拓端研究室(TRL)的研究员。

本文借鉴了作者最近为《R语言数据分析挖掘必知必会 》课堂做的准备。

​非常感谢您阅读本文,如需帮助请联系我们!

 
QQ在线咨询
售前咨询热线
15121130882
售后咨询热线
0571-63341498

关注有关新文章的微信公众号


永远不要错过任何见解。当新文章发表时,我们会通过微信公众号向您推送。

技术干货

最新洞察

This will close in 0 seconds