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

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

---

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

1. 输入为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 = 0
    samplecount = 0
    L = array("d",[X0])
    l = array("d")
    while True:
        X = L[samplecount]