模式识别与机器学习（Python实现）：如何用MCMC产生任意的概率分布随机数？Python简单实现

模式识别与机器学习（Python实现）：如何用MCMC产生任意的概率分布随机数？Python简单实现

欢迎大家来到安静到无声的《模式识别与人工智能（程序与算法）》，如果对所写内容感兴趣请看模式识别与人工智能（程序与算法）系列讲解 - 总目录，同时这也可以作为大家学习的参考。欢迎订阅，优惠价只需9.9元，请多多支持！

如何产生任意的概率分布？本文采用MCMC的方法进行演示。这个任务在进行算法仿真时具有很好的意义，首先我们需要设定一个所需的分布（需满足的要求是${\int }_{-\infty }^{+\infty }p(x)dx=1$），本文以指数分布为例，主要阐述实现的过程，如需更详细的解释，请参考相关的文献。

实现过程主要为6步。

1. 首先，给定一个初始状态$x_0$，(随机分布产生)。
2. 对于$t=1,2,3,....$，我们循环以下过程进行采样，可以以10000个数据为例。
3. 以$t$为例，在$t$时刻进行采样时，采样以高斯分布为例$Q\left(x|x_t\right)=y$，其中$x_t$放在高斯分布的均值处。方差为1。
4. 我们采样得到了$y$，再从$[0,1]$的均匀分布中采样得到$u$，即$u\sim U(0,1)$
5. 计算$\alpha$的值，$\alpha =f\left(y\right)q\left(x_t|y\right)$，其中$f\left(y\right)$是上文已知的一个概率分布。
6. 判断$u<\alpha$，若是令$x_{t+1}=y$，否则令$x_{t+1}=x_t$。

上面的方法可能会陷入局部最小值，所以可以采用metropolis hastings的采样方法，区别是计算$\alpha$的方法不同，即为$$\alpha =\min \left\{\frac{f(y)q(x_t|y)}{f(x_t)q(y|x_t)},1\right\}$$
由此可得表示代码：

import numpy as np
np.random.seed(0)
x_t = np.random.rand(1)  
import math
def function(x):
    Lambda = 2
    e = 2.718281828459
    if x > 0:
        y = Lambda * e ** (-Lambda * x)
    else:
        y = 0
    return(y)
list_x_t = [float(x_t)]
for i in range(100000):
    y=np.random.normal(loc=x_0 , scale=1.0, size=None)  #loc:均值 scale:标准差
    u = np.random.rand(1)
    alpha = min(float(function(y)* np.random.normal(loc=y , scale=1.0)/function(x_t)* np.random.normal(loc=x_t , scale=1.0)),1.)
    if u < alpha:
        x_t = y
    else:
        x_t = x_t
    list_x_t.append(float(x_t))
    

import matplotlib.pyplot as plt
 
plt.hist(list_x_t,bins=1000)
plt.title("Data distribution histogram")
plt.show()