实验五 常见分布的相关计算、随机抽样与模拟
【实验类型】验证性 【实验学时】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 (i in 1: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 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=\"Normal Distribution\") curve(dnorm(x), add = TRUE, col=\"red\绘制正态密度曲线 expr<-expression(paste(mu==0, \ legend(1, 0.35, legend = expr, col=\"red\添加给定内容的图例 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) # 有放回抽样 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 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 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\main=\"Random Walk in Two Dimensions\") x<-y<-50 # 两变量同时赋值 points(50, 50, pch=16, col=\"red\绘制初始点 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 <- 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 } else if (alternative == \"greater\"|alternative ==\"less\" ){ result$p.value<-p } else return(\"your input is wrong\") 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 } else if (alternative==\"two.sided\") { if(p>.5) p<-1-p p<-2*p result$p.value<-p } else return(\"your input is wrong\") 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) else if(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(alpha/2,df=n-1,lower.tail=FALSE), (n-1)*v/qchisq(alpha/2,df=n-1,lower.tail=TRUE)) result } x<-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.0] 5.4 (1) X<-c(25,28,23,26,29,22) Y<-c(28,23,30,35,21,27) var.test(X,Y) 结论：由于P值（=0.2）>0.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.570] 5.5 two.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=qnorm(1-alpha/2)*(sigma1/m+sigma2/n)^(1/2) xbar +c(-zstar, +zstar) } x<-c(628,583,510,554,612,523,530,615,573,603,334,564) y<-c(535,433,398,470,567,480,498,560,503,426,338,547) sigma1=2140 sigma2=3250 two.sample.ci(x,y,conf.level=0.95,sigma1,sigma2) two.sample.ci(x,y,conf.level=0.90,sigma1,sigma2) 结论：置信水平为0.95的置信上限为：31.29 置信水平为0.90的置信下限为：107.69 5.6 x<-c(15.2,14.5,15.5,14.8,15.1,15.6,14.7) y<-c(15.2,15.0,14.8,15.2,15.0,14.9,15.1,14.8,15.3) var.test(x,y) 由于P值（=0.04）<0.05，因此可以认为机床乙生产的滚珠的方差比机床甲生产的滚珠直径的方差小。 5.7 binom.test(224,400,conf.level=0.99) 因此可以以99%的置信水平认为顾客中喜欢A的人数比例p落在（0.4944，0.6241）中，其点估计为0.56 5.8 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 } size.norm2(10,0.05,2,100) 因此在0.95的置信度下至少要抽取99个产品。 5.9 size.bin<-function(d, p, conf.level=0.95) { alpha = 1 - conf.level ((qnorm(1-alpha/2))/d)^2*p*(1-p) } size.bin(0.01, 0.05, 0.90) 因此要抽取1285个样本验收可满足上述要求。 因篇幅问题不能全部显示，请点此查看更多更全内容