## MCMC_MH 对连续和离散函数应用实例

MCMC_MH 对连续和离散函数应用实例

MCMC_MH 对离散函数采样

离散函数采样来自作者kissmeanus
https://blog.csdn.net/kissmeanus/article/details/89073551

'''
Pi:输入概率分布，离散情况采用元素为概率值的数组表示
N:认为迭代N次马尔可夫链收敛
Nlmax:马尔可夫链收敛后又取的服从p分布的样本数
isMH:是否采用MH算法，默认为True
'''
import matplotlib.pyplot as plt
import numpy as np
from array import array


def mcmc(P,Q,N = 1000, Nlmax = 100000, isMH = False):
    T = N + Nlmax - 1
    result = [0 for i in range(T)]
    X0 = np.random.randint(len(P))
    result[0] = X0
    t = 0
    while t < T - 1:#还需要采样T-1次
        t = t+1
        x_cur = np.argmax(np.random.multinomial(1,Q[result[t-1]]))
        if isMH:
            a = (P[x_cur] * Q[x_cur][result[t-1]]) / (P[result[t-1]] * Q[result[t-1]][x_cur])
            acc = min(1,a)
        else:
            acc = min(1,P[x_cur] * Q[x_cur][result[t-1]])
        u = np.random.uniform(0,1)
        if u < acc:
            result[t] = x_cur
        else:
            result[t] = result[t-1]
    result_chuli = result[N:]
    return result_chuli

def count(a, n):
    l1 = array('d')
    l2 = array('d')
    l3 = array('d')
    for e in a:
        l1.append(e)
    for e in range(n):
        l2.append(l1.count(e))
    for e in l2:
        l3.append(e/sum(l2))
    return l2,l3
    
if __name__ == '__main__':
    Pi = [0.2,0.3,0.5]
    Q = np.array([[0.2,0.3,0.5],
                  [0.2,0.3,0.5],
                  [0.2,0.3,0.5]])
    a = mcmc(Pi,Q,isMH = True)
    
    #最后用饼图来表述所采样本与目标函数Pi之间的关系
    l1 = count(a, len(Pi))[0]#l1获得各个状态的次数
    l2 = ['state%d' % x for x in range(len(Pi))]
    explode_list = [0,0.05,0]
    plt.pie(l1, labels= l2, labeldistance= 0.1, explode=explode_list, autopct='%0.2f%%')
    plt.legend()
    plt.show()

MCMC_MH 对连续函数进行采样

'''
MCMC_MH 对连续函数进行采样
Pi:输入概率分布，离散情况采用元素为概率值的数组表示
N:认为迭代N次马尔可夫链收敛
Nlmax:马尔可夫链收敛后又取的服从p分布的样本数
isMH:是否采用MH算法，默认为True
连续函数采样中的状态转移矩阵Q
'''
import matplotlib.pyplot as plt
import numpy as np
import scipy.stats

def Pi_hanshu(x, loc= 0, scale= 1):
    '''
    目标函数
    '''
    Pi = scipy.stats.norm.pdf(x, loc=loc,scale=scale)
    return Pi
def mcmc(N=1000, Nlmax=100000, isMH=False):
    T = N + Nlmax - 1
    result = [0 for i in range(T)]
    X0 = scipy.stats.norm.rvs(loc=0, scale=1, size= 1)
    result[0] = X0[0]
    t = 0
    while t < T - 1:
        t += 1
        x_cur = (scipy.stats.norm.rvs(loc=result[t-1], scale= 1, size= 1))[0]
        if isMH:
            a = (Pi_hanshu(x_cur) * (scipy.stats.norm.pdf(x_cur,loc=result[t-1], scale= 1)))\
                / (Pi_hanshu(result[t-1]) * (scipy.stats.norm.pdf(result[t-1], loc=x_cur, scale= 1)))
            acc = min(1,a)
        else:
            acc = min(1,Pi_hanshu(x_cur) * scipy.stats.norm.pdf(x_cur, loc=result[t-1], scale= 1))
        u = np.random.uniform(0,1)
        if u < acc:
            result[t] = x_cur
        else:
            result[t] = result[t-1]
    result_chuli = result[N:]
    return result_chuli

def dia_scatter(a):
    x = [i for i in range(len(a))]
    plt.scatter(x, a, marker='x', c='r')
    plt.show()

if __name__ == '__main__':
    a = mcmc()
    
	最后用散点图来表示所采样本与目标函数N(0,1)之间的关系
    dia_scatter(a)

从图上可以看出，所采样本围绕x轴服从N（0,1）分布