R语言总和性试验.docx

中北大学理学院实验报告实验课程名称：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 D=A*B3、 E=A%*%B4、F G H=F%*%G3、函数solve()有两个作用；solve（A,b）可用于求解线性方程组Ax=b,solve(A)可用于求解矩阵A的逆，用两种方法编程求解方程组Ax=b的解。A A3,3=10 b=matrix(1:1,3,1) solve(A,b) ,11, -1.000000e+002, 1.000000e+003, 3.806634e-164、用三种方法求解它们的內积与外积。1、 x=c(1,2,3,4,5) y=c(2,4,6,8,10) a=t(x)%*%ye=x%*%t(y)2、 b=crossprod(x,y) f=outer(x,y)3、 c=x%*%y d=x%o%y5、编写一个用二分法求解非线性方程的根的函数，并求方程x3-x-1=0在区间1,2内的根，精度要求e=10-5.程序：fzero0)list(fail=finding root is fail!)elserepeatif(abs(b-a)eps)breakx-(a+b)/2if(f(a)*f(x)0) b-x else a f fzero(f,1,2,1e-5)$root1 1.324718$fun1 -1.405875e-05第三章：1、从1到100个自然数中随机不放回的抽取5个数，并求它们的和。 sum(sample(1:100,5)1 2052、从一副扑克牌（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-053、从正态分布N（100,100）中随机产生1000个随机数，（1）rnorm(1000,mean=100,sd=10)结果随机执行/产生1000个随机数作出这1000个正态随机数的直方图；hist( rnorm(1000,mean=100,sd=10)结果随机执行/产生对应直方图从这1000个随机数中随机有放回的抽取500个作出直方图（2） A 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.26114、模拟随机游动：从标准正态分布中产生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-testdata: rnorm(100, mean = 0, sd = 1)t = 0.31167, df = 99, p-value = 0.7559alternative hypothesis: true mean is not equal to 095 percent confidence interval: -0.1638460 0.2249099sample 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,col=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 plot(mydata$heightmydata$weight,xlab=Weight,ylab=height) lines(lowess(mydata$height,mydata$weight),lwd=2)(2)、绘制不同性别下体重对身高的散点图； mydata coplot(mydata$heightmydata$weight|mydata$性别,xlab=weight,ylab= height)（3）、绘制不同年龄段体重对身高的散点图； mydata coplot(mydata$heightmydata$weight|mydata$年龄,xlab=weight,ylab= height)（4）、绘制不同性别下不同年龄段体重对身高的散点图； mydata coplot(mydata$heightmydata$weight|mydata$性别*mydata$年龄,xlab=wei ght,ylab=height)第五章：1、设总体X是用无线电测距仪测量距离的误差，它服从（A,B）上的均匀分布在200次测量中，误差Xi的次数有ni次：求A,B的矩法估计值； x a b z=b/a A B x lambda lambda%输出矩估计值113、已知某种木材的横纹抗压力服从N(mu,sigma2),现在对十个试件作横纹抗压力实验数据如下：（1） 、求mu的置信水平为0.95的置信区间 x t.test(x)%方差未知时求均值的置信区间，调用t.test函数One Sample t-testdata: xt = 41.08, df = 9, p-value = 1.495e-11alternative hypothesis: true mean is not equal to 095 percent confidence interval: 432.3069 482.6931sample estimates:mean of x 457.5 置信区间为（432.3069， 482.6931）4、某卷烟厂生产两种卷烟A和B，现在分别对两种香烟的尼古丁含量进行六次测试，结果如下： A B var.test(A,B) F test to compare two variancesdata: A and BF = 0.3, num df = 5, denom df = 5, p-value = 0.2alternative hypothesis: true ratio of variances is not equal to 195 percent confidence interval: 0.04187 2.13821sample 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.levelzstar=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.39355、 为比较俩个小麦品种的产量，选择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.levelzstar=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-2140sigma2-3250two.sample.ci(x,y,conf.level=0.95,sigma1,sigma2)1 31.29 114.37two.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.levelzstar=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-2140sigma2-3250two.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 variancesdata: x and yF = 5.2, num df = 6, denom df = 8, p-value = 0.04alternative hypothesis: true ratio of variances is not equal to 195 percent confidence interval: 1.121 29.208sample 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 correctiondata: 224 out of 400, null probability 0.5X-squared = 5.5, df = 1, p-value = 0.02alternative hypothesis: true p is not equal to 0.595 percent confidence interval: 0.5098 0.6091sample 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)2size.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-testdata: vt = -5, df = 8, p-value = 0.001alternative hypothesis: true mean is not equal to 95095 percent confidence interval: 915.3 937.3sample 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-testdata: x and yt = 5.3, df = 21, p-value = 3e-05alternative hypothesis: true difference in means is not equal to 095 percent confidence interval: 0.7955 1.8212sample 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 X A miscellany aov.mis summary(aov.mis)#调用summary函数提取分析结果。以上结果中：Df表示自由度；sum sq表示平方和；mean sq表示均方和；F value表示F检验统计量的值，即为：F比；Pr(F)表示检验的p值；A就是因素A；Residuals为残差。可以看出P=0.0162 plot(miscellany$Xmiscellany$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)#按照缺省的“bonferroni”对p值进行调整。同时置信区间发：Tukey法TukeyHSD(aov(XA,miscellany)共有10个两两比较结果，在0.05的显著性水平下A5-A1，A5-A2，A5-A3，A5-A4的差异是显著的。其他结果均不显著。方差齐性检验：代码段： bartlett.test(X,A,data=miscellany)#方差齐性检验之“bartlett”检验。由p值0.13090.05,故：接受原假设，认为各处理组数据等方差。无交互作用双因子方差分析：例题8.2.1原来检验果汁中铅含量有三种方法，A1，A2，A3，现在研究出另一种检验方法A4,能否用A4代替前三种方法，要进行实验考察。实验数据如下：代码段：juice juice.aov summary(juice.aov)#提取结果。p值说明了因素B对含铅量有显著影响，而因素A对含铅量影响不显著。 bartlett.test(XA,data=juice)由于p=0.96590.05，所以接受原假设认为各处理组数据等方差。 bartlett.test(XB,data=juice)由于p=0.003766 rats op plot(TimeToxicant+Cure,data=rats)Hit to see next plot: with(rats,interaction.plot(Toxicant,Cure,Time,trace.label=Cure)#绘制交互效应图。 with(rats,interaction.plot(Cure,Toxicant,Time,trace.label=Toxicant)#绘制出交互效应图。无明显的相交情况，说明两个因素无交互作用。 rats.aov summary(rats.aov)#方差分析根据p值可知，因素toxicant和Cure对time的影响是高度显著的，而且交互作用对time的影响却是不显著的。第九章 回归分析(1）、 由专业知识可知，合金的强度Y与合金中碳含量X有关，现从生产中收集了一批数据(xi, yi) (i=l,2,n) ，见下表，试分析合金的强度Y与合金中碳含量X之间的关系。(1)完成一元线性回归的计算；(2)计算自变量x在区间0.10, 0.23内的回归方程的预测估计值、预测区间和置信区间(取=0.05)，并将数据点、预测估计曲线、预测区间曲线和置信区间曲线画在同一张图上。1) 回归估计x-c(0.10, 0.11, 0.12, 0.13, 0.14, 0.15, 0.16,0.17, 0.18, 0.20, 0.21, 0.23)y-c(42.0, 43.5, 45.0, 45.5, 45.0, 47.5, 49.0,53.0, 50.0, 55.0, 55.0, 60.0)lm.sol|t|) (Intercept) 28.493 1.580 18.04 5.88e-09 *x 130.835 9.683 13.51 9.50e-08 *Signif. codes: 0 * 0.001 * 0.01 * 0.05 . 0.1 1Residual standard error: 1.319 on 10 degrees of freedomMultiple R-squared: 0.9481,Adjusted R-squared: 0.9429 F-statistic: 182.6 on 1 and 10 DF, p-value: 9.505e-08计算结果表明，系数和方程均通过检验，回归方程为：Y=28.493+130.835X(2) 计算预测值并绘图new - data.frame(x = seq(0.10, 0.24, by=0.01)pp-predict(lm.sol, new, interval=prediction);pc-predict(lm.sol, new, interval=confidence);par(mai=c(0.8, 0.8, 0.2, 0.2)matplot(new$x, c