R语言总和性试验
中北大学理学院
实验报告
实验课程名称：R语言与统计分析实验类别：验证型
专业：应用统计学
班级： 13080441
学号： 1308044142
姓名：吴庚雷
中北大学理学院
R语言与统计分析综合实验
【实验类型】验证性
【实验目的】
（1）掌握利用R语言实现数据处理并进行严格的统计分析；
（2）学会运用R语言进行程序的编写；
（3）熟练掌握R语言绘图功能；
（4）掌握R语言统计分析中的“参数估计”，“假设检验”，“方差分析”，“回归分析”，等基本分析函数。
【实验要求】
（1）实验过程要求用R软件完成；
（2）实验结果逐个导入Word文档，并按问题作出解释；
（3）实验报告按照既定格式书写。
【实验仪器与软件平台】
计算机 R软件
【实验前的预备知识】
1、实验室电脑要求安装有R软件；
2、上实验课程的学生要对涉及到的统计概念有所了解；
3、要求学生事先查阅并熟悉R的相关命令。
【实验内容】
第二章：
1、
用rep()构造一个向量x，它由3个3,4个2,5个1构成。
x<-rep(c(3,2,1),c(3,4,5))
2、
由1.2...16构成两个方阵，其中矩阵A按列输入，矩阵B按行输入，并计算以下：
A<-matrix(1:16,4,4)
B<-matrix(1:16,4,4,byrow=TRUE)
1、
C=A+B
2、
>D=A*B
3、
>E=A%*%B
4、
F<-A[-3,]
[,1] [,2] [,3] [,4]
[1,] 1 5 9 13 [2,] 2 6 10 14 [3,] 4 8 12 16
>G<-B[,-3]
[,1] [,2] [,3]
[1,] 1 2 4
[2,] 5 6 8
[3,] 9 10 12
[4,] 13 14 16
>H=F%*%G
3、
函数solve()有两个作用；solve（A,b）可用于求解线性方程组Ax=b,solve(A)可用于求解矩阵A的逆，用两种方法编程求解方程组Ax=b的解。
A<-matrix(1:9,3,3,byrow=TRUE)
>A[3,3]=10
>b=matrix(1:1,3,1)
>solve(A,b)
[,1]
[1,] -1.000000e+00
[2,] 1.000000e+00
[3,] 3.806634e-16
4、
用三种方法求解它们的內积与外积。
1、
>x=c(1,2,3,4,5)
>y=c(2,4,6,8,10)
>a=t(x)%*%y
e=x%*%t(y)
2、
>b=crossprod(x,y)
>>f=outer(x,y)
3、
>c=x%*%y
>d=x%o%y
5、
编写一个用二分法求解非线性方程的根的函数，并求方程x^3-x-1=0在区间[1,2]内的根，精度要求e=10^-5.
程序：
fzero<-function(f,a,b,eps=1e-5)
{
if(f(a)*f(b)>0)list(fail="finding root is fail!")
else{repeat{if(abs(b-a) break x<-(a+b)/2 if(f(a)*f(x)<0) b<-x else a<-x} list(root=(a+b)/2,fun=f(x))}} >f<-function(x) x^3-x-1 >fzero(f,1,2,1e-5) $root [1] 1.324718 $fun [1] -1.405875e-05 第三章： 1、 从1到100个自然数中随机不放回的抽取5个数，并求它们的和。 >sum(sample(1:100,5)) [1] 205 2、 从一副扑克牌（52张）中随机抽取5张求以下概率 抽到的是10，J,Q,K,A; （1）>4/choose(52,5) [1] 1.539077e-06 抽到的是同花顺 （2）9*choose(4,1)/choose(52,5) [1] 1.385169e-05 3、 从正态分布N（100,100）中随机产生1000个随机数， （1）rnorm(1000,mean=100,sd=10) 结果随机执行//产生1000个随机数 作出这1000个正态随机数的直方图； hist( rnorm(1000,mean=100,sd=10)) 结果随机执行//产生对应直方图 从这1000个随机数中随机有放回的抽取500个作出直方图（2）>A<- rnorm(1000,mean=100,sd=10) >sample(A,500,replace=TRUE) 或者sample(rnorm(1000,mean=100,sd=10),500) 结果随机执行//产生500个对应随机数 >hist(sample(A,500,replace=TRUE)) 结果随机执行//产生对应直方图 比较它们的样本均值和样本方差 （3）>mean(rnorm(1000,mean=100,sd=10)) [1] 100.0266 >var(rnorm(1000,mean=100,sd=10)) [1] 98.5499 >mean( sample(rnorm(1000,mean=100,sd=10),500)) [1] 99.9142 >var( sample(rnorm(1000,mean=100,sd=10),500)) [1] 106.2611 4、 模拟随机游动：从标准正态分布中产生1000个随机数，并用函数cumsum（）作出累积和，最后使用命令plot()作出随机游动示意图 >rnorm(1000,mean=0,sd=1) //产生1000个标准正态分布随机数 cumsum(rnorm(1000,mean=0,sd=1)) //累积和 plot(cumsum(rnorm(1000,mean=0,sd=1))) //随机执行产生随机游动图 5、 从标准正态分布中随机产生100个随机数，由此数据求总体均值的95%的置信区间 t.test(rnorm(100,mean=0,sd=1)) //t.test(c(数据))专用于求出95%的置信区间 结果如下： One Sample t-test data: rnorm(100, mean = 0, sd = 1) t = 0.31167, df = 99, p-value = 0.7559 alternative hypothesis: true mean is not equal to 0 95 percent confidence interval: -0.1638460 0.2249099 sample estimates: mean of x 0.03053192 第四章： 1、模拟得到1000个参数为0.3的伯努利分布随机数，并用图示表示出来 > plot(rbinom(1000,1,0.3)) 2、 用命令rnorm()命令产生1000个均值为10，方差为4的正态分布随机数，用直方图呈现数据分布并添加核密度曲线； >x=rnorm(1000,mean=10,sd=2)%生成正态分布随机数 > hist(x,xlim=c(min(x),max(x)),probability=T,nclass=max(x)-min(x)+1, col='lightblue',main='Normal distribution,mean=10,sd=2')%绘制直方图 >lines(density(x,bw=1),col='red',lwd=3)%添加核密度曲线 7、 假定某校100名女生的血清总蛋白含量服从均值为75，标准差为3的分布，并假定数据由下面命令产生： >options(digits=4) >x=rnorm(100,75,3)%生成随机数 计算样本均值，标准差，以及五数概括； > summary(x)%计算均值和五数 Min. 1st Qu. Median Mean 3rd Qu. Max. 65.7 72.4 74.2 74.7 76.7 82.3 > sd(x)%标准差 [1] 3.338 绘制出直方图，核密度估计曲线，和QQ图； >hist(x,xlim=c(min(x),max(x)),probability=T,nclass=max(x)-min(x)+1,c ol='lightgreen',main='Normal distrubution,mean=75,sd=3')%绘制直方图 > lines(density(x,bw=1),col='red',lwd=3)%添加核密度曲线 >qqnorm(x,main="Normality Check via QQ Plot") > qqline(x,col='red')%绘制QQ图 > qqline(x,col='red')%绘制QQ图 根据数据绘制出茎叶图和框须图； > stem(x)%绘制茎叶图 The decimal point is at the | 68 | 344 70 | 16889937 72 | 0112233667902334556777789 74 | 0012244556688990023588 76 | 001223667779011234459 78 | 000255790158 80 | 2378906 82 | 20 > boxplot(x)%绘制框须图 8、 某校测得20名学生的4项指标：性别，年龄，身高，体重，具体数据如课本表4.1所示： (1)、绘制体重对身高的散点图； >mydata<-read.delim("clipboard") >plot(mydata$height~mydata$weight,xlab="Weight",ylab="height") > lines(lowess(mydata$height,mydata$weight),lwd=2) (2)、绘制不同性别下体重对身高的散点图； >mydata<-read.delim("clipboard") > coplot(mydata$height~mydata$weight|mydata$性别,xlab="weight",ylab= " height") （3）、绘制不同年龄段体重对身高的散点图； >mydata<-read.delim("clipboard") > coplot(mydata$height~mydata$weight|mydata$年龄,xlab="weight",ylab= "height") （4）、绘制不同性别下不同年龄段体重对身高的散点图； >mydata<-read.delim("clipboard") > coplot(mydata$height~mydata$weight|mydata$性别*mydata$年龄,xlab="w eight",ylab="height") 第五章： 1、 设总体X是用无线电测距仪测量距离的误差，它服从（A,B）上的均匀分布在200次测量中，误差Xi的次数有ni次：求A,B的矩法估计值； >x<-rep(c(3,5,7,9,11,13,15,17,19,21),c(21,16,15,26,22,14,21,22,18,2 5))%读入数据 > a<-2*mean(x) > b<-12*var(x) > z=b/a > A<-(a-z)/2%A的矩估计值 A=4.035345 > B<-(a-A)%B的矩估计值 B=20.58466 2、为检验某自来水消毒设备的效果，从消毒后的水中随机抽取50L，化验每L 水中大肠杆菌的个数（假设其服从泊松分布），化验结果如表中所示，试问平均每升水中大肠杆菌的个数为多少时，才能使上述情况概率达到最大； > x<-rep(c(0,1,2,3,4,5,6),c(17,20,10,2,1,0,0))%读取数据 > lambda<-mean(x)%定义泊松分布的矩估计 > lambda%输出矩估计值 [1]1 3、 已知某种木材的横纹抗压力服从N(mu,sigma^2),现在对十个试件作横纹抗压力实验数据如下： （1）、求mu的置信水平为0.95的置信区间 > x<-c(482,493,457,471,510,446,435,418,394,469)%读取数据 > t.test(x)%方差未知时求均值的置信区间，调用t.test函数 One Sample t-test data: x t = 41.08, df = 9, p-value = 1.495e-11 alternative hypothesis: true mean is not equal to 0 95 percent confidence interval: 432.3069 482.6931 sample estimates: mean of x 457.5 置信区间为（432.3069， 482.6931） 4、 某卷烟厂生产两种卷烟A和B，现在分别对两种香烟的尼古丁含量进行六次测试，结果如下： > A<-c(25,28,23,26,29,22) > B<-c(28,23,30,35,21,27) （1）、试问两种卷烟中的尼古丁含量的方差是否相等？ > var.test(A,B) F test to compare two variances data: A and B F = 0.3, num df = 5, denom df = 5, p-value = 0.2 alternative hypothesis: true ratio of variances is not equal to 1 95 percent confidence interval: 0.04187 2.13821 sample estimates: ratio of variances 0.2992 (2)、试求两种卷烟的尼古丁平均含量差的95%的置信区间； 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) } x<-c(25,28,23,26,29,22) y<-c(28,23,30,35,21,27) sigma1=sqrt(var(x)) sigma2=sqrt(var(y)) two.sample.ci(x,y,conf.level=0.95,sigma1,sigma2) [1] -4.0602 0.3935 5、 为比较俩个小麦品种的产量，选择22块条件相似的试验田，采用相同的耕作方法做试验，结果播种甲品种的12块试验田的单位面积产量和播种乙的12块试验田的单位面积产量如表所示，假定每个品种的单位面积产量均服从正太分布，甲品种产量方差为2140，乙品种产量方差为3250.试求这俩个品种平均面积产量差的置信水平为0.95的置信上限和置信水平为0.90的置信下限； 程序如下：编写函数求产量差的直线区间 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) } 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) [1] 31.29 114.37 two.sample.ci<-function(x,y,conf.level=0.90,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) } 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.90,sigma1,sigma2)%调用函数进行双样本 已知方差估计均值 [1] 37.97 107.69 因此置信水平为0.95的置信上限为：114.37；置信水平为0.90的置信下限为：37.97 (6)、 有两台机床生产同一型号的滚珠根据以往经验得知，滚珠的直径都服从正态分布，现在从这两台机床生产的滚珠中随机抽取7个和9个，测得其直径的数据如下：试问机床乙生产的滚珠的方差是否比甲的方差小？ 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)%调用函数进行双样本方差比估计 F test to compare two variances data: x and y F = 5.2, num df = 6, denom df = 8, p-value = 0.04 alternative hypothesis: true ratio of variances is not equal to 1 95 percent confidence interval: 1.121 29.208 sample estimates: ratio of variances 5.216 (7)、 某公司对本公司生产的两种自行车型号A,B的销售情况进行调查，随机选取400人询问他们对A,B的选择，其中有224人喜欢A，试求顾客中喜欢A的人数比例p的置信水平为99%的区间估计； prop.test(224,400,correct=TRUE)%调用函数进行单总体比率p的区间估计 1-sample proportions test with continuity correction data: 224 out of 400, null probability 0.5 X-squared = 5.5, df = 1, p-value = 0.02 alternative hypothesis: true p is not equal to 0.5 95 percent confidence interval: 0.5098 0.6091 sample estimates: p 0.56 所以我们以0.95的置信水平认为比例p落在（0.5098，0.6091）中，点估计为0.56 （8）、 某公司生产一批新产品，产品总体服从正态分布，现要估计这批产品的平均重量，最大允许误差为1，样本标准差s=10,试问在95%的置信水平下至少要抽取多少个产品？ %定义函数 size.norm1<-function(d,var,conf.level){ alpha=1-conf.level ((qnorm(1-alpha/2)*var^(1/2))/d)^2 } size.norm1(1,100,conf.level=0.95)%调用函数 [1] 384.1 所以应该抽取384个产品 (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)%调用函数 [1] 1285 需要抽取1285个样本 第六章： (1)、 v<-c(914,920,910,934,953,940,912,924,930) t.test(v,mu=950) One Sample t-test data: v t = -5, df = 8, p-value = 0.001 alternative hypothesis: true mean is not equal to 950 95 percent confidence interval: 915.3 937.3 sample estimates: mean of x 926.3 由结果可知：p值为0.001显著小于0.01，所以应该拒绝原假设，即：认为这批枪弹的初速显著降低。 (3)、 下面给出了俩种的计算器充电以后所能使用的时间的观测值： 设俩样本独立且数据所属的俩个总体的密度函数至多差一个平移量，试问能否认为型号A的计算器平均使用时间比B长（0.01）？ x<-c(5.5,5.6,6.3,4.6,5.3,5.0,6.2,5.8,5.1,5.2,5.9) y<-c(3.8,4.3,4.2,4.9,4.5,5.2,4.8,4.5,3.9,3.7,3.6,2.9) t.test(x,y,var.equal=TRUE) Two Sample t-test data: x and y t = 5.3, df = 21, p-value = 3e-05 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: 0.7955 1.8212 sample estimates: mean of x mean of y 5.500 4.192 由结果可知：因为p值显著小于0.01，因此拒绝原假设，即： (4)、 测得俩批电子器件的样本电阻 设这俩批器材的电阻值分别服从正态分布N和M，样本独立。 (1）试样本俩个总体的方差是否相等（0.01）？ (2）试样本俩个总体的均值是否相等（0.05）？ x<-c(0.140,0.138,0.143,0.142,0.144,0.137) y<-c(0.135,0.140,0.142,0.136,0.138,0.130) t.test(x,y,var.equal=TRUE) Two Sample t-test data: x and y t = 1.9, df = 10, p-value = 0.09 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: -0.0007721 0.0084388 sample estimates: mean of x mean of y 0.1407 0.1368 由结果可知：p值为0.09，大于显著性水平0.05，说明应该接受原假设，即：均值显著相等 var.test(x,y) F test to compare two variances data: x and y F = 0.44, num df = 5, denom df = 5, p-value = 0.4 alternative hypothesis: true ratio of variances is not equal to 1 95 percent confidence interval: 0.06196 3.16425 sample estimates: ratio of variances 0.4428 由结果可知：p值为0.4，大于显著性水平0.01，说明应该接受原假设，即：方差显著相等 (5)、 有人称某地成年人中大学毕业生比例不低于30%，为检验之，随机调查该地15名成年人，发现有3名大学毕业生，有0.95的置信度，问该人的看法是否成立？binom.test(c(12,3),p=0.3)%单样本比率p的近似检验 Exact binomial test data: c(12, 3) number of successes = 12, number of trials = 15, p-value = 9e-05 alternative hypothesis: true probability of success is not equal to 0.3 95 percent confidence interval: 0.5191 0.9567 sample estimates: probability of success 0.8 由结果可知：p值显著小于显著性水平0.05，所以应该拒绝原假设，即：此人看法不成立 第八章 单因素方差分析： 例题8.1.1 以淀粉为原料生产葡萄糖的过程中，残留许多糖蜜，可作为生产酱色的原料，在生产酱色的过程之前应该尽可能的除去杂质，以保证酱色质量。为此对于除杂方法进行选择，在实验中用五种不同除杂方式，每种作4次试验，数据如下： 代码段： > X<-c(25.6,22.2,28.0,29.8,24.4,30.0,29.0,27.5,25.0,27.7,23.0,32.2,2 8.8,28.0,31.5,25.9,20.6,21.2,22.0,21.2) > A<-factor(rep(1:5,each=4))#定义各因子变量的数据读入分段。 > miscellany<-data.frame(X,A)#将定义完成各因子变量的分段混合后赋给“miscellany”混合变量。 > aov.mis<-aov(X~A,data=miscellany)#进行方差分析。 > summary(aov.mis)#调用summary函数提取分析结果。 以上结果中：Df表示自由度；sum sq表示平方和；mean sq表示均方和；F va lue表示F检验统计量的值，即为：F比；Pr(>F)表示检验的p值；A就是因素A；Residuals为残差。 可以看出P=0.0162<0.05，说明有理由拒绝原假设，即，有显著性差异。 > plot(miscellany$X~miscellany$A)#绘制以上数据的箱线图。 均值的多重比较： 代码段： P.adjust.methods#调用处调整p值得方法列表。 Pairwise.t.test(X,A,p.adjust.method=”none”)//应当注意“none”的符号问题。 >pairwise.t.test(X,A,p.adjust.method='none') #得到多重比较p值，“none”指的是不作任何调整,默认值按照holm调整。 Pairwise.t.test(X,A,p.adjust.method=’holm’)//应当注意首字母不能大写> pairwise.t.test(X,A,p.adjust.method='holm') #按照缺省的“holm”对p值进行调整。 >pairwise.t.test(X,A,p.adjust.method='bonferroni')#按照缺省的“bonfer roni”对p值进行调整。 同时置信区间发：Tukey法 >TukeyHSD(aov(X~A,miscellany)) 共有10个两两比较结果，在0.05的显著性水平下A5-A1，A5-A2，A5-A3，A5-A4的差异是显著的。其他结果均不显著。 方差齐性检验： 代码段： >bartlett.test(X,A,data=miscellany)#方差齐性检验之“bartlett”检验。 由p值0.1309>0.05,故：接受原假设，认为各处理组数据等方差。 无交互作用双因子方差分析： 例题8.2.1 原来检验果汁中铅含量有三种方法，A1，A2，A3，现在研究出另一种检验方法A4,能否用A4代替前三种方法，要进行实验考察。实验数据如下： 代码段： >juice<-data.frame(X=c(0.05,0.46,0.12,0.16,0.84,1.30,0.08,0.38,0.4,0. 10,0.92,1.57,0.11,0.43,0.05,0.10,0.94,1.10,0.11,0.44,0.08,0.03,0.93, 1.15),A=gl(4,6),B=gl(6,1,24))#数据输入与整理。 > juice.aov<-aov(X~A+B,juice)#无交互双因素方差分析。 > summary(juice.aov)#提取结果。 p值说明了因素B对含铅量有显著影响，而因素A对含铅量影响不显著。 > bartlett.test(X~A,data=juice) 由于p=0.9659>0.05，所以接受原假设认为各处理组数据等方差。 > bartlett.test(X~B,data=juice) 由于p=0.003766<0.05所以拒绝原假设，认为各处理组数居方差不等。 有交互作用的方差分析： 代码段： > rats<-data.frame(Time=c(0.31,0.45,0.46,0.43,0.82,1.10,0.88,0.72,0. 43,0.45,0.63,0.76,0.45,0.71,0.66,0.62,0.38,0.29,0.40,0.23,0.92,0.61, 0.49,1.24,0.44,0.35,0.31,0.40,0.56,1.02,0.71,0.38,0.22,0.21,0.18,0.2 3,0.30,0.37,0.38,0.29,0.23,0.25,0.24,0.22,0.30,0.36,0.31,0.33),Toxic ant=gl(3,16,48,labels=c("M","N","K")),Cure=gl(4,4,48,labels=c("Q","W ","R","T")))#定义数据。 > op<-par(mfrow=c(1,2)) > plot(Time~Toxicant+Cure,data=rats) Hit