随机采样之接受拒绝采样

之前提到的逆变换采样(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…

C++ 并发专题 - 条件变量的使用

一&#xff1a;概述&#xff1a; 在 C 中&#xff0c;条件变量&#xff08;std::condition_variable&#xff09;是一种用于线程间同步的机制&#xff0c;主要用于在多线程环境中让一个线程等待某个条件满足后再继续执行。条件变量通常配合互斥锁&#xff08;std::mutex&#…

【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…

获取服务器相关信息

获取主机厂商型号信息&#xff0c;如有需要&#xff0c;可以根据获取到的厂商型号&#xff0c;去官网获取产品规格书 sudo dmidecode -t system# dmidecode 3.2 Getting SMBIOS data from sysfs. SMBIOS 3.2.0 present.Handle 0x0001, DMI type 1, 27 bytes System Informatio…

经典的ORACLE 11/12/19闪回操作

1、闪回表 SQL> show parameter recycle NAME TYPE VALUE ------------------------------------ ----------- ------------------------------ buffer_pool_recycle string db_recycle_cache_size …

PH热榜 | 2024-11-07

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

注册了个域名Wordpress.cn.com

wordpress.cn.com 这个域名一看是很山寨的&#xff01; 我玩wordpress比较多&#xff0c;所以就想有一个wordpress之类的域名拿着耍耍&#xff0c;就像我的typecho.pro之类的那样&#xff0c;所以就想着 wp/wordpress之类的类型选一个&#xff0c;看了不少国别和新后缀&…

浏览器内置对象XMLHttpRequest

XMLHttpRequest 是浏览器提供的一个内置对象&#xff0c;用于在客户端和服务器之间进行异步通信。它是一种由浏览器提供的 JavaScript API&#xff08;应用程序编程接口&#xff09;&#xff0c;开发者可以通过 JavaScript 代码来使用它。通过 XMLHttpRequest&#xff0c;开发者…

(蓝桥杯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…

Spring JDBC模板

Spring JDBC模板&#xff08;JdbcTemplate&#xff09;是Spring框架提供的一个简化JDBC操作的工具类。它封装了JDBC的常见操作&#xff0c;如查询、更新、插入和删除等&#xff0c;简化了数据库访问代码&#xff0c;减少了样板代码。下面是一个详细的示例&#xff0c;展示如何使…

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…