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 (i in 1: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{50<X≤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);x sort(x) #默认升序，不处理缺失数据sort(x,decreasing=TRUE) #降序sort(x,st=TRUE) #处理缺失数据sort(x,decreasing =TRUE,st=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", 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) # 有放回抽样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<fun(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(z<fun(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="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 <- y + yi # 下一步行走位置if (x>100 | x<0 | y>100 | y<0) 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}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.1x<-c(175 ,176 ,173,175,174,173,173,176,173,179 )result<-z.test(x,10,1.5,0.05)result$conf.intx<-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<-vchi2<-(n-1)*v/varresult$chi2<-chi2p<-pchisq(chi2,n-1)if(alternative == "less"|alternative=="greater"){result$p.value<-p} else if (alternative=="two.sided") {if(p>.5)p<-1-pp<-2*presult$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.levelzstar= 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/n2xbar=xbar1-xbar2alpha = 1 - conf.levelzstar=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)^2t1<-qt(alpha/2,n0,lower.tail=FALSE)n1<-(t1*s/d)^2while(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.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/varresult$chi2<-chi2;p<-pchisq(chi2,n-1)result$p.value<-pif(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，因此认为两种卷烟中尼古丁含量的方差相等。