统计

正太散布参考资料 定义一维正太散布规范正太散布 二项分布参考资料 伯努利散布参考资料 均匀分布参考资料 beta散布

蒙特卡罗积分引言motivation数值办法求积分在统计计算中具备非常重要的位置，很多需要求边缘概率密度散布的状况（尤其是算贝叶斯后验）须要用到数值积分。 R中一般一维数值积分办法罕用的数值积分办法在不论是本科还是研究生的数值剖析中都进行了具体的介绍，包含基于梯形公式，辛普森公式，基于Romberg积分等办法，这篇文章中临时不进行具体介绍。值得一提的是，R中实现了两种一维的数值积分法：area（from MASS）和integrate。 试验1 利用R内置函数计算一维数值积分待求的积分如下： $$f(x)=\int\_{0}^{\infty}x^{\lambda-1}exp(-x)dx$$ 采纳integrate函数进行数值计算并和利用gamma函数求得的解析解进行比照，代码如下： > ch=function(la){integrate(function(x){x^(la-1)*exp(-x)},0,Inf)$val}> plot(lgamma(seq(0.01,10,le=50)),log(apply(as.matrix(seq(0.01,10,le=50)),1,ch)),xlab = "log(integrate(f))",ylab = expression(log(Gamma(lambda))))失去的后果如下图所示，能够看出在$\lambda \in [0.01,10]$的范畴内integrate函数求得数值积分和解析解并无显著差别。 第二个要进行数值积分的函数如下式所示，这个函数能够简略认为是当先验为均匀分布时10个柯西散布样本对模型参数的边缘后验（未进行归一化）。 $$m(x)=\int\_{-\infty}^{\infty}\sum^{10}*{i=1} \frac{1}{\pi} \frac{1}{1+(x*{i}-\theta)} dx$$ 设$\theta =350$，代码如下，从运行论断来看显然integrate函数算的有错，因为把积分下限从正无穷改为400后后果居然变小了。 > cac=rcauchy(10)+350> like=function(theta){+ u=dcauchy(cac[1]-theta)+ for(i in 2:10)+ u=u*dcauchy(cac[i]-theta)+ return(u)+ }> integrate(like,-Inf,Inf)5.559304e-44 with absolute error < 1.1e-43> integrate(like,-Inf,400)5.882426e-31 with absolute error < 1.1e-30从新生成10个样本，这次设置$\theta =0$，而后将integrate函数的后果和area函数的后果比照，代码如下。 > library(MASS)> cac=rcauchy(10)> nin=function(a){integrate(like,-a,a)$val}> nan=function(a){area(like,-a,a)}> x=seq(1,1000,le=10000)> y=log(apply(as.matrix(x),1,nin))> z=log(apply(as.matrix(x),1,nan))> plot(x,y,ylim = range(cbind(y,z)),lwd=2)> lines(x,z,col='sienna',lwd=6)试验后果如图所示，图中彩色线为积分限-integrate函数积分值的曲线，很显然谬误，黄色线为积分限-area函数积分值的曲线，其趋势根本正确。 经典蒙特卡罗积分从上文的试验能够看出，当积分的模式变得比较复杂时R内建的一维数值积分算法可能会不稳固。这种惯例数值积分办法在计算高维积分时也面临效率和精度问题。因而，面对统计推断中随处可见的高维简单积分，蒙特卡罗办法被用于解决这个问题。利用蒙特卡罗法进行数值积分的原理能够写成如下公式，其中$\mathbb{E}_{f}[h(X)]$示意蒙特卡罗生成的$X \sim f(x)$的大量样本的$h(x)$均值的实践表达式，将样本均值$\hbar_{n}$作为$\mathbb{E}_{f}[h(x)]$的估计量，就能够间接计算出在$ f(x)$撑持集上对$h(x)f(x)$进行积分的后果。 $$\begin{aligned}\mathbb{E}_{f}[h(X)]&=\int_{\mathcal{X}}^{}h(x)f(x)dx\hbar_{n} &=\frac{1}{n}\sum^{n}_{j=1}h(x_{j})\end{aligned}$$ 在利用场景中，最简略间接的办法就是将被积函数作为$h(x)$，将$f(x)$设为积分区间上的均匀分布$U(a,b)$，生成大量的$X \sim U(a,b)$样本点并用$\hbar_{n}$作为积分后果的估计值。 试验2 利用经典蒙特卡罗办法计算数值积分待求的积分如下： ...

随机变量生成这部分包含： 采纳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）这种办法的流程如下所示： 取得想要生成的散布的累积散布函数F(x)并推导出其反函数$F^{-1}(x)$；利用内建函数生成大量遵从均匀分布U(0,1)的样本点；对上一步生成的样本点进行$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}^{*}\\$$ ...