【随机数生成方法】各类随机数生成法及其R语言实现

本文旨在介绍如何生成满足特定分布的随机数，并给出了每类方法的大致原理和实例

方法一，使用R语言概率函数

比较常见的发布如正态分布，t、F分布、二项分布等R都提供了直接的概率函数，下表展示了R中可用的单变量概率函数。
其中在分布名前加p累积分布函数（cdf），加d代表概率密度函数，加q表示分位数函数，加r表示抽样函数。

例如正态分布，相应操作如下

rnorm(100,0,1) #抽取100个服从N(0,1)的随机数
pnorm(0,0,1) #输出0.5
qnorm(0.5,0,1) # 输出0
dnorm(0,0,1) # 输出0.3989423

方法二，逆变化法

对于某些比较特殊的分布，R中可能没有提供现成的概率函数，例如双参数指数函数。此时可使用逆变换法生成符合该分布的随机数。逆变换法的基本原理是：对于连续性随机变量X，其分布函数F(x)也是一个随机变量，服从参数为0，1的均匀分布U(0,1)，则依据概率论知识，若u~U(0,1)，则F-1(u) ~F(x),F-1表示F的反函数。所以，逆变换法概括如下：

1.生成随机变量u~U(0,1)

2、计算F-1(u)

3、令x = F-1(u) ,得到服从目标分布的随机数

实例演示

用逆变换法生成1000个服从参数为1的指数分布随机数，并于R语言内置函数比较。

 n <- 1000
 k <- 1
 x <- numeric(n)
 while (k <= n)
 {
   u <- runif(1,0,1) # 步骤1
   x[k] <- -log(1-u) / 1 #步骤2，3
   k <- k + 1
 }
 print(x)
 plot(quantile(x,seq(0,0.8,0.01)),qexp(seq(0,0.8,0.01)),
 xlab="逆变换法分位数",
 ylab= "实际分位数",main = "QQPlot",pch=16)
 abline(a=0,b=1)
 

通过逆变换法生成的随机数和实际分布随机数高度一致，说明方法是有效的。
离散发布下的随机数生成： 以上的逆变换法涉及到了求分布函数的反函数，此时自然衍生出一个问题，如果离散状况下分布函数不连续，则无法求得其反函数。回忆连续情形下的逆变换法的步骤2，其本质是求X的u分位数。分位数的基本意义见下，所以离散情形下的逆变换算法为：

1、生成随机变量u~U(0,1)

2、从支撑集中解出u分位数x，则x即为满足该分布的随机数

例如逆变换法生成100个服从参数为0.4的两点分布

n <- 100
p <- 0.4
u <- runif(n)
x <- as.integer(u>0.6)

方法三，接受拒绝法

在逆变换法中需要求累积分布函数的反分位数函数，而当随机变量的分布函数表达式比较复杂时很难求得其反函数，此时可以使用随机模拟的方式生成随机数。接受拒绝法本质是蒙特卡洛模拟，比如我们计算一个圆的面积可以随机生成在圆外接正方形内的点，当随机点的数量足够多时，圆的面积所占比例就等于在圆内随机点的比例。那么与其类似，如果把圆看出密度曲线，其基本原理就很清晰了，下图中红色的点落入了分布曲线下方，表示接受该点服从F(x),灰色的点在密度曲线外，表示拒绝。这是一个关于接受拒绝法的直观理解。不难看出u是产生于矩形内的随机点，其服从均匀分布，u被称为辅助分布，又称proposal分布。当然proposal分布并不一定要是均匀分布，此时形状将不是矩形。

接受拒绝法的步骤如下：

1、找到一个随机变量Y，其密度函数g(t)满足:f(t)/g(t)<c

2、生成一个随机变量服从U(0,1)

3、若u < f(y)/(cg(y)),则令x=y，否则返回2

例如生成100个参数为（2，2）的贝塔分布，取proposal分布为U(0,1),c=6,此时f(x)/cg(x) = x(1-x)

n <- 100
k <- 0
y <- numeric(n)
while (k < n){
u <- runif(1)
x <- runif(1)  # 这里不能写成x<-u
if (x* (1-x) > u) {
k <- k + 1
y[k] <- x
}
}

方法四，卷积与混合

这种方式一般用于模拟不常见的特定分布，例如多峰正态分布、混合伽马分布，使用场景有限，后续如有需要再补充。

方法五，生成多元分布

以上特定分布随机数生成方法均是单变量分布，若要生成多元分布要借助谱分解、Choleski分解或奇异值分解。下面演示生成生成100个均值向量为0，协方差矩阵如下的二元正态分布

sigma <- matrix(c(1,0.9,0.9,1),nc=2)
n <- 100; u <- 0
z <- matrix(rnorm(200),nc=2,byrow = T)
a <- chol(sigma) # Choleski分解
x <- z %*% a + matrix(rep(u,2*n),nc=2)
print(x)
print(cor(x))

其样本协方差矩阵如下，非常接近总体。