应用统计分析实验R软件2013ppt课件

1、应用统计分析实验 R软件,SPSS：这是一个很受欢迎的统计软件 容易操作， 输出漂亮， 功能齐全， 价格合理。 对于非统计工作者是很好的选择,SAS：这是功能非常齐全的软件； 美国政府政策倾斜(“权威性”) 许多美国公司使用。 价格不菲,每年交费.即使赠送,条件苛刻 尽管现在已经尽量“傻瓜化”，仍然需要一定的训练才可以进入,R软件：免费的,志愿者管理的软件。 编程方便，语言灵活，图形功能强大 有不断加入的各个方向统计学家编写的统计软件包。也可以自己加入自己算法的软件包. 这是发展最快的软件，受到世界上统计师生的欢迎。是用户量增加最快的统计软件。 对于一般非统计工作者来说，主要问题是它没有“傻瓜

2、化,Minitab：这个软件是很方便的功能强大而又齐全的软件，也已经“傻瓜化”，在我国用的不如SPSS与SAS那么普遍。 Eviews：这是一个主要处理回归和时间序列的软件。 GAUSS：这是一个很好用的统计软件，许多搞经济的喜欢它。主要也是编程功能强大。目前在我国使用的人不多。 MATLAB：这也是应用于各个领域的以编程为主的软件，在工程上应用广泛。但是统计方法不多,R的历史,S语言在1980年代后期在AT它是一个由志愿者组成的工作努力的国际团队,下载R软件 ,学习网站 /pages/newhtm/r/s

3、chtml,基本语法 向量 矩阵 list与data.frame 读写数据文件 控制语句与自定义函数,一. R软件的使用,基本语法,1. 变量使用即定义: 变量名区分大小写, 也可用中文命名 变量赋值可采用4种形式：=，, assign() 变量类型自动由变量赋值确定,a=10; bc ; assign(“d”,40) 中国=“中华人民共和国” #生成字符串变量,2.注释符号 # 语句连接符,3. 算术运算符: +,-,*,/,(乘方),% (模), %/% (整除,4.常用的数学函数有:abs , sign , log , log2, log10 , sqrt , exp , sin , c

4、os , tan , acos , asin, atan , cosh , sinh， tanh,5. 查看帮助 help(round) ?abs,向量,向量的赋值(一维数组, 下标从1开始) a=c(d1,d2,d3,) 间隔为1的等差序列: a:b 指定间隔的等差序列: seq(from,to,by) seq(length, from, by) 重复序列： rep(vec, times) rep(vec,times,len,each,a=c(3,5,8,10); b=1:10; c=seq(1,10,2); d=seq(-pi,pi, 0.2) e=rep(a,3); f=rep(a, 2

5、, each=3,随机向量 rnorm(10) #10个服从标准正态分布的随机数,a=1:5 (b=1:5) #同上，只不过显示出来 a2 #取出a中第二个元素 ac(2,4)=c(4,8) #修改a中第2、4个元素分别为4、8 a-5 #扣除第5个元素取出来 a3 #判断a中元素是否小于3 1 TRUE TRUE FALSE FALSE FALSE aa3 # 取出a中小于3的元素 a6=12 a=a-c(1,3,5) #去掉第1、3、5元素,2.向量的下标运算,3.向量的长度,length(a,matrix(data=NA, nrow=1, ncol=1, byrow=FALSE,A=ma

6、trix(1:10, 2,5) B=matrix(1:10,2,5,byrow=TRUE) #按行放置元素,注意：默认是按列放置元素,1 ,2 ,3 ,4 ,5 1, 1 2 3 4 5 2, 6 7 8 9 10,1 ,2 ,3 ,4 ,5 1, 1 3 5 7 9 2, 2 4 6 8 10,矩阵(二维数组,x=matrix(rnorm(24),4,6) x2,1 #第2,1元素 xc(2,1), #第2和第1行 x,c(1,3) #第1和第3列 xx,10,1 #第1列大于0的元素 x,-c(1,3) #没有第1、3列的x x-2,-c(1,3) #没有第2行、第1、3列的x,1.矩阵的

7、元素访问,2. 矩阵的维数问题 dim(A) #获得维数,返回向量 nrow(A) ,ncol(A) #获得行数和列数,rownames(A), colnames(A) #访问各维名称,3. 向量和数组/矩阵的转化: 只要定义向量的维数即可实现向量和数组转化,c=1:12; a=matrix(c, nrow=2,ncol=6) dim(c)=c(3,4) b=as.vector(c) A=diag(c(1,4,5) #以向量为对角元生成对角矩阵 a=diag(A) #获取矩阵的对角元,3. 矩阵运算 +,-,*,/ 分别是矩阵内部元素的四则运算 向量矩阵间：向量按列匹配与矩阵运算 例如： A=

8、matrix(1:6, nrow=3) B=matrix(10:15,nrow=3) C=c(100,200) 则： A+B A*B A+C,11 17 10 52 101 204 13 19 22 70 202 105 15 21 36 90 103 206,A %*% B #乘法 t(mat) #转置 det(mat)#行列式 solve(mat)#逆矩阵 eigen(mat) #求特征值与特征向量,4.多维数组,a=array(data=1:24,dim=c(2,3,4) a,1,列表和数据框,1. 列表是一种特殊的对象集合，各元素类型任意 生成： list(name1=value1,

9、, namen=valuen) 访问/修改：对象名下标 或 对象名$namei,stu=list(age=10, name=“Tom”, interests=c(“swimming”,”drawing”) stu2 stu$name stu$name=“john” names(stu) #得到所有的对象名,2. 数据框是R的一种数据结构，以矩阵形式保存数据各列类型可以不同，每列为一变量，每行为样品 各列长度相等 data.frame(,stu=data.frame(name=c(Tom, Rose), age=c(30,32) names(stu) #得到所有的变量名 colnames(stu

10、) #列名，同上 rownames(stu) #得到行名 attach(x) #把数据框中的变量链接到内存中 x=data.frame(matrix(1:6,nrow=2) #矩阵转化为数据框,x=data.frame(id=101:120,score=round(rnorm(20,70,10) #取出前两行数据 x1:2, #选出score60的数据 xx260,读写数据文件,读表格形式文本文件 read.table(file, header=FALSE,sep=“ “,) header表示第一行是否有变量名，sep表示数据分割的字符,D1=read.table(“e:test1.txt”,

11、header=TRUE) D2=read.table(“clipboard”) #从剪贴板中获得文本数据,2. 计算的数据保存成文本文件或CSV文本文件 write.table(data, file=“”) #纯文本格式,write.table(D1,file=“e:file_1.txt”,1. 分支结构 单分支：if(条件) 语句 if(条件) 复合语句; 双分支：if(条件) 语句1 else 语句2 多分支：if() 语句1 else if () 语句2 else if() 语句n else 语句n+1,x=1; if(x0) a=10;b=10 else a=20; b=20,控制语句

12、控制语句与自定义函数,2. 循环结构 for (name in express) expr; while(condition) expr; repeat exprs; if(达到中止条件) break,例子：计算11+22+33+1010 i=1;j=1; k=1; s1=0; s2=0; s3=0; for( i in 1:10) s1=s1+ii; #用for while(j10) break,fun-function(arg1,arg2,) 注意：最后一句表达式的值为返回值,myfun=function(k ) sum=0; for(i in 1:k) sum=sum+ii sum,自定

13、义函数,统计量 数据中心化与标准化 和分布有关 数据的图形表示 统计图形,二. 数据描述性统计,统计量,sum(),max(),min() mean() #平均值 median() #中位数 var() #方差 sd() #标准差 cov() #协方差阵，参数为矩阵或数据框 cor() #相关系数，参数为矩阵或数据框 summary() table() #列联表，统计不同值出现的个数,x = c(10,90,71:78) mean(x) mean(x, trim = 0.10,A = cbind(c(90,60,70),c(85,66,95) apply(A,1,mean) # 对每一行取平均

14、 apply(A,2,mean) #对每一列取平均 cov(A) cor(A,数据中心化与标准化,scale(x, center = TRUE, scale = TRUE,dnorm(x, mean = 0, sd = 1) #计算正态分布的密度函数 pnorm(p, mean = 0, sd = 1) #计算正态分布的分布函数 qnorm(q, mean = 0, sd = 1) #计算正态分布的分位数 rnorm(n, mean = 0, sd = 1) #计算n个正态分布的随机数 其它分布的关键词： unif , exp,chisq, t, f, binom,pois,得到参数为1的指数

15、分布在2的密度函数值 dexp(2,1 ) #产生5个均匀分布U(0,10)的随机数 runif(5,0,10,和分布相关,抽样样本,x=sample(1:20,10) #从1，20中随机地不放回抽取10个值作为样本 y=sample(1:20,10, replace =TRUE) #从1，20中随机地有放回抽取10个值作为样本,数据的图形表示,x=seq(-10,10,2); y=rnorm(11) plot(x, y) plot(x, y,main=散点图,xlab=横坐标x,ylab=纵坐标y) text(x,y, 1:length(x),-1) # 写入点序号,第三个参数默认如此 pl

16、ot(x,y,type = o, col = red,pch=2,lty=3,1.散点图 plot(x,y,type=“p” #散点图，默认 type=“l #绘实线 type=“o” #实线通过所有的点 type=“n” #不画点 参数lty表示线的类型 0=blank, 1=solid (default), 2=dashed, 3=dotted, 4=dotdash, 5=longdash, 6=twodash,低水平作图函数有 points(x=,y=)#增加点（x,y) lines(x=,y=) #增加（x,y)连成的线 abline( a=,b=) #增加线y=a+bx text(x

17、=, y=, labels = ) #增加文本 legend(x, y = , legend, col ) #增加标注,x=seq(-10,10,2); y=rnorm(11) par(mfrow=c(3,1) #准备画3个图 plot(x, y) plot(x, y,main=散点图,xlab=横坐标x,ylab=纵坐标y) text(x,y, 1:length(x),-1) # 写入点序号 plot(x,y,type = o, col = red,pch=2,lty=3,2.散点图矩阵 pairs(x,A=matrix(rnorm(72,70,20),nrow=24) pairs(A,3.

18、 多组散点图,matplot(x,y,x=seq(0,2*pi,0.2) matplot(x,cbind(sin(x),cos(x),type=c(p, l),col=2:3,4.星座图 5.脸谱图,library(aplpack) faces(x,stars(x,例1. 在同一图中绘制标准正态分布和自由度为5的t分布的概率密度图,x=seq(-4,4,0.01) plot(x,dnorm(x),type=l,ylab=f(x),main=概率密度图) lines(x,dt(x,5),col=red,lty=2) legend(2,0.4,c(N(0,1),t(5),lty=1:2,col=c

19、(black,red,x=seq(-4,4,0.01) matplot(x,cbind(dnorm(x),dt(x,5),type=l,ylab=f(x),main=概率密度图) legend(2,0.4,c(N(0,1),t(5),lty=1:2,col=1:2,lwl.andrews=function(X) t=seq(-pi,pi,0.2) n=nrow(X) f=matrix(0,nrow=length(t),ncol=n) for(i in 1:n) f,i=Xi,1/sqrt(2) for(j in 2:ncol(X) if (j%2=0) f,i=f,i+Xi,j*sin(j/2

20、*t) else f,i=f,i+Xi,j*cos(j%/%2*t) matplot(t,f,col=1:n,type=o,main=调和曲线图),例2：绘制调和曲线图,实验一 R软件的使用,求矩阵B的行列式、逆矩阵、特征根和特征向量及ABA,实验目的：熟练使用R软件，掌握对向量、矩阵和数据框的操作和对数据的描述性统计,学号 性别 数学 语文 英语 物理 1 3001 男 80 76 81 90 2 3002 男 55 67 79 63 3 3003 女 27 52 53 43 4 3004 女 62 57 66 69 5 3005 女 94 80 76 71 6 3006 男 97 59 6

21、1 74 7 3007 男 79 62 78 88 8 3008 男 70 92 78 83 9 3009 女 67 66 65 55 10 3010 男 97 46 55 75 11 3011 男 84 82 88 92 12 3012 男 40 57 67 53 13 3013 男 72 70 80 84 14 3014 男 43 31 56 60 15 3015 男 95 75 58 74 16 3016 男 62 85 91 83 17 3017 女 71 76 68 74 18 3018 男 80 40 52 71 19 3019 女 96 75 77 80 20 3020 女 82

22、 51 73 77 21 3021 女 70 73 85 77 22 3022 男 61 40 47 56 23 3023 女 64 75 66 71 24 3024 女 60 21 48 65 25 3025 男 81 88 80 84 26 3026 女 89 54 65 71 27 3027 男 55 38 58 62 28 3028 女 60 54 68 62 29 3029 男 65 63 65 75 30 3030 男 69 55 73 75 31 3031 男 80 29 61 68 32 3032 男 91 59 73 85 33 3033 男 69 73 75 78 34 3

23、034 女 90 88 73 87 35 3035 女 64 70 72 77,2. 某班35人期末成绩见右图 绘制四科成绩的散点图矩阵、雷达图和脸谱图 绘制前5名学生各科成绩的轮廓图 绘制数学与物理成绩的散点图，其中男生与女生用不同的颜色和形状的点表示,计算四科成绩的平均分和相关矩阵 分性别统计四科成绩的平均分,例3. 在同一图中绘制二项分布B(10,0.2)、泊松分布P(2)的概率分布图形,三. 回归分析,与线性回归相关常用函数全局挑选最佳自变量组合逐步回归,lm( ) #建立线性回归模型 summary( ) #提取模型资料 predict( ) #用模型作预测 plot() #绘制模型

24、诊断的几种图形,与线性回归相关常用函数,lm(yx) # y为向量, x可为向量与矩阵 lm(yx+0) #回归方程没有常数项 lm(yx1+x2+I(x12) #自变量有三项x1,x2,x12,lm(yx1+x2,data=X) # X为数据框, y,x1,x2为其中变量名 lm(y.,data=X) #自变量为其它所有变量 lm(y.-x2,data=X) #自变量为除x2的其它所有变量,例1：某气象站收集了15年关于年初的最高温度X(从1月11日算起)与秋季第一次强冷空气出现日期Y(从9月11日算起)的数据,flm=lm(yx) summary(flm) plot(x,y,main=一元

25、回归图) abline(flm) points(x,flm$fitted,col=red,pch=3,常数项,回归系数,回归系数的标准差,t检验值,t检验P值。标注*，则说明此项回归系数显著不为零，即有意义,模型残差的标准差的估计量,决定系数R2,修正后的决定系数,回归方程的F检验统计量值及检验P值,Residuals: Min 1Q Median 3Q Max -7.3217 -1.6211 -0.1418 1.5485 6.0382 Coefficients: Estimate Std. Error t value Pr(|t|) (Intercept) 3.36014 2.06876 1

26、.624 0.128 x 0.82006 0.08895 9.219 4.58e-07 * - Signif. codes: 0 * 0.001 * 0.01 * 0.05 . 0.1 1 Residual standard error: 3.243 on 13 degrees of freedom Multiple R-squared: 0.8673, Adjusted R-squared: 0.8571 F-statistic: 84.99 on 1 and 13 DF, p-value: 4.584e-07,若需编程，很多信息存放在lm和summary对象中 flm$coef #得到回归

27、系数 flm$resi #得到残差向量 flm$fitt #得到拟合向量 summary(flm)$sigma #模型残差的标准差的估计量 summary(flm)$r.squared #决定系数R2 summary(flm)$adj.r.squared #修正后的决定系数,预测，注意第二个参数为数据框类型，成员名称必须与原始的自变量一致 predict(flm,data.frame(x=c(14,26) predict(flm,data.frame(x=c(14,26),interval=prediction) #区间预测 fit lwr upr 1 14.84093 7.472051 22

28、.2098 2 24.68160 17.389301 31.9739,例2:回归模拟,x=runif(30,0,10) #均匀分布的随机数 x=sort(x) y=3*x-0.2*x2+rnorm(x) # y与x的准确关系 plot(x,y) l1=lm(yx);l2=lm(yx+0); l3=lm(yx+I(x2); l4=lm(yx+I(x2)+0) abline(l1,col=red,lty=1); abline(a=0,b=l2$coef,col=green,lty=2) lines(x,l3$fit,col=blue,lty=3) lines(x,l4$fit,col=orange

29、,lty=4) legend(bottomright, legend=c(线性, 无常数项线性, 二次, 二次无常数项),col=c(red, green, blue, orange),lty=1:4,实验二：最小二乘估计的求解及检验,全局挑选最佳自变量组合,library(leaps) r=summary(regsubsets(y.,data=X) n=nrow(X) aic=n*log(r$rss/n)+2*apply(r$which,1,sum) #计算AIC统计量 data.frame(r$outmat,r$rss,r$rsq,r$adjr2,r$cp,AIC=aic,例：水泥放热数据

30、挑选最佳自变量组合,x1 x2 x3 x4 y 1 7 26 6 60 78.5 2 1 29 15 52 74.3 3 11 56 8 20 104.3 4 11 31 8 47 87.6 5 7 52 6 33 95.9 6 11 55 9 22 109.2 7 3 71 17 6 102.7 8 1 31 22 44 72.5 9 2 54 18 22 93.1 10 21 47 4 26 115.9 11 1 40 23 34 83.8 12 11 66 9 12 113.3 13 10 68 8 12 109.4,残差平方和,决定系数R2,修正后的决定系数,Cp统计量,AIC统计量,

31、step(lm(y.,data=X) step(lm(y1,data=X),+x1+x2+x3+x4,逐步回归,实验三：回归分析中的自变量的挑选,四.判别分析,距离判别 Bayes判别 fisher判别,计算马氏距离 mahalanobis(x, center, cov, inverted=FALSE, .,距离判别,某总体的中心即均值向量,某总体的协方差阵,若设为TRUE,表明cov已求逆。默认为False,dda1=function(x,newx=NULL) #x 为训练样本，最后一列为类别, newx为待判样本 x=as.matrix(x); p=ncol(x)-1; #变量个数 k=m

32、ax(x,p+1); #类别个数 n=nrow(x); #已知数据的个数 if(is.null(newx)|ncol(newx)!=p) newx=x,1:p #如果待判样本newX为空，则为训练样本 m=nrow(newx); #待判数据的个数 md=matrix(-1,m,k); #距离矩阵k列，分别为到k个总体的距离 colnames(md)=paste(Dis,1:k,sep=); newClass=rep(-1,m); #新的类别 cov_each=matrix(0,nrow=p,ncol=p) ; for(i in 1:k) temp=xx,p+1=i,1:p; center=ap

33、ply(temp,2,mean); #计算当前总体的中心 cov_each=cov(temp) md,i=mahalanobis(newx,center,cov_each);# 计算待判样本到当前总体的马氏距离 newClass=apply(md,1,which.min); #距离矩阵md每一行最小的列就是归属类 list(matrix_Distance=md,newClass=newClass),假设各总体协方差阵不相等时的距离判别,例1,湿度差 压温差 q 1 -1.9 3.2 1 2 -6.9 10.4 1 3 5.2 2.0 1 4 5.0 2.5 1 5 7.3 0.0 1 6 6.

34、8 12.7 1 7 0.9 -15.4 1 8 -12.5 -2.5 1 9 1.5 1.3 1 10 3.8 6.8 1 11 0.2 6.2 2 12 -0.1 7.5 2 13 0.4 14.6 2 14 2.7 8.3 2 15 2.1 0.8 2 16 -4.6 4.3 2 17 -1.7 10.9 2 18 -2.6 13.1 2 19 2.6 12.8 2 20 -2.8 10.0 2,假设协方差阵不等时的距离判别结果,MarsDistances ,1 ,2 1, 0.23479212 2.1460083 2, 3.03833480 6.9243049 3, 0.4579097

35、3 7.3283556 4, 0.40796928 6.6425393 5, 1.15825305 13.3381431 6, 2.44342495 9.2742507 7, 5.21569514 31.8233372 8, 4.51109115 31.4452037 9, 0.02172825 3.5632768 10, 0.51080152 2.9935950 11, 0.31493244 0.4242599 12, 0.55187015 0.1087923 13, 2.70320062 1.9090285 14, 0.67515753 1.5358775 15, 0.07154715 4

36、.3721831 16, 0.90447066 4.1014802 17, 1.61226774 0.4918710 18, 2.57485181 1.7141432 19, 1.92540870 2.3430943 20, 1.56318137 0.9992699,newClass 1 1 1 1 1 1 1 1 1 1 1 2 2 1 1 1 2 2 1 2,实验目的：掌握距离判别准则和R软件编程 实验任务：编写假设各总体协方差阵相等时的距离判别算法，并对p182数据进行距离判别,实验四：距离判别,Bayes判别,Library(MASS) qda() #Quadratic discrim

37、inant analysis. #其实质是正态总体时的最大后验概率准则,Library(MASS) #默认先验概率prior 为各组数据所占比例 qd=qda(g.,data=X)# X为数据框类型，含有类别g qd=qda(x,g) #x中不含类别，g为类别,z=predict(qd, test) #test为待判样本,若省略，则对训练样本判别 z$class # 类别结果 z$posterior # 后验概率 table(X$g,z$class) #若对训练样品判别，可输出判对情况,library(MASS) qd=qda(g.,data=X) predict(qd)$class 1 2

38、1 1 1 1 1 1 1 1 2 2 2 2 2 1 1 2 2 2 2 plot(X,1:2,pch=X,3,col=X,3) text(X,1,X,2,adj=-0.5, col=X,3,湿度差 压温差 q 1 -1.9 3.2 1 2 -6.9 10.4 1 3 5.2 2.0 1 4 5.0 2.5 1 5 7.3 0.0 1 6 6.8 12.7 1 7 0.9 -15.4 1 8 -12.5 -2.5 1 9 1.5 1.3 1 10 3.8 6.8 1 11 0.2 6.2 2 12 -0.1 7.5 2 13 0.4 14.6 2 14 2.7 8.3 2 15 2.1 0.

39、8 2 16 -4.6 4.3 2 17 -1.7 10.9 2 18 -2.6 13.1 2 19 2.6 12.8 2 20 -2.8 10.0 2,舍一验证 qd=qda(g.,data=X,CV=TRUE,class 1 2 2 1 1 1 1 1 1 1 2 2 2 2 2 1 1 2 2 2 2,ld=lda(G.,data=X) ld$scaling #投影方向 z=predict(ld,test) z$class z$posterior #后验概率 z$x # X* ld$scaling中心化后的矩阵,费希尔判别 数据降维后，在各组协方差阵相等下进行的最大后验概率准则判别,湿度

40、差 压温差 q 1 -1.9 3.2 1 2 -6.9 10.4 1 3 5.2 2.0 1 4 5.0 2.5 1 5 7.3 0.0 1 6 6.8 12.7 1 7 0.9 -15.4 1 8 -12.5 -2.5 1 9 1.5 1.3 1 10 3.8 6.8 1 11 0.2 6.2 2 12 -0.1 7.5 2 13 0.4 14.6 2 14 2.7 8.3 2 15 2.1 0.8 2 16 -4.6 4.3 2 17 -1.7 10.9 2 18 -2.6 13.1 2 19 2.6 12.8 2 20 -2.8 10.0 2,ld$scaling LD1 湿度差 -0.

41、0682677 压温差 0.1562112 predict(ld)$class 1 1 2 1 1 1 2 1 1 1 1 2 2 2 2 1 2 2 2 2 2,library(MASS) class1=predict(qda(class.,data=train)$class #Bayes判别 class2=predict(lda(class.,data=train)$class #Fisher线性判别 newtrain=predict(lda(class.,data=train)$x class3=predict(qda(newtrain, train$class)$class #Fish

42、er投影后再做Bayes判别 rbind(class1, class2, class3) table(train$class, class1) table(train$class, class2) table(train$class, class3,class1 2 1 1 1 1 1 1 1 1 2 2 2 2 2 1 1 2 2 2 2 class2 1 2 1 1 1 2 1 1 1 1 2 2 2 2 1 2 2 2 2 2 class3 1 2 1 1 1 2 1 1 1 2 2 2 2 2 1 2 2 2 2 2,2281342011 2451341040 2001671227 1

43、7015078 1001672014 225125714 130100612 15011776 1201331026 160100510 185115519 17012564 16514253 135108212 10011772,例2,1 ,2 ,3 ,4 ,5 ,6 ,7 ,8 ,9 ,10 ,11 ,12 ,13 ,14 ,15 class1 1 1 1 1 1 2 2 2 2 2 3 3 3 3 3 class2 1 1 1 3 1 2 2 3 2 2 2 3 3 3 3 class3 1 1 1 1 1 2 2 3 2 2 2 3 3 3 3,实验五 判别分析实验目的：掌握Bayes判别和Fisher判别的基本原理及R软件实现实验任务,费希尔于1936年发表的鸢尾花数据是对3种鸢尾花：刚毛鸢尾花