MCMC算法—MH算法的Python实现（一）

针对离散数据概率分布的MCMC算法python实现

对于mcmc算法，如何选择状态转移矩阵对实验结果是否有影响，设计下面几组实验

输入为p=[0.9,0.05,0.05]，取q为[1/3,1/3,1/3],输出10000000个样本，观察样本的概率分布，代码如下：

#coding=utf-8
'''
Created on 2016年9月14日
基本MCMC算法以及M-H算法的python实现
@author: whz
p:输入的概率分布，离散情况采用元素为概率值的数组表示
N:认为迭代N次马尔可夫链收敛
Nlmax:马尔可夫链收敛后又取的服从p分布的样本数
isMH:是否采用MH算法，默认为True
'''
from __future__ import division
import matplotlib.pyplot as plt
import matplotlib as mpl
import numpy as np
from array import array
def mcmc(p,N=10000,Nlmax=10000,isMH=True):A = np.array([[1.0 / len(p) for x in range(len(p))] for y in range(len(p))], dtype=np.float64) X0 = np.random.randint(len(p))count = 0samplecount = 0L = array("d",[X0])l = array("d")while True:X = L[samplecount]cur = np.argmax(np.random.multinomial(1,A[X]))count += 1if isMH:a = (p[cur]*A[cur][X])/(p[X]*A[X][cur])alpha = min(a,1)else:alpha = p[cur]*A[cur][X]u = np.random.uniform(0,1) if u<alpha:samplecount += 1L.append(cur)if count>N:l.append(cur)if len(l)>=Nlmax:breakelse:continueLa = np.frombuffer(L)la = np.frombuffer(l)return La,la
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,l2p = np.array([0.9,0.05,0.05])
a = mcmc(p,Nlmax=10000000)[1]
l1 = ['state%d'% x for x in range(len(p))]
plt.pie(count(a,len(p))[0],labels=l1,labeldistance=0.3,autopct='%1.2f%%')
plt.title("sampling")
plt.show()

图1可见该样本收敛于0.7691,0.1154,0.1155，与我们的输入的差距相差甚远

2. 把输入更改成p=[0.8,0.1,0.1]，以同样的迭代次数，输出同样数量样本可得到
图2也是同样有很大差距，可以用方差来衡量差距

3. 把输入更改成p=[0.6,0.2,0.2]，以同样的迭代次数，输出同样数量样本可得到
图3同样还是有着很大的差距

4. 最后把输入更改成p=[1/3,1/3,1/3]，以同样的迭代次数，输出同样数量样本可得到
图4这次准确率却很高

由以上四组实验，猜测使构造的转移矩阵的每一行与输入的分布p(x)相等时，可以达到较高的精度，所以，如下更改代码：

#coding=utf-8
'''
Created on 2016年9月15日
基本MCMC算法以及M-H算法的python实现(离散数据的情况）
@author: whz
p:输入的概率分布，离散情况采用元素为概率值的数组表示
构造的转移矩阵每一行都与输入的p相同
N:认为迭代N次马尔可夫链收敛
Nlmax:马尔可夫链收敛后又取的服从p分布的样本数
isMH:是否采用MH算法，默认为True
'''
from __future__ import division
import matplotlib.pyplot as plt
import numpy as np
from array import array
def mcmc(p,N=10000,Nlmax=10000,isMH=True):A = np.array([p for y in range(len(p))], dtype=np.float64) X0 = np.random.randint(len(p))count = 0samplecount = 0L = array("d",[X0])l = array("d")while True:X = L[samplecount]cur = np.argmax(np.random.multinomial(1,A[X]))count += 1if isMH:a = (p[cur]*A[cur][X])/(p[X]*A[X][cur])alpha = min(a,1)else:alpha = p[cur]*A[cur][X]u = np.random.uniform(0,1) if u<alpha:samplecount += 1L.append(cur)if count>N:l.append(cur)if len(l)>=Nlmax:breakelse:continueLa = np.frombuffer(L)la = np.frombuffer(l)return La,la
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,l2if __name__ == '__main__':    p = np.array([0.9,0.05,0.05])a = mcmc(p)[1]l1 = ['state%d'% x for x in range(len(p))]plt.pie(count(a,len(p))[0],labels=l1,labeldistance=0.3,autopct='%1.2f%%')plt.title("sampling")plt.show()

得到如下：

与图1相比，很明显，图5准确率提高很多，依照同样程序，分别更改p，同样也可以得到对应于图2、3、4的准确率高的结果。

结论

要将状态转移矩阵的每一行都设计成与输入概率分布相一致，这样该算法才能得到更好的结果。
该算法该实现的适用范围：输入概率分布为离散的，有限的。
思考：虽然无限的离散型的应该也能按照上述的模式计算，但是对于无穷矩阵的计算，本小白还没研究，等着日后思考；而对于连续分布，则没有所谓的概率分布，而替代以概率密度，而且，由于在实现过程中使用了随机数生成，所以可以轻易通过变换随机数的方式生成的连续性分布的样本如果使用mcmc算法未免有点大炮轰蚊子，但若是考虑那些无法显示表出概率密度的连续型分布，我实在是还没想到好的解决办法，待后续解决。

MCMC算法—MH算法的Python实现（一）相关推荐

PRML第十一章读书笔记——Sampling Methods 拒绝采样/重要性采样/采样重要性重采样/数据增广IP算法、Metropolis算法/MH算法/吉布斯、切片采样、混合MC、估计配分函数
(终于把第十章读完了,这一章应该相对轻松.但这两天状态有待调整,所以没咋认真读) 目录 11.1 Basic Sampling Algorithms P526 标准概率分布 P528 拒绝采样 P53 ...
MCMC笔记Metropilis-Hastings算法（MH算法）
1 前言我们在MCMC笔记:齐次马尔可夫链_UQI-LIUWJ的博客-CSDN博客中介绍了平稳条件,当马尔可夫链达到平稳状态时(也就是各个状态之间的转移概率已经和时间无关了),那我们可以通过此时的 ...
机器学习经典算法详解及Python实现--元算法、AdaBoost
http://blog.csdn.net/suipingsp/article/details/41822313 第一节,元算法略述遇到罕见病例时,医院会组织专家团进行临床会诊共同分析病例以判定结果. ...
机器学习算法清单！附Python和R代码
来源:数据与算法之美本文约6000字,建议阅读8分钟. 通过本文为大家介绍了3种机器学习算法方式以及10种机器学习算法的清单,学起来吧~ 前言谷歌董事长施密特曾说过:虽然谷歌的无人驾驶汽车和机器人 ...
《大厂算法面试题目与答案汇总，剑指offer等常考算法题思路，python代码》V1.0版...
为了进入大厂,我想很多人都会去牛客.知乎.CSDN等平台去查看面经,了解各个大厂在问技术问题的时候都会问些什么样的问题. 在看了几十上百篇面经之后,我将算法工程师的各种类型最常问到的问题都整理了出来, ...
相似图片检测：感知哈希算法之dHash的Python实现
原文:https://blog.csdn.net/haluoluo211/article/details/52769325 相似图片检测:感知哈希算法之dHash的Python实现某些情况下,我们需 ...
python算法与程序设计基础第二版-算法与程序设计基础(Python版) - 吴萍
基本信息书名:21世纪高等学校计算机基础实用规划教材:算法与程序设计基础(Python版) 定价:39.00元作者:吴萍21世纪高出版社:清华大学出版社出版日期:2015_2_1 ISBN:9 ...
机器学习算法一览（附python和R代码）
机器学习算法一览(附python和R代码) 来源:数据观时间:2016-04-19 15:20:43 作者:大数据文摘 "谷歌的无人车和机器人得到了很多关注,但我们真正的未来却在于能 ...
2021-03-15 数据挖掘算法—K-Means算法 Python版本
数据挖掘算法-K-Means算法 Python版本简介又叫K-均值算法,是非监督学习中的聚类算法. 基本思想 k-means算法比较简单.在k-means算法中,用cluster来表示簇:容易证明 ...
一文读懂FM算法优势，并用python实现
介绍我仍然记得第一次遇到点击率预测问题时的情形,在那之前,我一直在学习数据科学,对自己取得的进展很满意,在机器学习黑客马拉松活动中也开始建立了自信,并决定好好迎接不同的挑战. 为了做得更好,我购买了 ...

MCMC算法—MH算法的Python实现（一）

针对离散数据概率分布的MCMC算法python实现

结论

MCMC算法—MH算法的Python实现（一）相关推荐

最新文章

热门文章