对于一个随机变量X，对其直接进行抽样比较困难时，我们可以选择一个比较容易产生随机数且逼近f很好的一个分布f(y)来对其进行抽样，下面以贝塔分布为例进行舍选抽样
ps:实际上,应该是先找了一个f(y),使得cf(y)的图像范围可以包括p(x)的函数图像，这样我们对p(x)进行随机抽样，位于p(x)内的x保留，位于p(x)外的x舍去，留下的x就是我们需要的p(x)的抽样
代码如下：
gx=function(x){
if(x<=0.8)
y=3*x
else if(x>0.8)
y=12*(1-x)
}#定义gx
be <- function(n){
n<-1000
j<-k<-0#给定j和k的初始值，以便对其进行计数
num<-0
y<-numeric(n)
while(k < n){
r2<-runif(1)#产生r1服从均匀分布
j<-j+1
r1<-runif(1)#从g中产生一个随机数
if(r1<=0.8){
x<-(0.8*r1)^(1/2)
}else{
x<-1-(0.2*(1-r1))^(1/2)
}#步骤2
if((12*x^2*(1-x))/gx(x)>r2){#舍选规则
k<-k + 1
y[k]<-x
}
num<-num+1
}
re<-list(y=y,num=num)
return(re)
}
hist(be(1200)$y, prob = TRUE)
xx = seq(0,1, len=1000)
yy = 12*xx^2*(1-xx)
lines(xx, yy, col = "red") #验证
re = be(1200)
还可以对其计算循环次数：
#计算实际循环次数
gx=function(x){
if(x<=0.8)
y=3*x
else if(x>0.8)
y=12*(1-x)
}#定义gx
n<-1000
j<-k<-0#给定j和k的初始值，以便对其进行计数
y<-numeric(n)
while(k < n){
r2<-runif(1)#产生r1服从均匀分布
j<-j+1
r1<-runif(1)#从g中产生一个随机数
if(r1<=0.8){
x<-(0.8*r1)^(1/2)
}else{
x<-1-(0.2*(1-r1))^(1/2)
}#步骤2
if((12*x^2*(1-x))/gx(x)>r2){#舍选规则
k<-k + 1
y[k]<-x
}
}
j
算法的舍选效率为1000/j。
参考文献：
(31条消息) R语言-舍选法产生正态分布随机数_Yuan..的博客-CSDN博客_舍选抽样法产生随机数
R语言随机数产生方法 - 知乎 (zhihu.com)
超星-直播课堂 (chaoxing.com)
