随机采样之接受拒绝采样

之前提到的逆变换采样(Inverse Transform Sampling)是一种生成随机样本的方法,它利用累积分布函数(CDF)的逆函数来生成具有特定分布的随机变量。以下是逆变换采样的缺点:

  1. 计算复杂性:对于某些分布,找到累积分布函数(CDF)的逆函数可能是困难的,甚至是不可能的。
  2. 效率问题:对于具有重尾分布的随机变量,逆变换采样可能非常低效,因为CDF的逆可能需要大量的计算。
  3. 数值稳定性:在数值计算中,由于浮点数的精度限制,逆变换采样可能会引入误差,尤其是在CDF的值接近1时。

一、接受拒绝采样

接受-拒绝采样(Accept-Reject Sampling)方法是一种更为通用的采样方法,它可以用来生成具有任意分布的随机样本。这种方法不要求我们知道CDF的逆,而是利用一个简单的概率分布(称为提议分布)来生成样本,然后以一定的概率接受或拒绝这些样本。

接收-拒绝采样的基本步骤:

  1. 选择提议分布 g ( x ) g(x) g(x):选择一个容易从中抽样的分布 g ( x ) g(x) g(x),并且确保对于所有的 x x x,有 f ( x ) ≤ M ⋅ g ( x ) f(x) \leq M \cdot g(x) f(x)Mg(x),其中 f ( x ) f(x) f(x)是目标分布, M M M是一个正常数。

  2. 抽样:从提议分布 g ( x ) g(x) g(x)中抽取样本 x x x和从均匀分布 U ( 0 , 1 ) U(0, 1) U(0,1)中抽取样本 u u u

  3. 接受-拒绝条件:如果 u ≤ f ( x ) M ⋅ g ( x ) u \leq \frac{f(x)}{M \cdot g(x)} uMg(x)f(x),则接受 x x x作为目标分布 f ( x ) f(x) f(x)的一个样本;否则拒绝 x x x

接受拒绝采样可以使用下图进行表示(图片来源:【数之道】马尔可夫链蒙特卡洛方法是什么?十五分钟理解这个数据科学难点)。
在这里插入图片描述

二、接受拒绝采样证明

要证明接收-拒绝采样确实产生服从目标分布 f ( x ) f(x) f(x)的样本,我们需要证明对于所有的 x x x,有:
P ( X = x ) = f ( x ) (1) P(X=x) = f(x)\tag1 P(X=x)=f(x)(1)

其中 P ( X = x ) P(X=x) P(X=x)是样本 x x x被接受的概率。

证明:
  1. 接受概率:样本 x x x被接受的概率是 f ( x ) M ⋅ g ( x ) \frac{f(x)}{M \cdot g(x)} Mg(x)f(x),因为 u u u是从 U ( 0 , 1 ) U(0, 1) U(0,1)中抽取的。

  2. 联合概率:样本 x x x从提议分布 g ( x ) g(x) g(x)中抽取的概率是 g ( x ) g(x) g(x),并且 u u u [ 0 , f ( x ) M ⋅ g ( x ) ) [0, \frac{f(x)}{M \cdot g(x)}) [0,Mg(x)f(x))区间的概率是 f ( x ) M ⋅ g ( x ) \frac{f(x)}{M \cdot g(x)} Mg(x)f(x)。因此,联合概率是:

    P ( X = x , U ≤ f ( x ) M ⋅ g ( x ) ) = g ( x ) ⋅ f ( x ) M ⋅ g ( x ) = f ( x ) M (2) P(X=x, U \leq \frac{f(x)}{M \cdot g(x)}) = g(x) \cdot \frac{f(x)}{M \cdot g(x)} = \frac{f(x)}{M}\tag2 P(X=x,UMg(x)f(x))=g(x)Mg(x)f(x)=Mf(x)(2)

  3. 边缘概率:现在我们需要计算 X X X的边缘概率 P ( X = x ) P(X=x) P(X=x),即样本 x x x被接受的总概率。由于 u u u是均匀分布的,我们可以将联合概率在 u u u的所有可能值上积分:

    P ( X = x ) = ∫ 0 1 P ( X = x , U = u ) d u = ∫ 0 1 f ( x ) M d u = f ( x ) M ⋅ ∫ 0 1 d u = f ( x ) M (3) P(X=x) = \int_0^1 P(X=x, U=u) \, du = \int_0^1 \frac{f(x)}{M} \, du = \frac{f(x)}{M} \cdot \int_0^1 du = \frac{f(x)}{M}\tag3 P(X=x)=01P(X=x,U=u)du=01Mf(x)du=Mf(x)01du=Mf(x)(3)

  4. 归一化常数:由于 M M M是使得 f ( x ) ≤ M ⋅ g ( x ) f(x) \leq M \cdot g(x) f(x)Mg(x)对所有 x x x成立的最小常数,我们可以将上式中的 M M M移到 f ( x ) f(x) f(x)的定义中,从而得到:

    P ( X = x ) = f ( x ) (4) P(X=x) = f(x)\tag4 P(X=x)=f(x)(4)
    这就证明了接收-拒绝采样确实产生了服从目标分布 f ( x ) f(x) f(x)的样本。

三、接受拒绝采样模拟

借用作者anshuai_aw1的例子,设我们需要采样的pdf为:
f ( x ) = 0.3 exp ⁡ ( − ( x − 0.3 ) 2 ) + 0.7 exp ⁡ ( − ( x − 2 ) 2 / 0.3 ) (5) f(x)=0.3 \exp \left(-(x-0.3)^{2}\right)+0.7 \exp \left(-(x-2)^{2} / 0.3\right)\tag5 f(x)=0.3exp((x0.3)2)+0.7exp((x2)2/0.3)(5)
其归一化常数为 Z = 1.2113 Z = 1.2113 Z=1.2113, 参考分布为 g ( x ) = N ( μ = 1.4 , σ 2 = ( 1. 2 2 ) ) g(x) =N(\mu=1.4,\sigma^2=(1.2^2)) g(x)=N(μ=1.4,σ2=(1.22)), M = 2.5 M=2.5 M=2.5, 以确保 M ⋅ g ( x ) ≥ f ( x ) M \cdot g(x) \geq f(x) Mg(x)f(x)。采样的代码如下:

import numpy as np
import matplotlib.pyplot as pltdef f(x):return (0.3*np.exp(-(x-0.3)**2) + 0.7* np.exp(-(x-2.)**2/0.3))/1.2113
x = np.arange(-4.,6.,0.01)
plt.plot(x,f(x),color = "red")size = int(1e+07)
mu = 1.4
sigma = 1.2
M = 2.5x = np.random.normal(loc = mu,scale = sigma, size = size)
g_x = 1/(np.sqrt(2*np.pi)*sigma)*np.exp(-0.5*(x-mu)**2/sigma**2)
u = np.random.uniform(low = 0, high = M*g_x, size = size)  #在[0,M*g_x]中均匀采样
fx =  0.3*np.exp(-(x-0.3)**2) + 0.7* np.exp(-(x-2.)**2/0.3)
sample = x[u <= fx] # u < fx(x)
plt.hist(sample,bins=150, density=True, edgecolor='black')
plt.show()

结果如下,其中红色曲线的是公式(5)所示pdf的图像,蓝色区域是采样结果,可见采样结果跟真实分布几乎一致。
在这里插入图片描述

参考资料:

[1]【数之道】马尔可夫链蒙特卡洛方法是什么?十五分钟理解这个数据科学难点
[2] 逆采样(Inverse Sampling)和拒绝采样(Reject Sampling)原理详解

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.mzph.cn/bicheng/59812.shtml

如若内容造成侵权/违法违规/事实不符,请联系多彩编程网进行投诉反馈email:809451989@qq.com,一经查实,立即删除!

相关文章

用 Python 爬取淘宝商品价格信息时需要注意什么?

用 Python 爬取淘宝商品价格信息时&#xff0c;需要注意以下方面&#xff1a; 一、法律和道德规范&#xff1a; 遵守法律法规&#xff1a;网络爬虫的行为应在法律允许的范围内进行。未经淘宝平台授权&#xff0c;大规模地爬取其商品价格信息并用于商业盈利等不当用途是违法的…

免费数据集网站

1、DataSearch https://datasetsearch.research.google.comhttp://DataSearch 2、FindData findata-科学数据搜索引擎https://www.findata.cn/ 3、Kaggle Kaggle: Your Machine Learning and Data Science CommunityKaggle is the world’s largest data science community …

在 FPGA 中实现 `tanh` 和 Softplus 函数

✅作者简介&#xff1a;2022年博客新星 第八。热爱国学的Java后端开发者&#xff0c;修心和技术同步精进。 &#x1f34e;个人主页&#xff1a;Java Fans的博客 &#x1f34a;个人信条&#xff1a;不迁怒&#xff0c;不贰过。小知识&#xff0c;大智慧。 &#x1f49e;当前专栏…

基于java+SpringBoot+Vue的旅游管理系统设计与实现

项目运行 环境配置&#xff1a; Jdk1.8 Tomcat7.0 Mysql HBuilderX&#xff08;Webstorm也行&#xff09; Eclispe&#xff08;IntelliJ IDEA,Eclispe,MyEclispe,Sts都支持&#xff09;。 项目技术&#xff1a; Springboot mybatis Maven mysql5.7或8.0等等组成&#x…

【Python TensorFlow】入门到精通

TensorFlow 是一个开源的机器学习框架&#xff0c;由 Google 开发&#xff0c;广泛应用于机器学习和深度学习领域。本篇将详细介绍 TensorFlow 的基础知识&#xff0c;并通过一系列示例来帮助读者从入门到精通 TensorFlow 的使用。 1. TensorFlow 简介 1.1 什么是 TensorFlow…

数据库管理-第258期 23ai:Oracle Data Redaction(20241104)

数据库管理258期 2024-11-04 数据库管理-第258期 23ai&#xff1a;Oracle Data Redaction&#xff08;20241104&#xff09;1 简介2 应用场景与有点3 多租户环境4 特性与能力4.1 全数据编校4.2 部分编校4.3 正则表达式编校4.4 随机编校4.5 空值编校4.6 无编校4.7 不同数据类型上…

基于SpringBoot的医药管理系统+LW示例参考

1.项目介绍 系统角色&#xff1a;管理员、收银员功能模块&#xff1a;管理员&#xff08;收银员信息管理、药品管理、药品类别、出库信息管理、入口信息。药品库存图表&#xff09;、收银员&#xff08;药品库存图表、会员积分信息等&#xff09;技术选型&#xff1a;SpringBo…

PH热榜 | 2024-11-07

DevNow 是一个精简的开源技术博客项目模版&#xff0c;支持 Vercel 一键部署&#xff0c;支持评论、搜索等功能&#xff0c;欢迎大家体验。 在线预览 1. SWE-Kit 标语&#xff1a;打造你自己的“德文”——一个像软件工程师一样的智能助手&#xff01; 介绍&#xff1a;SWE-K…

(蓝桥杯C/C++)——基础算法(下)

目录 一、时空复杂度 1.时间复杂度 2.空间复杂度 3.分析技巧 4.代码示例 二、递归 1.递归的介绍 2.递归如何实现 3.递归和循环的比较 4.代码示例 三、差分 1.差分的原理和特点 2.差分的实现 3.例题讲解 四、枚举 1.枚举算法介绍 2.解空间的类型 3. 循环枚举解…

echarts功能五 --geo地理组件、VisualMap图例组件

利用geoJson文件生成geo地理组件&#xff0c;如下图所示&#xff1a; 三个颜色区域分别代表了3个区域图层&#xff1b;淡蓝色代表了线条&#xff1b;正中心是geo地理组件&#xff1b;右下角展示图例&#xff0c;是VisualMap视觉映射组件。 共包含以下功能&#xff1a; &#…

WordCloudStudio:AI生成模版为您的文字云创意赋能 !

在信息泛滥的时代&#xff0c;如何有效地将文字内容变成生动的视觉元素&#xff1f;WordCloudStudio为您提供了答案。无论您是市场营销专家、教育工作者、数据分析师&#xff0c;还是创意设计师&#xff0c;WordCloudStudio都能帮助您轻松创建引人注目的文字云。更重要的是&…

25-RVIZ CARLA插件

RVIZ插件(RVIZ plugin)提供了一个基于RVIZ(RVIZ) ROS包的可视化工具。 用RVIZ运行ROS桥接 RVIZ插件需要一个名为ego_vehicle的自车。要查看ROS-bridge使用RVIZ的示例&#xff0c;请在运行CARLA服务器的情况下执行以下命令&#xff1a; 1. 启用RVIZ启动ROS桥接&#xff1a; # …

FP7209单节锂电升压恒流80V,PWM控制调光调色应急电源驱动方案,支持LED开路保护、LED短路保护、开关NMOS过电流保护、过温保护、过热保护

FP7209是针对LED驱动器的升压拓扑开关调节器。它提供了内置的门驱动销&#xff0c;用于驱动外部N-MOSFET。误差放大器的非反相输入端连接到一个0.25V的参考电压。如UVP、OVP、OCP等&#xff0c;保护系统电路有三个功能。LED电流可以通过一个连接到DIM针脚的外部信号来调整。DIM…

JS常用数组方法 reduce filter find forEach

文章目录 reduce应用&#xff1a;数据扁平化 filterfind从数组 [1,2,3,4,5,6] 中找出值为 2 的元素 forEach用于遍历&#xff0c;forEach 方法没有返回值&#xff0c;它总是返回 undefined。 reduce 数组变量名.reduce((sum,value) > { // 向sum变量上累加值 // 一定要retur…

精选报告| 2024年,5份必读的“虚仿教育”行业报告合集

以3D/XR应用为主的虚拟仿真实验教学课程&#xff0c;在教育信息化建设过程中已成为必选的技术方案。通过构建虚拟教育环境&#xff0c;允许学习者在数字空间中进行互动学习&#xff0c;这种方法在基础教育、职业培训、远程教育等关键教育领域已经展现出前所未有的变革潜力&…

Ethernet 系列(8)-- 基础学习::ARP

目录 1. ARP的目的&#xff1a; 1.1 什么是ARP 1.2 什么时候用ARP 2. ARP如何工作&#xff1a; 2.1 主机-主机的直接通信 2.2 主机-路由-主机的间接通信 3. ARP header&#xff1a; 1. ARP的目的&#xff1a; 1.1 什么是ARP: ARP-地址解析协议&#xff0c;是第3层地址&#xff…

uniapp组件实现省市区三级联动选择

1.导入插件 先将uni-data-picker组件导入我们的HBuilder项目中&#xff0c;在DCloud插件市场搜索uni-data-picker 点击下载插件并导入到我们的项目中 2.组件调用 curLocation &#xff1a;获取到的当前位置&#xff08;省市区&#xff09; <uni-data-picker v-slot:defa…

软件分享丨火绒应用商店

【资源分享】 资源名&#xff1a;火绒应用商店 官方网址&#xff1a;点击跳转 火绒应用商店是由火绒安全推出的一款独立软件。它提供了海量的应用程序&#xff0c;涵盖办公、社交、游戏、视频、工具等多种领域和类别&#xff0c;方便用户轻松找到所需的应用并进行一键下载安装…

信息化运维方案,实施方案,开发方案,信息中心安全运维资料(软件资料word)

1 编制目的 2 系统运行维护 2.1 系统运维内容 2.2 日常运行维护方案 2.2.1 日常巡检 2.2.2 状态监控 2.2.3 系统优化 2.2.4 软件系统问题处理及升级 2.2.5 系统数据库管理维护 2.2.6 灾难恢复 2.3 应急运行维护方案 2.3.1 启动应急流程 2.3.2 成立应急小组 2.3.3 应急处理过程 …

鸿蒙ArkTS中的布局容器组件(Column、Row、Flex、 Stack、Grid)

在鸿蒙ArkTS中&#xff0c;布局容器组件有很多&#xff0c;常见的有&#xff1a;   ⑴ Column&#xff1a;&#xff08;垂直布局容器&#xff09;&#xff1a;用于将子组件垂直排列。   ⑵ Row&#xff1a;&#xff08;水平布局容器&#xff09;&#xff1a;用于将子组件水…