模式识别与机器学习(Python实现):如何用MCMC产生任意的概率分布随机数?Python简单实现

模式识别与机器学习(Python实现):如何用MCMC产生任意的概率分布随机数?Python简单实现


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

如何产生任意的概率分布?本文采用MCMC的方法进行演示。这个任务在进行算法仿真时具有很好的意义,首先我们需要设定一个所需的分布(需满足的要求是 ∫ − ∞ + ∞ p ( x ) d x   = 1 \int_{ - \infty }^{ + \infty } {p(x)dx{\text{ }} = 1} +p(x)dx =1),本文以指数分布为例,主要阐述实现的过程,如需更详细的解释,请参考相关的文献。
在这里插入图片描述
实现过程主要为6步。

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

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

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()

在这里插入图片描述

  • 6
    点赞
  • 7
    收藏
    觉得还不错? 一键收藏
  • 打赏
    打赏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

CV视界

如果感觉有用,可以打赏哦~

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值