R语言实验报告.doc

.一、 实验目的1.用 R 生成服从某些具体已知分布的随机变量二、 实验内容在 R 中各种概率函数都有统一的形式，即一套统一的前缀+分布函名： d 表示密度函数（density）； p 表示分布函数（生成相应分布的累积概率密度函数）；q 表示分位数函数，能够返回特定分布的分位（quantile）； r 表示随机函数，生成特定分布的随机数（random）。1、 通过均匀分布随机数生成概率分布随机数的方法称为逆变换法。对于任意随机变量 X，其分布函数为 F，定义其广义逆为：F-(u)=infx;F(x)u若uu (0,1)，则F-（u）和X的分布一样Example 1 如果XExp（1）（服从参数为 1 的指数分布），F（x）=1-e-x。若u=1-e-x并且uu(0,1)，则X=-logUExp(1)则可以解出x=-log(1-u)通过随机数生成产生的分布与本身的指数分布结果相一致R 代码如下：nsim = 104U = 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,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、 某些随机变量可由指数分布生成。若X i Exp(1) 独立同分布的随机变量，那么从X i出发可以得到以下三个标准分布Example 2 生成自由度为 6 的2分布。部分内容来源于网络，有侵权请联系删除！R 代码如下： nsim = 104U=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 from Exp,ylab=,xlab=,ncl=150,col=grey,xlim=c(0,8) d = density(Y)plot(d)3、 一种正态分布随机变量模拟使用 Box-Muller 算法得到 N(0,1)随机变量。这种方法与基于中心极限定理的近似算法相比较，Box-Muller 算法是精确的，它由两个均匀分布产生两个独立的正态分布，其仅有的缺点是必须计算 log，cos，sin。具体解释如下：如果两个随机变量U1, U2独立同分布于u(0,1)，那么由此可以产生两个独立的正态分布X1, X2。Example 3 两次产生相同的 10000 个服从正态分布的随机数，作其中一个的概率直方图，并添加正态分布的密度函数线。nsim = 104 set.seed(22)x1 - rnorm(nsim,mean = 0, sd = 1) set.seed(22)x2 - rnorm(nsim,mean = 0, sd = 1)这里的 x1，x2 数值完全相同。4、 设随机变量 X 的分布列PX=xi=pi, 记p(0)=PX0=0, , p(i)=PXxi=。设 r 是0,1区间上的均匀分布的随机数。当且仅当p(i-1)rp(i)时，令X=xi，则Pp(i-1)r size = 1;p = 0.5 rbinom(10,size,p)1 1 0 0 0 0 1 1 0 1 0接下来生成服从 B（10，0.5）的二项分布 size = 10; p = 0.5 rbinom(5,size,p)1 6 6 5 8 5由此可见，随着实验次数 n 的增大，二项分布越来越接近正态分布。R代码如下size = 1;p = 0.5 rbinom(10,size,p)#生成 5 个服从 B（10，0.5）的二项分布随机数size = 10; p = 0.5rbinom(5,size,p)par(mfrow = c(1,3) p=0.25for(n in c(10,20,50)x - rbinom(100,n,p)hist(x,prob = T,main = paste(n = ,n)xvals = 0:npoints(xvals,dbinom(xvals,n,p),type = h,lwd = 3)par(mfrow = c(1,1)5服从 beta 分布的随机变量的一般算法可以基于 Accept-Reject method 来生成，用的 工具分布为均匀分布 U0,1，假设两个参数都大于 1（通用的 rbeta 函数无此约束）。上 界 M 为 beta 分布密度的最大值，Example 5 对于 =2.7， =6.3 ，最大值 M=2.67，求出 beta 分布R代码Nsim=2500 a=2.7;b=6.3 M=2.67u=runif(Nsim,max=M) y=runif(Nsim) x=yudbeta(y,a,b)这里的 x 服从 beta 分布Accept-Reject algorithm的几个关键要点：1. 仅要求比率f/M，所有算法并不依赖于正则化常数。2. f/Mg的上界不必是紧的，该算法一人有效，即