大作业四：三种抽样算法及其效率比较

一、题目要求

设随机变量\(X\)的分布密度为\(p\left( x \right)=\frac{1}{2}+x,0\le x\le 1\)，分别给出用逆变换法，舍选法和复合法生成上述密度函数随机数的算法，并比较这三种算法的效率。

二、算法分析

2.1 逆变换法

2.1.1 算法要求

由随机变量\(X\)的密度函数\(p\left( x \right)=\frac{1}{2}+x\)，可以得到其分布函数： \[ F\left( x \right)=\int_{0}^{x}{p(x)dx}=\frac{{{x}^{2}}+x}{2}\tag{1} \] 该分布函数在定义域区间\(0\le x\le 1\)内连续且单调递增，必然存在反函数：

\[ {{F}^{-1}}(\mu )=\sqrt{2\mu +\frac{1}{4}}-\frac{1}{2}\tag{2} \] 满足逆变换法要求，于是可用逆变换法生成服从密度函数\(p\left( x \right)\)的随机变量\(X\)。

2.1.2 算法步骤

Step 1: 输入所要生成的随机数个数\(n\)；

Step 2: 生成\(n\)个服从\(U\left( 0,1 \right)\)的随机数，记做\(miu\)；

Step 3: 对\(miu\)做变换：\(\sqrt{2*miu+\frac{1}{4}}-\frac{1}{2}\)，记做\(y\)；

Step 4: 输出：数据列举、算法效率、ks检验；

Step 5: 绘制抽样密度曲线与理论密度曲线的对比图。

2.1.3 算法效率

定义抽样效率为：最终被接受数的占所有抽取随机数的比例。

由于按逆变换法抽取的每一个数都可以被选作所需要的随机数，因此逆变换法理论抽样效率是\(1\)。

2.2 舍选法

2.2.1 算法要求

\(p\left( x \right)\)是定义于有界闭区间\(\left[ 0,1 \right]\)的；有常数上界\(p\left( 1 \right)=\frac{3}{2}\)的概率密度函数。取\(a=0,b=1,c=\frac{3}{2}\)，因此可用简单舍选抽样法抽取服从密度函数\(p\left( x \right)\)的随机变量\(X\)。

2.2.2 算法步骤

Step 1: 输入所要生成的随机数个数\(n\)；

Step 2: 生成长度为\(n\)的零向量\(y\)；

Step 3: 计数变量\(i\leftarrow 0\)；

　　　循环次数\(j\leftarrow 0\)；

Step 4: 当\(i\le n\)时，生成服从\(U\left( 0,1 \right)\)的随机数\(x,r\)，循环次数\(j\leftarrow j+1\)；

　　　如果\(\frac{3}{2}r\le p(x)\)，接受\(x\)，并把它放置在\(y\)中，计数变量\(i\leftarrow i+1\)；

Step 5: 输出：数据列举、算法效率、ks检验；

Step 6: 绘制抽样密度曲线与理论密度曲线的对比图。

2.2.3 算法效率

定义抽样效率为：最终被接受数的占所有抽取随机数的比例。

由于简单舍选抽样法预先抽取的随机数可以视为在矩形闭区域内均匀分布。且x方向与y方向相互独立。因此容易得到抽样效率可以看作概率密度曲线与x轴所围面积与矩形闭区域面积的比值。由理论概率密度图像易得，简单舍选抽样法理论抽样效率\(\frac{2}{3}\)。

图一：舍选法采样图示

2.3 复合法

2.3.1 算法要求

对\(p(x)\)做变换：

\[ p\left( x \right)=\frac{1}{2}+x=\frac{1}{2}\times 1+\frac{1}{2}\times 2x\triangleq \frac{1}{2}{{p}_{1}}\left( x \right)+\frac{1}{2}{{p}_{2}}\left( x \right)\tag{3} \] 其中：

离散型随机变量\(I\)的取值： \[ p\left( I=i \right)=\frac{1}{2},i=1,2\tag{4} \] 并且： \[ \begin{equation} \int_{0}^{1}{{{p}_{1}}\left( x \right)dx}=\int_{0}^{1}{1dx}=1\\ \int_{0}^{1}{{{p}_{2}}\left( x \right)dx}=\int_{0}^{1}{2xdx}=1\tag{5} \end{equation} \]

可见\({{p}_{1}}\left( x \right),{{p}_{2}}\left( x \right)\)均为概率密度函数，不妨设\({{Z}_{1}},{{Z}_{2}}\)分别是服从概率密度为\({{p}_{1}}\left( x \right),{{p}_{2}}\left( x \right)\)的随机变量。显然，有

\[ \begin{align} & {{Z}_{1}}\sim U\left( 0,1 \right) \\ & {{Z}_{2}}\sim \max \left\{ {{\xi }_{1}},{{\xi }_{2}} \right\}\left( {{\xi }_{1}},{{\xi }_{2}}\sim U\left( 0,1 \right) \right) \\ \end{align}\tag{6} \]

那么：

\[ X=\left\{ \begin{align} & {{Z}_{1}}{\quad}ifI=1\\ & {{Z}_{2}}{\quad}ifI=2 \\ \end{align} \right.\tag{7} \]

所以，概率密度函数\(p\left( x \right)\)可以写做概率密度函数\({{p}_{1}}\left( x \right)与{{p}_{2}}\left( x \right)\)的线性叠加。因此可用复合法抽取服从密度函数\(p\left( x \right)\)的随机变量\(X\)。

2.3.2 算法步骤

Step 1: 输入所要生成的随机数个数\(n\)；

Step 2: 生成行数为3，列数为\(n\)的服从\(U\left( 0,1 \right)\)的随机数矩阵；

Step 3: 按列判断，如果这一列第一个数字小于\(\frac{1}{2}\)，接受这一列第二个数字数字；

　　　  否则接受这一列第二个数字和第三个数字中较大的那个。

Step 4: 输出：数据列举、算法效率、ks检验；

2.3.3算法效率

定义抽样效率为：一次性抽取随机数中最终被接受数的占所有抽取随机数的比例。

复合法抽样分两步进行，第一步，抽取服从\(U\left( 0,1 \right)\)的随机数，若该数小于\(\frac{1}{2}\)，选取重新抽取的一个服从\(U\left( 0,1 \right)\)的随机数；若该数大于\(\frac{1}{2}\)，选取重新抽取的两个个服从\(U\left( 0,1 \right)\)的随机数中较大的那个。这两种抽样，第一种抽样效率\(\frac{1}{2}\times 1\)，第二种抽样效率\(\frac{1}{2}\times \frac{1}{2}\)，总共\(\frac{3}{4}\)。因此，复合法理论抽样效率\(\frac{3}{4}\)。

三、函数运行

3.1 函数说明

• Fun1( )函数实现逆变换法生成随机数

  输入：n 数据类型：int 含义：需要生成随机数的个数

  输出：列表：数据列举：生成随机数中的前十个数字

  　　　　　   算法效率：该算法的效率

  　　　　　　ks检验：该算法ks检验结果

  　　　图像：抽样密度曲线与理论密度曲线的对比图 　　　 　　　

• Fun2( )函数实现舍选法生成随机数

  输入：n 数据类型：int 含义：需要生成随机数的个数

  输出：列表：数据列举：生成随机数中的前十个数字

  　　　　　   算法效率：该算法的效率

  　　　　　　ks检验：该算法ks检验结果

  　　　图像：抽样密度曲线与理论密度曲线的对比图 　　　 　　　

• Fun3( )函数实现复合法生成随机数

  输入：n 数据类型：int 含义：需要生成随机数的个数

  输出：列表：数据列举：生成随机数中的前十个数字

  　　　　　   算法效率：该算法的效率

  　　　　　　ks检验：该算法ks检验结果

  　　　图像：抽样密度曲线与理论密度曲线的对比图

3.2 运行结果

• Fun1( )函数运行结果

y1 <- Fun1(10000)

## $数据列举
##           [,1]      [,2]      [,3]      [,4]      [,5]
## [1,] 0.2572231 0.2961799 0.9988311 0.6889095 0.7652938
## [2,] 0.6340863 0.8052790 0.9171550 0.7012357 0.9671554
## 
## $算法效率
## [1] 1
## 
## $ks检验
## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  y
## D = 0.010331, p-value = 0.2362
## alternative hypothesis: two-sided

Plot(y1)

图二：逆变换法抽样结果

• Fun2( )函数运行结果

y2 <- Fun2(10000)

## $数据列举
##           [,1]      [,2]      [,3]      [,4]      [,5]
## [1,] 0.4960783 0.7287884 0.9950708 0.7869499 0.8272841
## [2,] 0.7656270 0.9699430 0.3936638 0.2597989 0.3514632
## 
## $算法效率
## [1] 0.6586314
## 
## $ks检验
## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  y
## D = 0.0064549, p-value = 0.799
## alternative hypothesis: two-sided

Plot(y2)

图三：舍选法抽样结果

• Fun3( )函数运行结果

y3 <- Fun3(10000)

## $数据列举
##           [,1]      [,2]      [,3]      [,4]      [,5]
## [1,] 0.1856155 0.4926821 0.5986530 0.4019989 0.4627393
## [2,] 0.7500650 0.2561250 0.7994553 0.5179717 0.6634379
## 
## $算法效率
## [1] 0.7497
## 
## $ks检验
## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  y
## D = 0.0080103, p-value = 0.5425
## alternative hypothesis: two-sided

Plot(y3)

图四：复合法抽样结果

由ks检验结果，三种算法均可以很好的生成给定服从给定概率密度分布的随机数，并且其抽样效率与理论相符，其中抽样效率最高的是逆变换法，达到100%。

四、附录

4.1 Fun1( )函数代码

#Fun1( )逆变换法
Fun1 <- function(n){
  miu <- runif(n)
  y <- sqrt(2*miu+1/4)-1/2
  
  L <- list(
    '数据列举' = matrix(y[1:10],ncol = 5,nrow = 2),
    '算法效率' = n/n,
    'ks检验' = ks.test(y,function(t){1/2*t+1/2*t^2})
  )
  print(L)#输出列表
  return(y)
 
  
}

4.2 Fun2( )函数代码

#Fun2( )舍选法
Fun2 <- function(n){
  
  y <- rep(0,n)
  i <- 0#计数n
  j <- 0#循环次数
  
  while(i<=n){
    x <- runif(1)
    r <- runif(1)
    j <- j+1
    
    if(3/2*r<=1/2+x){
      y[i] <- x
      i <- i+1
    }
    
  }
  
  L <- list(
    '数据列举' = matrix(y[1:10],ncol = 5,nrow = 2),
    '算法效率' = n/j,
    'ks检验' = ks.test(y,function(t){1/2*t+1/2*t^2})
    )
  print(L)#输出列表
  
  return(y)
  
}

4.3 Fun3( )函数代码

#Fun3() 混合法
#混合
Fun3 <- function(n){
  
  data <- runif(3*n)
  data <- matrix(data,nrow = 3,ncol = n)

  y <- apply(data,2,function(x){ifelse(x[1]<1/2,x[2],max(x[2],x[3]))})#筛选
  j <- apply(data,2,function(x){ifelse(x[1]<1/2, 1,0)})#效率
  
  L <- list(
    '数据列举' = matrix(y[1:10],ncol = 5,nrow = 2),
    '算法效率' = (2*n-sum(j))/(2*n),
    'ks检验' = ks.test(y,function(t){1/2*t+1/2*t^2})
  )
  print(L)#输出列表
  
  return(y)
  
}

4.4 Plot( )函数代码

#画图函数
Plot <- function(y){
  ggplot() +
    geom_histogram(#绘制抽得的随机数直方图
      aes(x = y, y = ..density..),
      bins = 30,
      fill = 'cornflowerblue',
      color = 'darkgrey'
    ) +
    geom_line(#绘制抽样密度曲线
      aes(x = y,color = 'red'),
      stat = 'density',
      size = 1.2
    ) +
    geom_line(#绘制理论密度曲线
      aes(x = y,y = 1/2+y,color = 'green'),
      size = 1.2
    )  +
    scale_colour_manual(
      name = '图例',
      values = c('red' = 'red', 'green' = 'green'), 
      labels = c( "理论密度","抽样密度")
    )  +
    labs(#加题目,改x,y轴
      x = 'x轴', 
      y = '概率密度'
    ) +
    theme_light() + 
    theme(legend.position = c(0.1,0.9))
  
}