R语言实验五.docx
- 文档编号:8870705
- 上传时间:2023-05-15
- 格式:DOCX
- 页数:26
- 大小:513.98KB
R语言实验五.docx
《R语言实验五.docx》由会员分享,可在线阅读,更多相关《R语言实验五.docx(26页珍藏版)》请在冰点文库上搜索。
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维向量,且初值为0
for(iin1:
3)
p[i]<-sum(X==i-1)/length(X)#sum(X==i-1)表示对X取值为i-1的个数求和,X的长度为455
p
5.
#例5.3:
计算3σ原则对应的概率
x<-1:
3;p<-pnorm(x)-pnorm(-x);p
#例5.4:
令α=0.025,计算上α分位点zα.
alpha<-0.025;z<-qnorm(1-alpha);z
6.
#例5.5:
计算P{X≤160},其中X~U[150,200]。
p<-punif(160,min=150,max=200);p
#例5.6:
已知X~E(1/241),计算P{50 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);p1 p2<-1-ppois(1,lambda=1000*0.001);p2 7. x<-c(0: 10,50);x xm<-mean(x);xm c(xm,mean(x,trim=0.10))#去掉两端各10%数据后取平均 8. x<-c(12,9,11,5,1,4,8,3,2,10,6,7);x y<-1: 12;y var(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);x sort(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+t2>t3)/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(y k/n*M*(b-a)}#k/n表示梯形面积与矩形面积的比值 #例5.13 fun<-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(z k/n*M*(b-a)*(d-c)#k/n表示曲顶柱体占长方体体积比值 } #例5.14 fun<-function(x,y)sqrt((x^2+y^2<=1)*(1-x^2-y^2))#定义被积函数,定义域内函数值不变,定义域外函数值为0 quad2(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(x>100|x<0|y>100|y<0)break}}#判断结束条件 RandomWalk()#调用函数 第二部分、教材例题: 1. #例5.1.2 X<-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.4 X<-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) lambda lambda<-1/sd(X)#使用二阶矩进行矩法估计 lambda #3 #例5.1.5 f<-function(P)(P^517)*(1-P)^483 optimize(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<-mean result$z<-z result$p.value<-p if(alternative=="two.sided"){ p<-2*p result$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)/sqrt(n), mean+sigma*qnorm(1-alpha/2,mean=0,sd=1, lower.tail=TRUE)/sqrt(n)) result } #求置信区间的程序: conf.int<-function(x,n,sigma,alpha){ options(digits=4) mean<-mean(x) c(mean-sigma*qnorm(1-alpha/2,mean=0,sd=1, lower.tail=TRUE)/sqrt(n), mean+sigma*qnorm(1-alpha/2,mean=0,sd=1, lower.tail=TRUE)/sqrt(n)) } #例5.2.1 x<-c(175,176,173,175,174,173,173,176,173,179) result<-z.test(x,10,1.5,0.05) result$conf.int x<-c(175,176,173,175,174,173,173,176,173,179) conf.int(x,10,1.5,0.05) #5 #不知道方差,用函数t.test()来求置信区间 x<-c(175,176,173,175,174,173,173,176,173,179) t.test(x) t.test(x)$conf.int#提取出置信区间的部分 #6 #函数chisq.var.test()可以用来求σ2置信区间 chisq.var.test<-function(x,var,alpha,alternative="two.sided"){ options(digits=4) result<-list() n<-length(x) v<-var(x) result$var<-v chi2<-(n-1)*v/var result$chi2<-chi2 p<-pchisq(chi2,n-1) if(alternative=="less"|alternative=="greater"){ result$p.value<-p }elseif(alternative=="two.sided"){ if(p>.5) p<-1-p p<-2*p result$p.value<-p }elsereturn("yourinputiswrong") result$conf.int<-c( (n-1)*v/qchisq(alpha/2,df=n-1,lower.tail=F), (n-1)*v/qchisq(alpha/2,df=n-1,lower.tail=T)) result } #7 #编写函数求置信区间 two.sample.ci<-function(x,y,conf.level=0.95,sigma1,sigma2){ options(digits=4) m=length(x);n=length(y) xbar=mean(x)-mean(y);alpha=1-conf.level zstar=qnorm(1-alpha/2)*(sigma1/m+sigma2/n)^(1/2) xbar+c(-zstar,+zstar) } #例5.3.1 x<-c(628,583,510,554,612,523,530,615) y<-c(535,433,398,470,567,480,498,560,503,426) sigma1<-2140 sigma2<-3250 two.sample.ci(x,y,conf.level=0.95,sigma1,sigma2) #8 #例5.3.2 x<-c(628,583,510,554,612,523,530,615) y<-c(535,433,398,470,567,480,498,560,503,426) t.test(x,y,var.equal=TRUE) #9 #例5.3.3 x<-c(20.5,19.8,19.7,20.4,20.1,20.0,19.0,19.9) y<-c(20.7,19.8,19.5,20.8,20.4,19.6,20.2) var.test(x,y) #10 #例5.4.1 prop.test(38,200,correct=TRUE)#函数prop.test()对p进行估计与检验 binom.test(38,200)#函数binom.test()可以求其置信区间 #11 #例5.5.1 like<-c(478,246) people<-c(1000,750) prop.test(like,people) #自己编写没有修正的两比例之间的区间估计函数ratio.ci() ratio.ci<-function(x,y,n1,n2,conf.level=0.95){ xbar1=x/n1;xbar2=y/n2 xbar=xbar1-xbar2 alpha=1-conf.level zstar=qnorm(1-alpha/2)*(xbar1*(1-xbar1)/n1+xbar2*(1-xbar2)/n2)^(1/2) xbar+c(-zstar,+zstar) } ratio.ci(478,246,1000,750,conf.level=0.95) #12 #例5.6.1 #定义函数size.norm1()求样本容量 size.norm1<-function(d,var,conf.level){ alpha=1-conf.level ((qnorm(1-alpha/2)*var^(1/2))/d)^2 } size.norm1(2,500,conf.level=0.95) #13 #通过循环确定样本容量,定义函数size.norm2() size.norm2<-function(s,alpha,d,m){ t0<-qt(alpha/2,m,lower.tail=FALSE) n0<-(t0*s/d)^2 t1<-qt(alpha/2,n0,lower.tail=FALSE) n1<-(t1*s/d)^2 while(abs(n1-n0)>0.5){ n0<-(qt(alpha/2,n1,lower.tail=FALSE)*s/d)^2 n1<-(qt(alpha/2,n0,lower.tail=FALSE)*s/d)^2 } n1 } #例5.6.2 size.norm2(10,0.01,2,100) #14 #估计比例p时样本容量的确定,定义函数size.bin() size.bin<-function(d,p,conf.level=0.95){ alpha=1-conf.level ((qnorm(1-alpha/2))/d)^2*p*(1-p) } #例5.6.3 size.bin(0.03,0.9,0.95) 第三部分、课后习题: 5.1 x=seq(3,21,by=2) y=c(21,16,15,26,22,14,21,22,18,25) u=rep(x,y) u1=mean(u) s1=sd(u) a=u1-sqrt(3)*s1 a b=u1+sqrt(3)*s1 b 因此,α,β的矩估计值为: 2.217379、22.40262 5.2 x=seq(0,6,by=1) y=c(17,20,10,2,1,0,0) x1=x*y mu=mean(x1) mu 结论: 平均每升水中大肠杆菌个数为7.142857时,才能使上述情况的概率达到最大. 5.3 (1) x=c(482,493,457,471,510,446,435,418,394,469) t.test(x) 结论: 置信水平为0.95的置信区间为[432.3069,482.6931] (2) chisq.var.test<-function(x,var,alpha,alternative="two.sided"){ options(digits=4) result<-list() n<-length(x) v<-var(x) result$var<-v; chi2<-(n-1)*v/var result$chi2<-chi2; p<-pchisq(chi2,n-1) result$p.value<-p if(alternative=="less") result$p.value<-pchaisq(chi2,n-1,loer.tail=F) elseif(alternative=="two.sider") result$p.value<-2*min(pchaisq(chi2,n-1), pchaisq(chi2,n-1,lower,tail=F)) result$conf.int<-c( (n-1)*v/qchisq(a
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 语言 实验