R语言实验五

实验五常见分布的相关计算、随机抽样与模拟【实验类型】验证性【实验学时】2学时【实验目的】1、掌握常见分布的分布函数、密度函数（或分布列）及分位数的计算方法；2、掌握样本统计量的计算方法及所表达的意义；3、了解随机模拟的基本思想及其应用。【实验内容】1、组合数与组合方案的生成、概率的计算，2、常见分布的分布函数、密度函数（或分布列）以及分位数的计算；3、随机数的生成与随机模拟（蒙特卡洛仿真）。【实验方法或步骤】第一部分、课件例题：1.#从1～5个数中，随机取3个的全部组合combn(1:5,3)#共10种组合方案combn(1:5,3,FUN=mean)#对每种组合方案求均值2.choose(5,3)#从5个数里面选3个的组合数目choose(50,3)factorial(10)#计算10！3.#3.从一副完全打乱的52张扑克中任取4张，计算下列事件的概率#(1)抽取4张依次为红心A，方块A，黑桃A和梅花A的概率1/prod(49:52)#prod()表示连乘积#(2)抽取4张为红心A，方块A，黑桃A和梅花A的概率．1/choose(52,4)4.设在15只同类型的零件中有2只是次品，一次任取3只，以X表示次品的只数，求X的分布律．x-c(1,1,rep(0,13));x#样本空间(用1表示次品,0为正品)X-combn(x,3,FUN=sum)#从样本空间中任取3个元素的方案，并对每个方案求和，共455个数（取值0,1,2）p-numeric(3)#定义p为数值型的3维向量，且初值为0for(iin1:3)p[i]-sum(X==i-1)/length(X)#sum(X==i-1)表示对X取值为i-1的个数求和,X的长度为455p5.#例5.3：计算3σ原则对应的概率x-1:3;p-pnorm(x)-pnorm(-x);p#例5.4：令α=0.025，计算上α分位点zα．alpha-0.025;z-qnorm(1-alpha);z6.#例5.5：计算P{X≤160}，其中X~U[150,200]。p-punif(160,min=150,max=200);p#例5.6：已知X~E(1/241)，计算P{50X≤100}。p-pexp(100,rate=1/241)-pexp(50,rate=1/241);p#例5.7：已知X~B(80,0.01)，计算P{X≥4}。p-1-pbinom(3,size=80,prob=0.01);p#例5.8：已知X~B(180,0.01)，且P{X≤k}≥0.95，求k。n-180;p-0.01;k-qbinom(0.95,n,p);k#例5.9：已知X~B(1000,0.01)，求P{X≥2}。p1-1-pbinom(1,size=1000,prob=0.001);p1p2-1-ppois(1,lambda=1000*0.001);p27.x-c(0:10,50);xxm-mean(x);xmc(xm,mean(x,trim=0.10))#去掉两端各10%数据后取平均8.x-c(12,9,11,5,1,4,8,3,2,10,6,7);xy-1:12;yvar(x)#方差sd(x)#标准差var(x,y)#协方差cov(x,y)#协方差9.x-c(12,9,11,5,1,4,8,3,2,10,6,7,NA);xsort(x)#默认升序，不处理缺失数据sort(x,decreasing=TRUE)#降序sort(x,na.last=TRUE)#处理缺失数据sort(x,decreasing=TRUE,na.last=FALSE)10.#10.x-c(12,9,11,5,1,4,8,3,2,10,6,7,NA)median(x)median(x,na.rm=TRUE)11.分位数quantile(1:10)12.#用于求解样本k阶中心矩的程序moment.R为：moment-function(x,k,mean=0)sum((x-mean)^k)/length(x)skew-function(x,flag=1){#计算偏度系数的程序mu-mean(x)if(flag==1){m2-moment(x,2,mean=mu)m3-moment(x,3,mean=mu)Cs-m3/sqrt(m2^3)}else{n-length(x);S-sd(x)Cs-n/(n-1)/(n-2)*sum((x-mu)^3)/S^3}Cs}x-rnorm(10000)#服从正态分布的随机数skew(x)#默认使用第一种方法(有偏)计算偏度系数skew(x,2)#使用第二种方法(无偏)计算偏度系数13.x-rnorm(12)#服从正态分布的随机数Fn-ecdf(x)#经验分布函数op-par(mfrow=c(3,1),mgp=c(1.5,0.8,0),mar=0.1+c(3,3,2,1))#绘制多图，3行1列plot(Fn)#默认参数，不画竖线plot(Fn,verticals=TRUE)#有竖线plot(Fn,verticals=TRUE,do.points=FALSE)#不画点par(op)#多组图14.生成随机数x-rnorm(1000)par(mai=c(0.9,0.9,0.6,0.2))#图形边界参数hist(x,prob=TRUE,col=lightblue,#直方图main=NormalDistribution)curve(dnorm(x),add=TRUE,col=red,lwd=2)#绘制正态密度曲线expr-expression(paste(mu==0,,,sigma==1))legend(1,0.35,legend=expr,col=red,lwd=2)#添加给定内容的图例15.生成10个随机数，且每次的运算的结果均相同。set.seed(166)runif(10)set.seed(166)runif(10)16.bootstrap.ci-function(x,B,alpha=0.05){medians-numeric(B)#生成一个长度为B的向量for(iin1:B){sam-sample(x,replace=TRUE)#有放回抽样medians[i]-median(sam)#向量每个元素为中位数}quantile(medians,c(alpha/2,1-alpha/2))#百分位数函数，计算置信水平为1-α的置信区间}#例5.11：：用样本中位数计算总体中位数θ的置信区间x-c(136.3,136.6,135.8,135.4,134.7,135.0,134.1,143.3,147.8,148.8,134.8,135.2,134.9,146.5,141.2,135.4,134.8,135.8,135.0,133.7,134.4,134.9,134.8,134.5,134.3,135.2)#样本bootstrap.ci(x,B=1000,alpha=0.05)#置信水平为0.95例5.12的求解：train-function(n){t1-sample(c(0,5,10),size=n,replace=T,prob=c(0.7,0.2,0.1))#火车出发时刻t2-rnorm(n,mean=30,sd=2)#从A站到B站的运行时间t3-sample(c(28,30,32,34),size=n,replace=T,prob=c(0.3,0.4,0.2,0.1))#此人到达B站时刻sum(t1+t2t3)/n#成功(赶上火车)出现的频率}train(10000)#模拟10000次试验并计算概率#求解定积分的quad1()函数quad1-function(fun,a,b,M=1,n=1e6){#M为函数上界x-runif(n,min=a,max=b)y-runif(n,min=0,max=M)#生成矩形区域的随机数k-sum(yfun(x))#落在曲边梯形内的随机点个数k/n*M*(b-a)}#k/n表示梯形面积与矩形面积的比值#例5.13fun-function(x)exp(-x^2)#被积函数quad1(fun,a=-1,b=1)#计算定积分integrate(fun,-1,1)#使用R中的积分函数#求解二重积分的quad2()函数quad2-function(fun,a,b,c,d,M=1,n=1e6){#M为上界x-runif(n,min=a,max=b)y-runif(n,min=c,max=d)z-runif(n,min=0,max=M)#生成长方体内的随机数对k-sum(zfun(x,y))#落在曲顶柱体内的随机点个数k/n*M*(b-a)*(d-c)#k/n表示曲顶柱体占长方体体积比值}#例5.14fun-function(x,y)sqrt((x^2+y^2=1)*(1-x^2-y^2))#定义被积函数，定义域内函数值不变，定义域外函数值为0quad2(fun,a=-1,b=1,c=-1,d=1)#计算二重积分RandomWalk-function(n=10000){#编制函数par(mai=c(0.9,0.9,0.5,0.2))#图形参数plot(c(0,100),c(0,100),type=n,xlab=X,ylab=Y,main=RandomWalkinTwoDimensions)x-y-50#两变量同时赋值points(50,50,pch=16,col=red,cex=1.5)#绘制初始点for(iin1:n){xi-sample(c(1,0,-1),1)#从1,0,-1中抽样一次yi-sample(c(1,0,-1),1)lines(c(x,x+xi),c(y,y+yi))#绘制步骤的连线x-x+xi;y-y+yi#下一步行走位置if(x100|x0|y100|y0)break}}#判断结束条件RandomWalk()#调用函数第二部分、教材例题：1.#例5.1.2X-c(1,1,0,1,0,0,1,0,1,1,1,0,1,1,0,1,0,0,1,0,1,0,1,0,0,1,1,0,1,1,0,1)theta-mean(X)t-theta/(1-theta)t#2#例5.1.4X-c(0.59132754,0.12854935,0.46900228,0.29835980,0.24341462,0.06566637,0.40085536,2.99687123,0.05278912,0.09898594)lambda-1/mean(X)lambdalambda-1/sd(X)#使用二阶矩进行矩法估计lambda#3#例5.1.5f-function(P)(P^517)*(1-P)^483optimize(f,c(0,1),maximum=TRUE)#4#z.test()函数定义z.test-function(x,n,sigma,alpha,u0=0,alternative=two.sided){options(digits=4)result-list()mean-mean(x)z-(mean-u0)/(sigma/sqrt(n))p-pnorm(z,lower.tail=FALSE)result$mean-meanresult$z-zresult$p.value-pif(alternative==two.sided){p-2*presult$p.value-p}elseif(alternative==greater|alternative==less){result$p.value-p}elsereturn(yourinputiswrong)result$conf.int-c(mean-sigma*qnorm(1-alpha/2,mean=0,sd=1,lower.tail=TRUE)/sqr