统计建模与R软件实验报告

1、开课学院、实验室： 数学与统计学院 实验时间 ： 2013 年 3 月 日实验项目名 称多维数组和矩阵实验项目类型验证演示综合设计其他指导教师张应应成 绩一、实验目的1.了解R软件的基本功能以及基本操作2.掌握R软件的基本操作二、实验内容 生成一个5阶的Hilbert矩阵（1） 计算Hilbert矩阵H的行列式（2） 求H的逆矩阵 （3） 求H的特征值和特征向量。三、实验原理、方法（算法）、步骤 1.函数det(A)是求矩阵A的行列式的值 2.求矩阵A的逆，其命令形式为slove(A) 3.函数eigen（Sm）是求对称矩阵Sm的特征值与特征向量n<-5;x<-array(0,di

2、m=c(n,n)for (i in 1:n)for (j in 1:n)xi,j<-1/(i+j-1);xdet(x)solve(x)eigen(x)四、实验环境（所用软件、硬件等）及实验数据文件 软件： R 2.15.3 五、实验结果及实例分析> n<-5;x<-array(0,dim=c(n,n)> for (i in 1:n)+ for (j in 1:n)+ xi,j<-1/(i+j-1)+ + ;x ,1 ,2 ,3 ,4 ,51, 1.0000000 0.5000000 0.3333333 0.2500000 0.20000002, 0.5000

3、000 0.3333333 0.2500000 0.2000000 0.16666673, 0.3333333 0.2500000 0.2000000 0.1666667 0.14285714, 0.2500000 0.2000000 0.1666667 0.1428571 0.12500005, 0.2000000 0.1666667 0.1428571 0.1250000 0.1111111> det(x)1 3.749295e-12> solve(x) ,1 ,2 ,3 ,4 ,51, 25 -300 1050 -1400 6302, -300 4800 -18900 268

4、80 -126003, 1050 -18900 79380 -117600 567004, -1400 26880 -117600 179200 -882005, 630 -12600 56700 -88200 44100> eigen(x)$values1 1.567051e+00 2.085342e-01 1.140749e-02 3.058980e-04 3.287929e-06$vectors ,1 ,2 ,3 ,4 ,51, 0.7678547 0.6018715 -0.2142136 0.04716181 0.0061738632, 0.4457911 -0.2759134

5、0.7241021 -0.43266733 -0.1166927473, 0.3215783 -0.4248766 0.1204533 0.66735044 0.5061636584, 0.2534389 -0.4439030 -0.3095740 0.23302452 -0.7671911935, 0.2098226 -0.4290134 -0.5651934 -0.55759995 0.376245545分析：从实验结果来看。R软件在处理数据上相当准确，方便。教师签名年 月 日开课学院、实验室：数学与统计学院 实验时间 ： 2013 年 3 月 日实验项目名 称描述统计量实验项目类型验证演

6、示综合设计其他指导教师张应应成 绩一、实验目的1.掌握利用R软件描述统计量的方法2.学会利用R软件绘制数据的分布图3.掌握多元数据的数据特征与相关分析在R软件中的操作二、实验内容 用Pearson相关检验法检验习题3.7中的身高与体重是否相关3、 实验原理、方法（算法）、步骤Pearson相关性检验：利用统计量t服从自由度为n-2的t分布的性质，对数据X和Y的相关性进行检验。其中。四、实验环境（所用软件、硬件等）及实验数据文件 R 2.15.3 数据文件 3.7数据.txt，其内容如下：学号 姓名 性别 年龄 身高 体重01 Alice F 13 56.5 84.002 Becka F 13

7、65.3 98.003 Gail F 14 64.3 90.004 Karen F 12 56.3 77.005 Kathy F 12 59.8 84.506 Mary F 15 66.5 112.007 Sandy F 11 51.3 50.508 Sharon F 15 62.5 112.509 Tammy F 14 62.8 102.510 Alfred M 14 69.0 112.511 Duke M 14 63.5 102.512 Guido M 15 67.0 133.013 James M 12 57.3 83.014 Jeffrey M 13 62.5 84.015 John

8、M 12 59.0 99.516 Philip M 16 72.0 150.017 Robert M 12 64.8 128.018 Thomas M 11 57.5 85.019 William M 15 66.5 112.0五、实验结果及实例分析 student<-read.table("3.7数据.txt",header=T) attach(student)> cor.test(身高,体重) #Pearson相关性检验 Pearson's product-moment correlationdata: 身高 and 体重 t = 7.5549, d

9、f = 17, p-value = 7.887e-07alternative hypothesis: true correlation is not equal to 0 95 percent confidence interval: 0.7044314 0.9523101 sample estimates: cor 0.8777852 其p值7.887e-07<0.05，拒绝原假设，所以身高与体重相关教师签名年 月 日开课学院、实验室： 数学与统计学院 实验时间 ： 2013年 月 日实验项目名 称参数估计实验项目类型验证演示综合设计其他指导教师张应应成 绩一、实验目的1、学会利用R软

10、件完成一个和两个正态总体的区间估计。2、学会利用R软件完成非正态总体的区间估计。3、学会利用R软件进行单侧置信区间估计二、实验内容 正常人的脉搏平均每分钟72次，某一声测得10例四乙基铅中毒患者的脉搏数（次/min）如下：54 67 68 78 70 66 67 70 65 69已知人的脉搏次数服从正态分布，试计算这10名患者平均脉搏次数的点估计和95%的区间估计，并做单侧区间估计，试分析这10名患者的平均脉搏次数是否低于正常人的平均脉搏次数。三、实验原理、方法（算法）、步骤假设正态总体，为来自总体的一个样本，为置信度，为样本均值，为样本方差。分别讨论总体方差已知和未知情况下，均值的单侧置信区

11、间估计。若已知，得到，于是得到的置信度为的单侧置信区间分别为，因此，的置信度为的单侧置信下限、上限分别为：，若未知，得到，于是得到的置信度为的单侧置信区间分别为，因此，的置信度为的单侧置信下限、上限分别为，4、 实验环境（所用软件、硬件等）及实验数据文件 数据见实验内容，所用软件：R2.15.15、 实验结果及实例分析在R软件中运行代码：> x<-c(54,67,68,78,70,66,67,70,65,69)> t.test(x) #做单样本正态分布区间估计 One Sample t-testdata: x t = 35.947, df = 9, p-value = 4.9

12、38e-11alternative hypothesis: true mean is not equal to 0 95 percent confidence interval: 63.1585 71.6415 sample estimates:mean of x 67.4 #平均脉搏点估计为 67.4 ， 95%置信度的区间估计为 63.1585, 71.6415 。> t.test(x,alternative="less",mu=72) #做单样本正态分布单侧区间估计 One Sample t-testdata:

13、x t = -2.4534, df = 9, p-value = 0.01828alternative hypothesis: true mean is less than 72 95 percent confidence interval: -Inf 70.83705 sample estimates:mean of x 67.4 p-value = 0.01828<0.05，拒绝原假设，平均脉搏低于常人。教师签名年 月 日开课学院、实验室： 实验时间 ： 2013 年 月 日实验项目名 称假设检验实验项目类型验证演示综合设计其他指导教师张应应成 绩一、实验目的1、了解假设检验的基本概

14、念与步骤2、掌握几种重要的参数检验3、掌握若干重要的非参数检验二、实验内容 为研究某中心要对抗凝血酶活力的影响，随机安排新药组病人12例，对照组病人10例，分别测定其抗凝血酶活力，其结果如下：新药组：126 125 136 128 123 138 142 116 110 108 115 140对照组：162 172 177 170 175 152 157 159 160 162试分析新药组和对照粗病人的抗凝血酶活力有无差别（）3、 实验原理、方法（算法）、步骤 1.基本原理：解决一个具体的假设检验问题，一般要借助直观分析和理论分析思想。其基本原理是实际推断原理：“小概率事件在一次实验中几乎不可

15、能发生”，如果发生，就认为是不正常的，应该拒绝。  2.方法： （1）提出原假设Ho  （2）确定假设检验统计量Z，并在Ho成立的条件下，导出Z的分布  （3）确定拒绝域：由直观分析先确定拒绝的形式，然后由显著水平及Z的分布P确定拒绝域的临界值，进而确定拒绝域C  （4）根据具体的一次样本值做出推断 3.具体算法： 假设是来自总体的样本，是来自总体的样本，且两样本独立，其检验问题有：双边检验： ，单边检验I ，单边检验II ，方差未知，和分别是和的样本方差。由统计知识可知，当为真时，其中因此，当T满足（成为拒绝域）：双边检验： 单边检验I 单边

16、检验II 则认为不成立，此方法也称为t检验法。4、 实验环境（所用软件、硬件等）及实验数据文件 见实验内容 软件：R2.15.3R软件。五、实验结果及实例分析> a <- c(126,125,136,128,123,138,142,116,110,108,115,140)> b <- c(162,172,177,170,175,152,157,159,160,162)#正态性检验：> ks.test(a,"pnorm",mean(a),sd(a) One-sample Kolmogorov-Smirnov testdata: a D = 0.1

17、464, p-value = 0.9266alternative hypothesis: two-sided > ks.test(b,"pnorm",mean(b),sd(b) One-sample Kolmogorov-Smirnov testdata: b D = 0.2222, p-value = 0.707alternative hypothesis: two-sided #方差齐性检验：> var.test(a,b) F test to compare two variancesdata: a and b F = 1.9646, num df = 11

18、, denom df = 9, p-value = 0.32alternative hypothesis: true ratio of variances is not equal to 1 95 percent confidence interval: 0.5021943 7.0488630 sample estimates:ratio of variances 1.964622 #可认为a和b的方差相同。#选用方差相同模型t检验：> t.test(a,b,var.equal=TRUE) Two Sample t-testdata: a and b t = -8.8148, df =

19、20, p-value = 2.524e-08alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: -48.24975 -29.78358 sample estimates:mean of x mean of y 125.5833 164.6000 p-value = 2.524e-08<0.05，因而认为两者有显著差别。教师签名年 月 日开课学院、实验室： 数学与统计学院 实验时间 ： 2013 年 月 日实验项目名 称 回归分析实验项目类型验

20、证演示综合设计其他指导教师张应应成 绩一、实验目的1.掌握回归分析的原理以及在R软件上的使用2.掌握回归诊断在R软件上的使用3.掌握回归诊断的原理二、实验内容 为了估计山上积雪融化后对下游灌溉的影响，在山上建立一个 观测站，测量最大面积积雪深度X与当年灌溉面积Y，测得连续10年的数据10年中最大积雪深度与当年灌溉面积的数据年序X/mY/hm2序号X/mY/hm215.1190767.8300023.5128774.5194737.1270085.6227346.2237398.0311358.83260106.424931.试画相应的散点图，判断Y与X是否有线性关系2.求出Y关于X的一元线性回

21、归方程；3.对方程做显著性分析；4.先测得今年的数据是X=7m，给出今年灌溉面积的预测值和相应的区间估计（=0.05）三、实验原理、方法（算法）、步骤 回归分析研究的主要问题是： 1.确定Y与X1，X2.Xp间的定量关系表达式，即回归方程 2.对求得的回归方程的可信度进行检验 3.判断自变量Xj(j=1,2,.p)对Y有无影响 4.利用所求的的回归方程进行预测和控制，在比较严密的分析中，一般不轻易剔除变量，而是对显著差异变量进行分析。四、实验环境（所用软件、硬件等）及实验数据文件 软件： R 2.15.3 数据见实验内容 五、实验结果及实例分析#输入数据并运行得：x<-c(5.1,3.5

22、,7.1,6.2,8.8,7.8,4.5,5.6,8.0,6.4)y<-c(1907,1287,2700,2373,3260,3000,1947,2273,3113,2493)plot(x,y)分析结果：由散点图可得x,y线性相关lm.sol<-lm(y1+x)summary(lm.sol)Call:lm(formula = y 1 + x)Residuals: Min 1Q Median 3Q Max -128.591 -70.978 -3.727 49.263 167.228 Coefficients: Estimate Std. Error t value Pr(>|t

23、|) (Intercept) 140.95 125.11 1.127 0.293 x 364.18 19.26 18.908 6.33e-08 *-Signif. codes: 0 * 0.001 * 0.01 * 0.05 . 0.1 1 Residual standard error: 96.42 on 8 degrees of freedomMultiple R-squared: 0.9781, Adjusted R-squared: 0.9754 F-statistic: 357.5 on 1 and 8 DF, p-value: 6.33e-08分析结果：由上述结果可得y关于x的一元

24、线性回归方程为y=140.95+364.18x；并由F检验和t检验，可得回归方程通过了回归方程的显著性检验#对数据进行预测，并且给相应的区间估计new<-data.frame(x=7)lm.pred<-predict(lm.sol,new,interval="prediction",level=0.95)lm.pred fit lwr upr1 2690.227 2454.971 2925.484分析结果：预测值为2690.227，估计区间为2454.971 ,2925.484教师签名年 月 日开课学院、实验室：数学与统计 实验时间 ：2013年 04月 20日

25、实验项目名 称判别分析实验项目类型验证演示综合设计其他指导教师张应应成 绩一、实验目的1.在已知样本有多少类的情况下对样本进行分类；2.熟悉R软件的各种操作；3.利用R软件进行各种方式的判别分析。二、实验内容 根据经验 今天与昨天的温度差X1及今天的压温差（气压与温度之差）X2是预报明天下雨或不下雨的两个重要因素。现有一批已收集的数据资料，如表所示，金测得x1=8.1，x2=2.0试问预报明天下雨还是预报明天不下雨？分别用距离判别、Bayes判别、（考虑方差相同与方差不同两种情况）和Fisher判别来得到你所需要的结论。具体数据见统计建模与R软件P420表8.7.三、实验原理、方法（算法）、步

26、骤距离判别法：，判断其值是否大于0.Bayes判别法： ，判断值与的关系；Fisher判别法： 4、 实验环境（所用软件、硬件等）及实验数据文件 实验数据见实验内容 软件：R2.15.1五、实验结果及实例分析：# 距离判别classX1=data.frame(x1=c(-1.9,-6.9,5.2,5.0,7.3,6.8,0.9,-12.5,1.5,3.8),x2=c(3.2,10.4,2.0,2.5,0.0,12.7,-15.4,-2.5,1.3,6.8)classX2=data.frame(x1=c(0.2,-0.1,0.4,2.7,2.1,-4.6,-1.7,-2.6,2.6,-2.8),

27、x2=c(0.2,7.5,14.6,8.3,0.8,4.3,10.9,13.1,12.8,10.0)TstX=data.frame(x1=c(8.1),x2=c(2.0)# 对训练样本的回代情况# var.equal=T：有4个错判，错判率为4/20=0.2# var.equal=F：有5个错判，错判率为5/20=0.25source("discriminiant.distance.R")discriminiant.distance(classX1,classX2,var.equal=T)discriminiant.distance(classX1,classX2,var.

28、equal=F)# 对测试样本进行判别：均判为第1组discriminiant.distance(classX1,classX2,TstX,var.equal=T)discriminiant.distance(classX1,classX2,TstX,var.equal=F)# Bayes# TrnX1, TrnX2以矩阵的形式输入TrnX1=matrix(c(-1.9,-6.9,5.2,5.0,7.3,6.8,0.9,-12.5,1.5,3.8,3.2,10.4,2.0,2.5,0.0,12.7,-15.4,-2.5,1.3,6.8),ncol=2)TrnX2=matrix(c(0.2,-0

29、.1,0.4,2.7,2.1,-4.6,-1.7,-2.6,2.6,-2.8,0.2,7.5,14.6,8.3,0.8,4.3,10.9,13.1,12.8,10.0),ncol=2)TstX=data.frame(x1=c(8.1),x2=c(2.0)# 对训练样本的回代情况# var.equal=T：有4个错判，错判率为4/20=0.2# var.equal=F：有5个错判，错判率为5/20=0.25source("discriminiant.bayes.R")discriminiant.bayes(TrnX1,TrnX2,rate=1,var.equal=T)disc

30、riminiant.bayes(TrnX1,TrnX2,rate=1,var.equal=F)# 对测试样本进行判别：均判为第1组discriminiant.bayes(TrnX1,TrnX2,rate=1,TstX,var.equal=T)discriminiant.bayes(TrnX1,TrnX2,rate=1,TstX,var.equal=F)discriminiant.bayes(classX1,classX2,rate=1,TstX,var.equal=T)discriminiant.bayes(classX1,classX2,rate=1,TstX,var.equal=F)# F

31、isher# 对训练样本的回代情况# 有4个错判，错判率为4/20=0.2source("discriminiant.fisher.R")discriminiant.fisher(classX1,classX2)# 对测试样本进行判别：判为第1组discriminiant.fisher(classX1,classX2,TstX)# 三种方法均预报明天下雨colMeans(classX1)colMeans(classX2)x1 x2 -0.38 8.25 教师签名年 月 日开课学院、实验室： 数学与统计学院 实验时间 ： 2013年 4月 日实验项目名 称主成分分析实验项目类

32、型验证演示综合设计其他指导教师张应应成 绩一、实验目的1.掌握利用主成分分析的理论思想，将多指标化为少数几个综合的指标；2.熟练应用R软件来进行主成分分析；3.利用主成分分析的结果对具体问题作出合理的解释分析。二、实验内容 用主成分法探讨城市工业主体，表格里是某市工业部门13个行业，分别是冶金，电力，煤炭，化学，机械，建材，食品，森工，纺织，缝纫，皮革，造纸，文教艺术品，8个指标，分别是年末固定资产净值X1，职工人数X2，工业总产值X3，全员劳动生产率X4，百元固定原值实现产值X5，资金利税率X6，标准燃料消费量X7,和能源利用开发效果X8.的数据。 见统计建模与R软件 P475 表9.7(1

33、) 利用主成分分析方法分析确定8个指标的几个主成分，并对主成分进行解释；(2) 利用主成分得分对13个行业内进行排序和分类。3、 实验原理、方法（算法）、步骤基本原理：主成分分析是数学上对数据降维的一种方法。其基本思想是设法将原来众多的具有一定相关性的指标X1，X2，XP（比如p个指标），重新组合成一组较少个数的互不相关的综合指标Fm来代替原来指标。计算步骤：（1） 计算协方差矩阵；（2） 求出的特征值i及相应的正交化单位特征向量；（3） 选择主成分；（4） 选择主成分；（5）选择主成分.四、实验环境（所用软件、硬件等）及实验数据文件 数据见实验内容 R 2.15.35、 实验结果及实例分析（

34、1）利用主成分确定了8个指标的主成分，有4个，即主成分碎石图所示> industry<-data.frame(+X1=c(90342,4903,6735,49454,139190,12215,2372,11062,17111,1206,2150,5251,14341)，+X2=c(52455,1973,21139,36241,203505,16219,6572,23078,23907,3930,5704,6155,13203)，+X3=c(101091,2035,3767,81557,215898,10351,8103,54935,52108,6126,6200,10383,193

35、96),+X4=c(19272,10313,1780,22504,10609,6382,12329,23804,21796,15586,10870,16875,14691),+ X5=c(82.0,34.2,36.1,98.1,93.2,62.5,184.4,370.4,221.5,330.4,184.2,146.4,94.6),+ X6=c(16.1,7.1,8.2,25.9,12.6,8.7,22.2,41.0,21.5,29.5,12.0,27.5,17.8),+X7=c(197435,592077,726396,348226,139572,145818,20921,65486,6380

36、6,1840,8913,78796,6354),+X8=c(0.172,0.003,0.003,0.985,0.628,0.066,0.152,0.263,0.276,0.437,0.274,0.151,1.574) )> industry.pr<-princomp(industry,cor=T) > summary(industry.pr) #做主成分分析，得到4个主成分，累积贡献率达94.68%Importance of components: Comp.1 Comp.2 Comp.3 Comp.4 Comp.5Standard deviation 1.7620762 1

37、.7021873 0.9644768 0.80132532 0.55143824Proportion of Variance 0.3881141 0.3621802 0.1162769 0.08026528 0.03801052Cumulative Proportion 0.3881141 0.7502943 0.8665712 0.94683649 0.98484701 Comp.6 Comp.7 Comp.8Standard deviation 0.29427497 0.179400062 0.0494143207Proportion of Variance 0.01082472 0.00

38、4023048 0.0003052219Cumulative Proportion 0.99567173 0.999694778 1.0000000000> load<-loadings(industry.pr) #求出载荷矩阵> loadLoadings: Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6 Comp.7 Comp.8X1 -0.477 -0.296 -0.104 0.184 0.758 0.245X2 -0.473 -0.278 -0.163 -0.174 -0.305 -0.518 0.527X3 -0.424 -0.37

39、8 -0.156 -0.174 -0.781X4 0.213 -0.451 0.516 0.539 0.288 -0.249 0.220X5 0.388 -0.331 -0.321 -0.199 -0.450 0.582 0.233 X6 0.352 -0.403 -0.145 0.279 -0.317 -0.714 X7 -0.215 0.377 -0.140 0.758 -0.418 0.194 X8 -0.273 0.891 -0.322 0.122 Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6 Comp.7 Comp.8SS loadings 1.

40、000 1.000 1.000 1.000 1.000 1.000 1.000 1.000Proportion Var 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125Cumulative Var 0.125 0.250 0.375 0.500 0.625 0.750 0.875 1.000> plot(load,1:2) > text(load,1,load,2,adj=c(-0.4,-0.3)> screeplot(industry.pr,npcs=4,type="lines") #得出主成分的碎石图> biplot(industry.pr) #得出在第一，第二主成分之下的散点图> p<-predict(industry.pr) #预测数据，讲预测值放入p中&