本文旨在介绍如何生成满足特定分布的随机数,并给出了每类方法的大致原理和实例
方法一,使用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))
其样本协方差矩阵如下,非常接近总体。