R语言实验报告6

1、一、实验目的用R生成服从某些具体已知分布的随机变量二、实验内容在R中各种概率函数都有统一的形式，即一套统一的前缀 +分布函名: d表示密度函数(density)；p表示分布函数(生成相应分布的累积概率密度函数)；q表示分位数函数，能够返回特定分布的分位(quan tile)；r表示随机函数，生成特定分布的随机数(random)。R中的各种概率统计分布汉文名称英文名称R対应的名字附加参数B分布betabetashapel, shape2r nep二项式分布binomialbinomsize, prob柯西分布Cauchycauchylocation, scale卡方分布ch卜squaredchi

2、sqdf, nep指数分布exp on entialexprateF分布Ffdfl, dflz nepGamma(Y)分布gammagammashape, scale几何分布geometricgeomprob超几何分布hypergeometrichypermz n, k对数正态分布log no rmalInormmean log, sdlogLogistic 分布logisticlogislocation, scale负二项式分布negative binomialnbinomsize, prob正态分布nor malnormmean, sd泊松分布Poiss onpoislambdaWilco

3、xon 分布signed ranksignrankn|t分布Students ttdfz nep均匀分布uniformunifmin, max韦们分布Weibullweibullshape, scale:秩和分布:Wilcox onwilcoxm, n1、通过均匀分布随机数生成概率分布随机数的方法称为逆变换法。对于任意随机变量X,其分布函数为F,定义其广义逆为：F-(u)=infx;F(x)Mu若uu (0,1),则F- (u)和x 的分布一样Examplel如果XExp (1)(服从参数为1的指数分布)，F (x) =1-e-x。若u=1-e-x并且 uu(0,1)，则 X=-logUExp

4、(1)则可以解出x=-log(1-u)Exp from RExp from UniformExp from R通过随机数生成产生的分布与本身的指数分布结果相一致R代码如下:nsim = 10A4U = runif(nsim)X = -log(U)Y = rexp(nsim)X11(h=3.5)Xpar(mfrow=c(1,2),mar=c(2,2,2,2)hist(X,freq=F,main=Exp from Uniform,ylab=,xlab=,ncl=150,col=grey,xlim=c(0,8) curve(dexp(x),add=T,col=sienna,lwd=2)hist(Y,

5、freq=F,main=Exp from R,ylab=,xlab=,ncl=150,col=grey,xlim=c(0,8) curve(dexp(x),add=T,col=sienna,lwd=2)2、某些随机变量可由指数分布生成。若XiExp(1)独立同分布的随机变量，那么从Xi出 发可以得到以下三个标准分布aJ=1Example 2生成自由度为6的XaJ=1Example 2生成自由度为6的X分布。chisq from Expooin o og ochisq from Expooin o og odensity.default(x = Y)0102030R代码如下：nsim = 10A

6、4U=matrix(data=U,nrow=3)X=-log(U)X=2* apply(X,2,sum) Y=rchisq(nsim,df = 6)par(mfrow=c(1,2),mar=c(2,2,2,2)hist(X,freq=F,main=chisq fromExp,ylab=,xlab=,ncl=150,col=grey,xlim=c(0,8) d = density(Y) plot(d)3、一种正态分布随机变量模拟使用Box-Muller算法得到N(0,1)随机变量。这种方法与基于中心极限定理的近似算法相比较，Box-Muller算法是精确的，它由两个均匀分布产生两个独立的正态分布

7、，其仅有的缺点是必须计算log，cos，sin。具体解释如下：如果两个 随机变量U1, U2独立同分布于u(0,1),那么由此可以产生两个独立的正态分布, x2。Xi = y-21og(t7i) cos(27rt72) ,= /-21og(t/i) sin(27rL72)Example 3两次产生相同的10000个服从正态分布的随机数，作其中一个的概率直方 图，并添加正态分布的密度函数线。normal mu=0,sigma=1x1nsim = 10A4set.seed(22)x1 - rnorm(nsim,mean = 0, sd = 1) set seed(22)x2 - rnorm(nsi

8、m,mean = 0, sd = 1),x24、设随机变量X 的分布列 PX=x.=p.,记 p（o）=PXWO=O,p（i）=PXWx.=工 i p。1 iij=i i设r是0,1区间上的均匀分布的随机数。当且仅当p（i-i）rp（i）时，令X=xi，则Pp（i- 1）rsize=1;p=0.5rbinom(10,size,p)11000011010接下来生成服从B (10, 0.5)的二项分布size=10;p=0.5rbinom(5,size,p)166585由此可见，随着实验次数n的增大，二项分布越来越接近正态分布。R代码如下size = 1;p = 0.5 rbi nom(10,si

9、ze,p)#生成5个服从B(10，0.5)的二项分布随机数 size = 10;p = 0.5rbin om(5,size,p)par(mfrow = C(1,3)p=0.25for (n in c(10,20,50)x - rbin om(100, n,p) hist(x,prob = T,mai n = paste( n = ”,n) xvals = 0:npoin ts(xvals,db in om(xvals ,n, p),type = h,lwd = 3) par(mfrow = c(1,1)5.服从beta分布的随机变量的一般算法可以基于Accept-Reject method来生成， 用的工具分布为均勻分布U0,1，假设两个参数都大于 1 (通用的rbeta函数无 此约束)。上界M为beta分布密度的最大值，Example5 对于a=2.7, |3=6.3 ，最大值M=2.67，求出beta 分布R代码Nsim=2500 a=2.7;b=6.3M=2.67u=r un if(Nsim,max=M) y=r un if(Nsim) x=yuvdbeta(y,a,b) 这里的x服从beta分布Accept-Rejectalgorithm 的几个关键要点：仅要求比率f/M,所有