乐趣区

关于统计:Introducing-Monte-Carlo-Methods-with-R第二章笔记

jiezi

随机变量生成

这部分包含：

  • 采纳 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% 非常靠近。

退出移动版