随机变量生成

这部分包含:

  • 采纳R语言内建函数的办法;
  • 采纳基于变换(transformation)的办法;
  • 采纳承受-回绝(accept-reject)的办法;

基于R内建函数

hist(rgamma(200,2.5,4.5)) #依照gamma(2.5,4.5)生成100个样本点并画出hist直方图

hist(runif(200,min=2,max=5) )#依照在区间[2,5]上的均匀分布生成200个样本点并画出hist直方图

其余一些常见的R语言内建散布如下:

  • 高斯分布:rnorm(n,mean=0,std=1)
  • 指数分布:rexp(n,rate=1)
  • 泊松散布:rpois(n,lambda)
  • Weibull散布:rweibull(n,shape,scale=1)
  • 柯西散布:rcauchy(n,location=0,scale=1)
  • Beta散布:rbeta(n,shape1.shape2)
  • t散布:rt(n,df)
  • F散布:rf(n,df1,df2)
  • 卡方散布:rchisq(n,df)
  • 二项分布:rbinom(n,size,prob)
  • 超几何散布:rhyper(nn,m,n,k)
  • Logistic散布:rlogis(n,location=0,scale=1)
  • 负二项分布:rnbinom(n,size,prob)
  • Wilcoxon散布:rwilcox(nn,m,n)
  • 对数正态分布:rlnorm(n,meanlog=0,sdlog=1)

基于变换(transformation)

基于变换的办法指的是采纳先依照简略的容易生成的散布生成样本点,而后利用某种变换关系将遵从简略散布的样本遵从满足简单散布的样本。比方先采纳R内建函数runif生成100个[0,1]均匀分布的样本点,而后利用变换将其变为遵从较简单的指数分布的样本点。

逆变换法(inverse transformation)

这种办法的流程如下所示:

  1. 取得想要生成的散布的累积散布函数F(x)并推导出其反函数$F^{-1}(x)$;
  2. 利用内建函数生成大量遵从均匀分布U(0,1)的样本点;
  3. 对上一步生成的样本点进行$F^{-1}(x)$变换,其后果能够被认为是遵从想要生成的F(x)散布。

为了了解这种办法的原理,首先给出如下定理:

若X是一个间断随机变量,其累积散布函数为F(x),则$U=F(x)\sim U(0,1)$。

这个定理能够简略了解如下:累积散布函数F(x)是把所有X的取值一对一地映射到[0,1]区间上。离散来看,能够认为F(x)在数轴[0,1]上每个点均有一个样本,因而“线密度”是平均的。如果对于$U=F(x)∈[0,1]$进行抽样,抽中每个值的概率均等,因而U=F(x)遵从U(0,1)。

文献(Introduction to probability and mathematical statistics / Lee J. Bain, Max Engelhardt)对下面定理给出了如下一个比拟严格的证实:

$$\begin{aligned}F_{u}(U)&=P(U \leq u) \\&=P(F(X) \leq u)\\ &=P(F^{-1}(F^{}(X)) \leq F^{-1}(u))\\&=P(X \leq F^{-1}(u))\\&=F(F^{-1}(u))\\&=u\end{aligned}$$

因为$F_{u}(U)∈[0,1]$且$F_{u}(U)=u$,显然$U$遵从U(0,1)

如果间接利用内建函数生成了$u \sim U(0,1)$,则$F^{-1}(u)$和X同散布能够证实如下:

$$\begin{aligned}F_{u}(U)&=P(U \leq u) \\&=P(F^{-1}(X) \leq F^{-1}(x))\\ &=P(F^{}(F^{-1}(X)) \leq F^{}(F^{-1}(x)))\\&=P(X \leq x)\end{aligned}$$

上式中第一个等号是依据累积散布函数的定义,第二个等号是依据上文对于F(x)的定理,第三个等号是因为F(x)是枯燥增函数。

试验:利用逆变换办法生成指数分布

指数分布的CDF为$F(x)=1-e^{-x}$,CDF的反函数为$x=-ln(1-F(x))$。R代码如下:

> N=10000  #样本数量> U=runif(N) #生成的U(0,1)样本> X=-log(1-U) #进行逆变换> hist(X) #画出直方图> hist(rexp(N)) #画出R内建的指数分布直方图

其余基于变换的办法

当$X \sim exp(1)$时,能够依据如下公式生成几种常见的散布:

$$Y=2\sum_{j=1}^{\nu}X_{j} \sim \chi_{2 \nu }^2 \qquad \nu\in\mathbb{N}^{*}\\Y=\beta\sum_{j=1}^{\alpha}X_{j} \sim \Gamma(\alpha,\beta) \qquad \alpha\in\mathbb{N}^{*}\\Y=\frac{\sum_{j=1}^{a}X_{j}}{\sum_{j=1}^{a+b}X_{j}} \sim \beta(a,b) \qquad a,b\in\mathbb{N}^{*}\\$$

Box-Muller算法能够用于生成正态分布,依照以下公式生成的$X_{1}$和$X_{2}$为互相独立的正态分布:

$$U_{1},U_{2} \sim U(0,1)\\X_{1}=\sqrt{-2log(U_{1})}cos(2 \pi U_{2})\\X_{2}=\sqrt{-2log(U_{1})}sin(2 \pi U_{2})\\$$

混合办法(Mxiture representations)

混合办法指的是依照下式生成概率密度散布$f(x)$的办法:

$$f(x)=\int_{\gamma}g(x|y)p(y)dy$$

其中$p(y)$和$g(x|y)$均为容易生成的散布。

应用混合办法时,学生成$p(y)$再生成$g(x|y)$,就可能失去$X \sim f(x)$。

试验:采纳混合办法生成$X\sim \mathcal{N}eg(6,0.3)$

依据相干实践,$X\sim \mathcal{N}eg(n,p)$能够合成为如下两个散布:

$$Y \sim \Gamma(n,(1-p)/p)\\X|y \sim \mathcal{P}(y)$$

依照上式生成X并与R内建函数进行比照的代码如下:

> N=10000> n=6;p=0.3> y=rgamma(N,n,rate=p/(1-p))> x=rpois(N,y)> par(mfcol=c(1,2))> hist(x)> hist(rnbinom(N,n,p))

试验后果如下图所示,能够看出依照混合办法生成的样本点和内建函数生成$X\sim \mathcal{N}eg(6,0.3)$的散布十分靠近。

承受-回绝算法(accept-reject algorithm)

当曾经获取想要生成的散布的概率密度散布$f(x)$(称为指标散布,target density)时,还有一种通过简略散布$g(x)$(称为提议散布,candidate density)生成简单散布的办法:承受-回绝算法。利用这种算法前必须依照如下步骤先确定$g(x)$和常数M:

  • 依据$f(x)$的固有个性,确定一个满足如下条件的$g(x)$:

(1)$g(x)$必须足够容易和不便生成;
(2)当$g(x)>0$时,$f(x)>0$;
(3)存在常数M,使得$M \ge \frac{f(x)}{g(x)}$对任意x均成立;

  • 利用数学原理或者数值优化算法失去符合要求的常数M

当$g(x)$和M曾经确定后,依据如下步骤生成$X \sim f(x)$:

  1. 生成N个提议散布$y\sim g(y)$的样本点$Y_{N}$和均匀分布$u\sim U(0,1)$的样本点$U_{N}$;
  2. 对每个$n\leq N$,当$U_{n} \leq \frac{f(Y_{n})}{Mg(Y_{n})}$时,承受$Y_{n}$作为生成的X样本点,否则回绝$Y_{n}$。

阐明这个算法的有效性等于要证实$f_{y}(y|accepted)=f_{x}(y)$,证实如下:

$$\begin{aligned}f_{y}(y|accepted)&=\frac{f_{y,accetped}(y,accetped)}{P(accepted)} \\&=\frac{f_{y,accetped}(y,accepted)}{\int_{y}{f_{accetped}(accetped|y)}f_{y}(y)}\\ &=\frac{g(y) \frac{f_x(y)}{Mg(y)}}{\int_{y}{\frac{f_{x}(y)}{Mg(y)}g(y)dy}}\\&=\frac{\frac{1}{M}f_{x}(y)}{\frac{1}{M}}\\&=f_{x}(y)\end{aligned}$$

上述证实中$f_{y}(y)=g(y)$,$f_{x}(y)$其实就是想要生成的$f(x)$,第二个等号处利用了贝叶斯定理。

从上述证实中还能够看出:

$$\begin{aligned}P(accepted)&=\int_{y}{f_{accetped}(accetped|y)}f_{y}(y) \\&=\frac{1}{M}\end{aligned}$$

因而,在承受-回绝算法中如果一次性生成了N个提议散布的样本点,最终可能被承受为X的大概有$\frac{N}{M}$个样本点。

试验:承受-回绝算法生成$X \sim \beta(2.7,6.3)$

利用承受-回绝算法生成$X \sim \beta(2.7,6.3)$,采纳$y \sim U(0,1)$作为提议散布,在已知$\beta(2.7,6.3)$的概率密度函数状况下,利用优化算法失去M值为2.67。R语言代码如下:

> M=2.67> N=10000> u=runif(N,max = M)> y=runif(N) #生成提议散布> x=y[u<dbeta(y,2.7,6.3)] #进行承受-回绝,只保留y中合乎承受条件的index的值作为生成x的后果> par(mfcol=c(1,2))> hist(x)> hist(rbeta(10000,2.7,6.3))

试验后果如下图所示,能够看出依照承受-回绝算法生成的样本点和内建函数生成$X \sim \beta(2.7,6.3)$的散布十分靠近。另外,生成的X样本点共有3824个,承受率为38.2%,和实践承受率$\frac{1}{M}$即37%非常靠近。