逆变换法产生特定分布的随机数

 

学校的测试题,因为主要负责这部分,所以就顺带总结一下

学校的测试题,因为主要负责这部分,所以就顺带总结一下

逆变换法产生特定分布的随机数

原理

R中所有随机变量的产生都需要通过均匀伪随机数产生器(uniform random number generator)生成。0-1之间的均匀伪随机数可以使用runif(n)函数实现,生成a到b之间的均匀随机数则可通过runif(n, a, b)。通过均匀分布随机数生成任意概率分布随机数的方法称为逆变换法(Inverse Transform method)。

基于以下定理(Probability Intergral Transformation) 若$X$为连续性随机变量,其$cdf$为$F_X$,则$U=F_X^{-1}(X)$ ~ $U(0,1)$

若$X$为离散型随机变量,其可能的取值记为:

\(\cdots<x_{i-1}<x_{i}<x_{i+1}<\cdots\) 定义

\[F_{X}^{-1}(u)=\inf \left\{x : F_{X}(x) \geq u\right\}\]

则逆变换为$F_{X}^{-1}(u)=x_{i}$,其中$F_{X}\left(x_{i-1}\right)<u \leq F_{X}\left(x_{i}\right)$.从而产生随机数x的方式为:取$x=x_i$,$F_{X}\left(x_{i-1}\right)<u \leq F_{X}\left(x_{i}\right)$

证明

定义:

\[F_{X}^{-1}(u)=\inf \left\{x : F_{X}(x)=u\right\}, \quad 0<u<1\]

那么

\[\begin{aligned} P\left(F_{X}^{-1}(U) \leq x\right) &=P\left(\inf \left\{t : F_{X}(t)=U\right\} \leq x\right) \\ &=P\left(U \leq F_{X}(x)\right)=F_{U}\left(F_{X}(x)\right)=F_{X}(x) \end{aligned}\]

用逆变换法产生Gamma 分布随机数:

\[F(x)=\int_{0}^{x} \frac{\lambda e^{-\lambda y}(\lambda y)^{n-1}}{(n-1) !} d y\]

步骤

Step 1. 产生区间$(0,1)$上的$n$ 个均匀分布随机数$U_1,U_2,\dots,U_n$;

Step 2. 令

\[X=-\frac{1}{\lambda} \log U_{1}-\frac{1}{\lambda} \log U_{2}-\cdots-\frac{1}{\lambda} \log U_{n}=-\frac{1}{\lambda} \log \left(U_{1} U_{2} \cdots U_{n}\right)\]
b.gamma=function(m, n, lambda){
  X=0
  for (i in 1:m) {
    U=runif(n); X[i]=-1/ lambda*sum(log(U))
  }
  X
}
data=b.gamma(100000,3,0.2)
hist(data,prob=TRUE,breaks = 100,main = "Histogram of rGamma")
#hist(rgamma(100,3,0.2))
x = seq(0,max(data),length.out=1000)
curve(dgamma(x,3,0.2), col="red",add=TRUE)

效果还是可以的:

Histogram of rGamma

用逆变换法产生离散型随机变量随机数

在泊松分布分布律给定的情形下,运用R软件,利用逆变换法模拟产生随机数。泊松分布分布律如下:

\[p_{i}=P\{X=i\}=e^{-\lambda} \frac{\lambda^{i}}{i !}, \quad i=0,1, \cdots\]

步骤

Step 1. 产生一个区间$(0,1)$上的均匀分布随机数$U$;

Step 2. 令$i = 0$, $p = e^{-\lambda}$ , $F =p$;

Step 3. 如果$U < F$ ,令$X = i$并停止;

Step 4. $p \leftarrow \lambda p /(i+1), F \leftarrow F+p, i=i+1$

Step 5. 返回Step 3.

b.rpois=function(n,lambda){
  Y=0
  for(j in 1:n) {
    u=runif(1)
    i=0;p=exp(-lambda);F=p
    while (u>=F){
      p=lambda*p/(i+1);
      F=F+p; 
      i=i+1;
    }
    Y[j]=i;
  }
  Y
}
data=b.rpois(100,3)
hist(data,prob=TRUE)

用极坐标变换法(Box-Muller)产生正态分布随机数

正态分布密度函数:

\[f(x)=\frac{1}{\sqrt{2 \pi} \sigma} e^{-\frac{(x-\mu)^{2}}{2 \sigma^{2}}}\]

原理

当$x$和$y$是两个独立且服从[0,1)均匀分布的随机变量时,$\cos (2\pi x) \cdot \sqrt { - 2\ln (1 - y)}$和$\sin (2\pi x) \cdot \sqrt { - 2\ln (1 - y)} $是两个独立且服从正态分布的随机变量

本质上也是应用了反变换法

把反变换法推广到二维的情况,设$U_1,U_2$为$(0,1)$上的均匀分布随机变量,$(U_1,U_2)$的联合概率密度函数为$f(u_1,u_2)=1(0 \le u_1,u_2 \le 1)$,若有:

\[\left\{\begin{matrix} {U_1} = {g_1}(X,Y)\\ {U_2} = {g_2}(X,Y) \end{matrix}\right.\]

其中,$g_1,g_2$的逆变换存在,记为

\[\left\{\begin{matrix} x = {h_1}({u_1},{u_2})\\ y = {h_2}({u_1},{u_2}) \end{matrix}\right.\]

且存在一阶偏导数,设$J$为Jacobian矩阵的行列式

\[\left | J \right | =\begin{vmatrix} \frac{\partial x}{\partial u_1} & \frac{\partial x}{\partial u_2} \\ \frac{\partial y}{\partial u_1} & \frac{\partial y}{\partial u_2} \end{vmatrix} \ne 0\]

则随机变量$(X,Y)$的二维联合密度为(回顾直角坐标和极坐标变换):

\[f[h_1(u_1,u_2),h_2(u_1,u_2)] \cdot \left | J \right | = \left | J \right |\]

根据这个定理我们来证明一下,

\[\left\{\begin{matrix} X = \sqrt {- 2\ln U_1} \cos (2\pi U_1) \\ Y = \sqrt {- 2\ln U_1} \sin (2\pi U_2) \end{matrix}\right.\]

求反函数得

\[\left\{\begin{matrix} U_1 = e^{ - \frac{X^2 + Y^2}{2}} \\ U_2 = \frac{1}{2 \pi} \arctan \frac{Y}{X} \end{matrix}\right.\]

计算Jacobian行列式

\[\begin{align} \left | J \right | =\begin{vmatrix} \frac{\partial U_1}{\partial X} & \frac{\partial U_1}{\partial Y} \\ \frac{\partial U_2}{\partial X} & \frac{\partial U_2}{\partial Y} \end{vmatrix} & =-(\frac{1}{\sqrt{2 \pi}}e^{-\frac{1}{2}X^2})(\frac{1}{\sqrt{2 \pi}}e^{-\frac{1}{2}Y^2}) \end{align}\]

由于$X_1,X_2$为$(0,1)$上的均匀分布,概率密度函数均为$1$,所以$X,Y$的联合概率密度函数为

\[-(\frac{1}{\sqrt{2 \pi}}e^{-\frac{1}{2}X^2})(\frac{1}{\sqrt{2 \pi}}e^{-\frac{1}{2}Y^2})\]

所以$Y_1,Y_2$是两个独立且服从正态分布的随机变量

步骤

Step 1. 产生区间$(0, 1)$ 上的均匀分布随机数$U_1,U_2$

Step 2. 令$R^{2}=-2 \log \left(U_{1}\right), \quad \theta=2 \pi U_{2}$

Step 3. 令

\[X=R \cos \theta=\sqrt{-2 \log \left(U_{1}\right)} \cos \left(2 \pi U_{2}\right)\] \[Y=R \sin \theta=\sqrt{-2 \log \left(U_{1}\right)} \sin \left(2 \pi U_{2}\right)\]
b.rnorm=function(n){
  X=0; Y=0
  for (i in 1:n) {
    U1=runif(1); U2=runif(1)
    X[i] <- sqrt(-2*log(U1))*cos(2*pi*U2)
    Y[i] <- sqrt(-2*log(U1))*sin(2*pi*U2)
  }
  c(X, Y)
}

hist(b.rnorm(1000),prob=TRUE)
curve(dnorm, col="red",add=TRUE)