R语言实验五.doc

_实验五 常见分布的相关计算、随机抽样与模拟【实验类型】验证性【实验学时】2 学时【实验目的】1、掌握常见分布的分布函数、密度函数（或分布列）及分位数的计算方法；2、掌握样本统计量的计算方法及所表达的意义；3、了解随机模拟的基本思想及其应用。【实验内容】1、组合数与组合方案的生成、概率的计算，2、常见分布的分布函数、密度函数（或分布列）以及分位数的计算；3、随机数的生成与随机模拟（蒙特卡洛仿真）。【实验方法或步骤】第一部分、课件例题：1.#从15个数中，随机取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 (i in 1:3) pi-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 ：计算PX160，其中XU150,200。p-punif(160, min=150, max=200); p# 例5.6 ：已知XE(1/241)，计算P50X100 。p-pexp(100, rate=1/241)-pexp(50, rate=1/241); p# 例5.7 ：已知XB(80,0.01)，计算PX4 。p - 1-pbinom(3, size=80, prob=0.01); p# 例5.8 ：已知XB(180,0.01)，且PXk0.95，求k 。n - 180; p - 0.01; k - qbinom(0.95, n, p); k# 例5.9 ：已知XB(1000,0.01)，求PX2 。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(m23) else n - length(x); S - sd(x) Cs-n/(n-1)/(n-2)*sum(x-mu)3)/S3 Csx-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=Normal Distribution)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(i in 1:B) sam - sample(x, replace=TRUE) # 有放回抽样 mediansi - 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 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(yfun(x) # 落在曲边梯形内的随机点个数 k/n*M*(b-a) #k/n 表示梯形面积与矩形面积的比值#例5.13fun-function(x) exp(-x2) # 被积函数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(x2+y2=1)*(1-x2-y2) # 定义被积函数，定义域内函数值不变，定义域外函数值为 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=Random Walk in Two Dimensions)x-y-50 # 两变量同时赋值points(50, 50, pch=16, col=red, cex=1.5) # 绘制初始点for (i in 1: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 100 | x100 | 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)(P517)*(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-mean result$z-z result$p.value-p if(alternative=two.sided) p-2*p result$p.value-p else if (alternative = greater|alternative =less ) result$p.value-p else return(your input is wrong) result$- 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#求置信区间的程序:-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.1x-c(175 ,176 ,173,175,174,173,173,176,173,179 )result-z.test(x,10,1.5,0.05)result$x-c(175 ,176 ,173,175,174,173,173,176,173,179 )(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)$ #提取出置信区间的部分#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.5) p-1-p p-2*p result$p.value-p else return(your input is wrong) result$-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.1x-c(628,583,510,554,612,523,530,615)y-c(535,433,398,470,567,480,498,560,503,426)sigma1-2140sigma2-3250two.sample.ci(x,y,conf.level=0.95, sigma1,sigma2)#8#例5.3.2x-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.3x-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.1prop.test(38,200,correct=TRUE) #函数prop.test( )对p进行估计与检验binom.test(38,200) #函数binom.test( )可以求其置信区间#11#例5.5.1like-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)2size.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) n10.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.2size.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.3size.bin(0.03, 0.9, 0.95)第三部分、课后习题：5.1x=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)*s1ab=u1+sqrt(3)*s1b因此，,的矩估计值为：2.217379、22.402625.2x=seq(0,6,by=1)y=c(17,20,10,2,1,0,0)x1=x*ymu=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) else if(alternative=two.sider) result$p.value-2*min(pchaisq(chi2,n-1), pchaisq(chi2,n-1,lower,tail=F) result$-c( (n-1)*v/qchisq(alpha/2,df=n-1,lower.tail=FALSE), (n-1)*v/qchisq(alpha/2,df=n-1,lower.tail=TRUE) resultx-c(482,493,457,471,510,446,435,418,394,469)chisq.var.test(x,var(x),0.10,alternative=two.side)结论:的置信水平为0.90的置信区间为659.8,3357.05.4(1)X-c(25,28,23,26,29,22)Y0.05，因此认为两种卷烟中尼古丁含量的方差相等。#(2)X-c(25,28,23,26,29,22)Y-c(28,23,30,35,21,27)t.test(X,Y,var.equal = FALSE)结论：两种香烟的尼古丁平均含量差的95%置信区间为：-7.237,3.5705.5two.sample.ci=function(x,y,conf.level=0.95,sigma1,sigma2)options(digists=4) m=length(x);n=length(y) xbar=mean(x)-mean(y) alpha=1-conf.level zstar