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

- -话说直接copy是不道德的哈

开课学院、实验室： 数学与统计学院 实验时间 ： 2013 年 3 月 日

- -话说直接copy是不道德的哈

> n<-5;x<-array(0,dim=c(n,n)) > for (i in 1:n){ + for (j in 1:n){ + x[i,j]<-1/(i+j-1) +} + };x [,1] [,2] [,3] [,4] [,5] [1,] 1.0000000 0.5000000 0.3333333 0.2500000 0.2000000 [2,] 0.5000000 0.3333333 0.2500000 0.2000000 0.1666667 [3,] 0.3333333 0.2500000 0.2000000 0.1666667 0.1428571 [4,] 0.2500000 0.2000000 0.1666667 0.1428571 0.1250000 [5,] 0.2000000 0.1666667 0.1428571 0.1250000 0.1111111 > det(x) [1] 3.749295e-12 > solve(x) [,1] [,2] [,3] [,4] [,5] [1,] 25 -300 1050 -1400 630 [2,] -300 4800 -18900 26880 -12600 [3,] 1050 -18900 79380 -117600 56700 [4,] -1400 26880 -117600 179200 -88200 [5,] 630 -12600 56700 -88200 44100 > eigen(x) $values [1] 1.567051e+00 2.085342e-01 1.140749e-02 3.058980e-04 3.287929e-06 $vectors [,1] [,2] [,3] [,4] [,5] [1,] 0.7678547 0.6018715 -0.2142136 0.04716181 0.006173863 [2,] 0.4457911 -0.2759134 0.7241021 -0.43266733 -0.116692747 [3,] 0.3215783 -0.4248766 0.1204533 0.66735044 0.506163658 [4,] 0.2534389 -0.4439030 -0.3095740 0.23302452 -0.767191193 [5,] 0.2098226 -0.4290134 -0.5651934 -0.55759995 0.376245545

分析：从实验结果来看。R 软件在处理数据上相当准确，方便。

教师签名 年 月 日

- -话说直接copy是不道德的哈

开课学院、实验室：数学与统计学院 实验时间 ： 2013 年 3 月 日

- -话说直接copy是不道德的哈

10 11 12 13 14 15 16 17 18 19

Alfred Duke Guido James

M 14 M 14 M 15 M 12

69.0 112.5 63.5 102.5 67.0 133.0 57.3 83.0 62.5 84.0 59.0 99.5 72.0 150.0 64.8 128.0 57.5 85.0 66.5 112.0

Jeffrey M 13 John Philip Robert Thomas M 12 M 16 M 12 M 11

William M 15

五、实验结果及实例分析student<-read.table("3.7 数据.txt",header=T) attach(student) > cor.test(身高,体重) #Pearson 相关性检验

Pearson's product-moment correlation

data: 身高 and 体重 t = 7.5549, df = 17, p-value = 7.887e-07 alternative 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,拒绝原假设，所以身高与体重相关

教师签名 年 月 日

- -话说直接copy是不道德的哈

开课学院、实验室： 数学与统计学院 实验时间 ： 2013年 月 日

- -话说直接copy是不道德的哈

S S ta (n 1) , X ta (n 1) n n 四、实验环境(所用软件、硬件等)及实验数据文件

X

数据见实验内容，所用软件：R2.15.1 五、实验结果及实例分析 在 R 软件中运行代码：> x<-c(54,67,68,78,70,66,67,70,65,69) > t.test(x) #做单样本正态分布区间估计 One Sample t-test data: x t = 35.947, df = 9, p-value = 4.938e-11 alternative 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-test data: x t = -2.4534, df = 9, p-value = 0.01828 alternative 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,拒绝原假设，平均脉搏低于常人。

教师签名 年 月 日

- -话说直接copy是不道德的哈

开课学院、实验室： 实验时间 ： 2013 年 月 日

- -话说直接copy是不道德的哈

2 2 方差 12 2

未知， S12 和 S2 分别是 X 和 Y 的样本方差。由统计知识可知，当 H 0 为真时， X Y T ~ t (n1 n2 2) 1 ( n1 1) S 2 (n 1) S 2 1 1 2 2 S n 其中Sww n 1 2 n n 2 1 2 因此，当 T 满足(成为拒绝域) : 双边检验： T t (n1 n2 2)2

单边检验 I 单边检验 II

T t (n1 n2 2) T t (n1 n2 2)

则认为 H 0 不成立，此方法也称为 t 检验法。 四、实验环境(所用软件、硬件等)及实验数据文件 见实验内容 软件：R2.15.3 R 软件。 五、实验结果及实例分析> 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 test

data: a D = 0.1464, p-value = 0.9266 alternative hypothesis: two-sided

> ks.test(b,"pnorm",mean(b),sd(b))

One-sample Kolmogorov-Smirnov test

data: b D = 0.2222, p-value = 0.707 alternative hypothesis: two-sided ####方差齐性检验： > var.test(a,b)

- -话说直接copy是不道德的哈

F test to compare two variances

data: a and b F = 1.9646, num df = 11, denom df = 9, p-value = 0.32 alternative 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-test

data: a and b t = -8.8148, df = 20, p-value = 2.524e-08 alternative 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,因而认为两者有显著差别。

教师签名 年 月 日

- -话说直接copy是不道德的哈

开课学院、实验室： 数学与统计学院 实验时间 ： 2013 年 月 日

- -话说直接copy是不道德的哈

五、实验结果及实例分析####输入数据并运行得： x<-c(5.1,3.5,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(y~1+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|) (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 freedom Multiple 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 的一元线性回归方程为 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 upr 1 2690.227 2454.971 2925.484 分析结果：预测值为 2690.227,估计区间为[2454.971 ,2925.484]

教师签名 年 月 日

- -话说直接copy是不道德的哈

开课学院、实验室：数学与统计 实验时间 ：2013年 04月 20日

- -话说直接copy是不道德的哈

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.25 source("discriminiant.distance.R") discriminiant.distance(classX1,classX2,var.equal=T) discriminiant.distance(classX1,classX2,var.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.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)

- -话说直接copy是不道德的哈

)

## 对训练样本的回代情况 ## var.equal=T:有 4 个错判，错判率为 4/20=0.2 ## var.equal=F:有 5 个错判，错判率为 5/20=0.25 source("discriminiant.bayes.R") discriminiant.bayes(TrnX1,TrnX2,rate=1,var.equal=T) discriminiant.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)

## Fisher ## 对训练样本的回代情况 ## 有 4 个错判，错判率为 4/20=0.2 source("discriminiant.fisher.R") discriminiant.fisher(classX1,classX2)

## 对测试样本进行判别：判为第 1 组 discriminiant.fisher(classX1,classX2,TstX)

## 三种方法均预报明天下雨 colMeans(classX1) colMeans(classX2) x1 x2

-0.38 8.25 教师签名 年 月 日

- -话说直接copy是不道德的哈

开课学院、实验室： 数学与统计学院 实验时间 ： 2013年 4月 日

- -话说直接copy是不道德的哈

五、实验结果及实例分析(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,19396), +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,63806,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.5 Standard deviation 1.7620762 1.7021873 0.9644768 0.80132532 0.55143824 Proportion of Variance 0.3881141 0.3621802 0.1162769 0.08026528 0.03801052 Cumulative Proportion 0.3881141 0.7502943 0.8665712 0.94683649 0.98484701 Comp.6 Comp.7 Comp.8 Standard deviation 0.29427497 0.179400062 0.0494143207 Proportion of Variance 0.01082472 0.004023048 0.0003052219 Cumulative Proportion 0.99567173 0.999694778 1.0000000000 > load<-loadings(industry.pr) ####求出载荷矩阵 > load Loadings: Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6 Comp.7 Comp.8 X1 -0.477 -0.296 -0.104 0.184 0.758 0.245 X2 -0.473 -0.278 -0.163 -0.174 -0.305 -0.518 0.527 X3 -0.424 -0.378 -0.156 -0.174 -0.781 X4 0.213 -0.451 0.516 0.539 0.288 -0.249 0.220 X5 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.8 SS loadings 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 Proportion Var 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125 Cumulative 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") ####得出主成分的碎石图

- -话说直接copy是不道德的哈

> biplot(industry.pr)

####得出在第一，第二主成分之下的散点图

> p<-predict(industry.pr) ####预测数据，讲预测值放入 p 中 > order(p[,1]);order(p[,2]);order(p[,3]);order(p[,4]); ####将预测值分别以第一，第二，第三，第四主成分进行排序 [1] 5 1 3 2 4 6 13 11 9 7 12 10 8 [1] 5 8 4 9 10 1 13 12 7 11 6 2 3 [1] 8 1 5 3 9 12 7 10 2 6 11 4 13 [1] 11 6 5 7 10 13 12 9 1 8 3 2 4

- -话说直接copy是不道德的哈

> kmeans(scale(p),4) ####将预测值进行标准化，并分为 4 类 K-means clustering with 4 clusters of sizes 5, 1, 4, 3 Cluster means: Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6 1 0.5132590 -0.03438438 -0.3405983 -0.5130031 0.2355151 0.22441040 2 -2.5699693 -1.32913757 -0.4848689 -0.9460127 -0.9000187 -0.06497950 3 0.2381581 0.72871986 -0.2995918 0.3126036 -0.4744091 -0.19709710 4 -0.3163193 -0.47127333 1.1287426 0.7535380 0.5400265 -0.08956137 Comp.7 Comp.8 1 -0.38197798 -0.7474855 2 -0.67500209 0.4569548 3 0.09063069 0.9826915 4 0.74078975 -0.2167643 Clustering vector: [1] 4 3 3 4 2 1 1 1 1 3 1 3 4 Within cluster sum of squares by cluster: [1] 19.41137 0.00000 24.49504 16.61172 (between_SS / total_SS = 37.0 %) Available components: [1] "cluster" [6] "betweenss" "centers" "size" "totss" "withinss" "tot.withinss"

#######用 order()分别对 4 个主成分的预测值进行排序，结果是如下表(26) ,而利用 kmeans()进行动态排序得到如下分类： 第 1 类：建材(6) ,森工(7) ,食品(8) ,纺织(9) ,皮革(11) ; 第 2 类：机械(5) ; 第 3 类：电力(2) ,煤炭(3) ,缝纫(10)造纸(12) ; 第 4 类：冶金(1)化学(4) ,文教艺术用品(13) 。 成分 第一主成分： 第二主成分： 第三主成分： 第四主成分： 5 5 8 11 1 8 1 6 3 4 5 5 2 9 3 7 4 10 9 10 13 个行业排序结果 6 1 12 13 13 13 7 12 11 12 10 9 9 7 2 1 7 11 6 8 12 6 11 3 10 2 4 2 8 3 13 4

表(26)各行业按主成分得分进行排序结果

教师签名 年 月 日

本文来源：https://www.bwwdw.com/article/kr54.html