接受拒绝法Acceptance-Rejection Method

原理解释
具体步骤
连续情况举例
离散情况举例

计算机能够通过调用自带函数，实现对均匀分布等常见分布的采样。但在实际情况中分布通常由其概率密度函数（PDF）表示，所以在已知PDF的情况下，如何让计算机生成我们想要的观测值？
接受拒绝方法（Acceptance-Rejection Method）也称拒绝采样，该方法用于生成服从某个概率密度函数的随机数，是一种蒙特卡罗方法（MC/R语言：蒙特卡洛方法求积分)

原理解释

我们想要获取下图中红色密度曲线f(x)f(x)f(x) 的采样数据，可以先用一个矩形 g(x)g(x)g(x) 将密度曲线框在里面，并在矩形中随即投点，密度曲线下侧的点用绿色表示，上侧的点用蓝色表示：

不难看出，随机投点的这个动作生成的点，服从均匀分布，只需将密度曲线下侧的绿色点提取出来，即获得了服从红色密度曲线的采样数据。
采用均匀分布作为建议分布在某些情况下效率较低，从上图就可以看出，均匀分布中很多蓝色的点都被剔除了，也就是生成这些点是无效的，所以可以选择一些其他曲线来把密度曲线框起来，效率会提高一点，如下图中黄色密度曲线：

由此可知，建议分布g(x)g(x)g(x)只是获得目标密度函数曲线下均匀分布样本的一个辅助工具。目标分布f(x)与建议分布g(x)g(x)g(x)、常数c满足：
① 对建议分布g(x)g(x)g(x)的取样方便
② g(x)g(x)g(x)与f(x)f(x)f(x)形状类似，即：f(x)≤cg(x)f(x) \leq cg(x)f(x)≤cg(x)

具体步骤

已知密度函数 f(x)f(x)f(x)，建议函数 g(x)g(x)g(x)
step1：生成随机样本x服从g(x)g(x)g(x)，找c
step2: 生成服从均匀分布的随机数u
step3:若u≤f(x)cg(x)u \leq \frac{f(x)}{cg(x)}u≤cg(x)f(x) 则认为x可接受，否则返回step1，重新开始

连续情况举例

例一:
生成1000个服从f(x)=6x(1−x)，0<x<1f(x) = 6x(1-x)，0<x<1f(x)=6x(1−x)，0<x<1的随机数。

解析：
令g(x)g(x)g(x)为0到1之间的均匀分布，当0<x<10<x<10<x<1时，f(x)g(x)=6x(1−x)1≤6\frac{f(x)}{g(x)}=\frac{6x(1-x)}{1} \leq 6g(x)f(x)=16x(1−x)≤6，所以 c=6c=6c=6，
生成服从均匀分布的随机数u，当满足：f(x)cg(x)=6x(1−x)6×1>u\frac{f(x)}{cg(x)}=\frac{6x(1-x)}{6\times1}>ucg(x)f(x)=6×16x(1−x)>u
则接受，反之则拒绝.
代码实现:

n = 1000  # 题目要求的样本数目
k = 0  # 初始化满足要求的样本数目
j = 0  # 初始化循环次数
x = numeric(n) # 初始化样本列表
while(k < n){   # while次数不固定 for次数固定u = runif(1)  # 生成uj = j + 1  # 每循环一次+1y = runif(1)  # 生成服从g(x)的样本数据if(y * (1 - y) > u){  # 满足则接受，否则跳过，开始下一次循环k = k + 1  # 满足的样本数据个数+1x[k] = y  # 记录满足的样本数据}
}
x
hist(x,prob=T, ylim = c(0,1.6))  # 绘制直方图
xx = seq(0,1,0.01)
yy = 6 * xx *(1 - xx)  # 绘制真实的密度曲线
lines(xx,yy)

结果展示：

例二:
生成5000个服从f(x)=2x，0<x<1f(x) = 2x，0<x<1f(x)=2x，0<x<1的随机数。

g(x)g(x)g(x)为0到1之间的均匀分布，c=2c=2c=2
代码展示：

n = 5000
k = 0
j = 0
x = numeric(n)
while(k < n){   u = runif(1)  j = j + 1  y = runif(1)  if( y > u){  k = k + 1 x[k] = y  }
}
x
hist(x,prob=T, ylim = c(0,1.6))  # 绘制直方图
xx = seq(0,1,0.01)
yy = 6 * xx *(1 - xx)  # 绘制真实的密度曲线
lines(xx,yy)

结果展示：

离散情况举例

例一：
离散型随机变量X服从下列分布列，生成5000个服从该分布的随机数

xxx	0	1	2
p(x)p(x)p(x)	0.3	0.4	0.3

解析：
设建议函数g(x)g(x)g(x)服从下列分布列：

xxx	0	1	2
p(x)p(x)p(x)	0.25	0.25	0.25

因为c=max(f(x))g(x),u<f(x)g(x)×cc=\frac{max(f(x))}{g(x)},u<\frac{f(x)}{g(x) \times c}c=g(x)max(f(x)),u<g(x)×cf(x)，所以u<f(x)max(f(x))u< \frac{f(x)}{max(f(x))}u<max(f(x))f(x)

代码展示:

n = 5000
k = 0  # 初始化
j = 0  # 循环次数计数
x = numeric(n) # n个0
p = c(0.3, 0.4, 0.3)  # (x)
while(k < n){   # while次数不固定 for次数固定u = runif(1)j = j + 1y = sample(0:2, 1, replace = T)if(u < p[y+1]/(0.4)){k = k + 1x[k] = y}
}
x
random_p = table(x)/n  # 计算每个数字出现的频率
rbind(random_p,p)  # 将其与实际分布列对比

结果展示：

R语言：接受拒绝法（Acceptance-Rejection Method）生成随机数相关推荐

R语言：逆变换法生成随机数
逆变换法生成随机数: 一.概念解释 1.PDF 2.PMF 3.CDF 二.连续型情况举例三.离散型情况举例一.概念解释 1.PDF probability density function 概率 ...
R语言贝叶斯方法在生态环境领域中的高阶技术
贝叶斯统计学即贝叶斯学派是一门基本思想与传统基于频率思想的统计学即频率学派完全不同的统计学方法,它在统计建模中具有灵活性和先进性特点,使其可以轻松应对复杂数据和模型结构. 然而,很多初学者在面对思想. ...
R语言使用edit函数在Rsudio中生成数据编辑器(在windows中生成编辑器)、在编辑器中输出需要的数据生成最终的dataframe
R语言使用edit函数在Rsudio中生成数据编辑器(在windows中生成编辑器).在编辑器中输出需要的数据生成最终的dataframe 目录
R语言使用rgamma函数生成符合Gamma分布的随机数、使用plot函数可视化符合Gamma分布的随机数(Gamma Distribution)
R语言使用rgamma函数生成符合Gamma分布的随机数.使用plot函数可视化符合Gamma分布的随机数(Gamma Distribution) 目录
R语言科学计数法详解：digits和scipen设置
控制R语言科学计算法显示有两个option: digitis和scipen.介绍的资料很少,而且有些是错误的.经过翻看R语言的帮助和做例子仔细琢磨,总结如下: 默认的设置是: getOption(&q ...
r语言三倍标准差法去除异常值，再计算平均值标准差
博主自己没能找到好的函数去除异常值,于是自己写好了一个简单实用的包.可以通过三倍标准差法删去每一行的异常值,然后计算出平均值标准差. 函数总共四个参数: file= 要计算的文件路径,在工作目录可以 ...
R语言函数：length计算长度、seq生成数据序列、rep将数据对象重复N遍复制、cut将连续变量分割为多水平的因子变量、pretty将连续变量x分成n个区间创建合适的断点、cat数据对象拼接
R语言函数:length函数计算数据对象的长度.seq函数生成数据序列(sequenceÿ
r语言插补法_R语言用多重插补法估算相对风险
在这里,我将用R中的一个小模拟示例进行说明.首先,我们使用X1和X2双变量法线和Y模拟大型数据集,其中Y遵循给定X1和X2的逻辑模型. 首先,我们模拟一个非常大的完整数据集: #simulate完整数 ...
R语言-随机前沿分析法--SFA
3.1介绍生产函数模型: lnqi=x'i*b+vi-ui (随机生产前沿函数) qi:产出变量向量 x'i:投入变量向量 b:变量参数估计 vi:统计噪声的对 ...

R语言：接受拒绝法（Acceptance-Rejection Method）生成随机数

接受拒绝法Acceptance-Rejection Method

原理解释

具体步骤

连续情况举例

离散情况举例

R语言：接受拒绝法（Acceptance-Rejection Method）生成随机数相关推荐

最新文章

热门文章