==统计编程的框架与R语言统计分析基础.pdf

《==统计编程的框架与R语言统计分析基础.pdf》由会员分享，可在线阅读，更多相关《==统计编程的框架与R语言统计分析基础.pdf（99页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、 统计编程的框架 与 R 语言统计分析基础 （北京地区第一次COS沙龙活动讲稿）中国人民大学统计学院 谢益辉2007 年 5 月 13 日 2 目 录 一、道冲，而用之或不盈一、道冲，而用之或不盈.4 1为什么需要统计编程？.5 2哪里需要编程？.6 3从一般计算机编程说起.10 4数值计算与统计编程.12 5统计编程的流程.16 二、欲善其事，先利其器二、欲善其事，先利其器.17 1什么是利器？.18 3 2为什么用R？.19 3谁不适合用R？.20 三、一箪食，一瓢饮，在陋巷三、一箪食，一瓢饮，在陋巷.21 1获取与安装.23 2R数据结构.26 3R统计分析.55 4R图形.76 4 一

2、、道冲，而用之或不盈一、道冲，而用之或不盈 道冲，而用之或不盈。渊兮，似万物之宗。挫其锐，解其纷，和其光，同其尘。湛兮，似或存。吾不知谁之子，象帝之先。老子 5 1为什么需要统计编程？为什么需要统计编程？?统计人的责任：了解算法的黑箱（我不入黑箱谁入黑箱？）?黑箱：忽悠的源泉（你快回来，我一人忽悠不来！）?统计人的自由：灵活计算与输出（我的地盘我做主！）?我的地盘儿，我也不知道听谁的。?我的地盘儿，听你的。?统计人的愿景：走在时代的前列（直挂云帆济沧海！）?吃罢早饭吃中饭，吃罢中饭吃夜饭，夜饭吃过困觉哉，困觉起来吃早饭。6 2哪里需要编程？哪里需要编程？?没错，很多情况下确实不需要（即便是用R

3、也不需要）?但是，无论如何你的GUI都不可能把所有的方法都囊括到菜单和按钮中的?正态性的检验有多少种方法？Jarque-Bera 检验?SPSS 做一个 t 分布的 KS 检验？（只有正态、泊松、指数等）?前沿：Copulas？Machine Learning？?等等 Jarque-Bera 检验 JB 统计量渐近服从分布 2(2)7 8 一个随意的初级 R 程序：n=length(x)S=mean(x-mean(x)3)/(n-1)/n*var(x)1.5 K=mean(x-mean(x)4)/(n-1)/n*var(x)2 chi=n/6*(S2+(K-3)2/4)pvalue=1-pch

4、isq(chi,2)R 源代码（tseries 包；jarque.bera.test()函数）jarque.bera.test function(x)if(NCOL(x)1)stop(x is not a vector or univariate time series)if(any(is.na(x)stop(NAs in x)DNAME-deparse(substitute(x)9 n-length(x)m1-sum(x)/n m2-sum(x-m1)2)/n m3-sum(x-m1)3)/n m4-sum(x-m1)4)/n b1-(m3/m2(3/2)2 b2-(m4/m22)STATI

5、STIC-n*b1/6+n*(b2-3)2/24 names(STATISTIC)-X-squared PARAMETER-2 names(PARAMETER)-df PVAL-1-pchisq(STATISTIC,df=2)METHOD 1:10#井号是 R 的注释符号 1 1 2 3 4 5 6 7 8 9 10 10:1 1 10 9 8 7 6 5 4 3 2 1 seq(1,10,0.5)#等差数列 1 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 30 11 6.0 6.5 7.0 7.5 8.0 8.5 9.0 9.5 10.0 rep(2,5)

6、#重复相同的对象 1 2 2 2 2 2 rep(1:3,times=3)#观察与下例的不同 1 1 2 3 1 2 3 1 2 3 rep(1:3,each=3)1 1 1 1 2 2 2 3 3 3 rep(1:3,1:3)31 1 1 2 2 3 3 3 产生结构化数据 向量很容易用函数 c()产生：x=c(9,1,1)x 1 9 1 1 (x=c(Xie,Yi,Hui)#为什么用括号 32 1 Xie Yi Hui 矩阵用 matrix()产生：matrix(1:10,2)#注意：默认按列排列 ,1,2,3,4,5 1,1 3 5 7 9 2,2 4 6 8 10 matrix(1:1

7、0,nrow=2,ncol=5,byrow=T),1,2,3,4,5 33 1,1 2 3 4 5 2,6 7 8 9 10 数据框用 data.frame()产生 x=data.frame(1:5,4:8)#把若干个向量合成数据框 x X1.5 X4.8 1 1 4 2 2 5 34 3 3 6 4 4 7 5 5 8 x=cbind(x,c(A,B,C,D,E)#绑上一列字符 x X1.5 X4.8 c(A,B,C,D,E)1 1 4 A 2 2 5 B 35 3 3 6 C 4 4 7 D 5 5 8 E dimnames(x)#看一下 x 的行列名 1 1 1 2 3 4 5 2 36

8、 1 X1.5 2 X4.8 3 c(A,B,C,D,E)colnames(x)#只看列名 1 X1.5 2 X4.8 3 c(A,B,C,D,E)37 colnames(x)=c(X1,X2,X3)#改列名 x#这样看起来就舒服多了 X1 X2 X3 1 1 4 A 2 2 5 B 3 3 6 C 4 4 7 D 5 5 8 E 38 因子用 factor()产生 列表用 list()产生 时间序列用 ts()产生 数据还可以从外部文件读入，甚至是剪贴板（clipboard），或者 SQL、Access 数据库（RODBC 包）例：D:x.txt 文件 V1 V2 V3 V4 V5 39 1

9、 5 9 13 17 2 6 10 14 18 3 7 11 15 19 4 8 12 16 20 x=read.table(d:x.txt,header=T)x V1 V2 V3 V4 V5 1 1 5 9 13 17 40 2 2 6 10 14 18 3 3 7 11 15 19 4 4 8 12 16 20 运算运算 算术运算：+,-,*,/,%（余数）,%/%（整数商）,（乘方）5%2 41 1 1 5%/%2 1 2 25 1 32 逻辑运算：&,|,!（且、或、非）；,=,=,=（“小于一个负数”如何表示？注意-是赋值符号，x-9 与 42 x=9 等价！正确写法是加上空格 x-

10、9 或者 x T&F 1 FALSE T|F 1 TRUE !T 1 FALSE 43 1=T 1 TRUE 2=T 1 FALSE 0=F 1 TRUE 1:10 5 1 FALSE FALSE FALSE FALSE FALSE TRUE TRUE TRUE TRUE TRUE 44 下标的使用（获取元素）下标的使用（获取元素）向量、因子、时间序列 xi；矩阵 xi,j xi,x,j，数据框同矩阵；数组就是根据维度多打几个逗号而已 xi,j,k,；列表要用双重中括号 xi 特殊的使用：在中括号中用逻辑表达式提取元素或给特定元素赋值 xx 5#提出大于 5 的元素 1 6 7 8 9 10

11、45 xx 5&x x=matrix(1:20,4,5)xx=2&x x ,1,2,3,4,5 1,1 5 9 13 17 46 2,2 6 10 14 18 3,3 7 11 15 19 4,4 8 12 16 20 xx=2&x x ,1,2,3,4,5 1,1 NA NA NA 17 2,NA NA NA NA 18 47 3,NA NA NA NA 19 4,NA NA NA 16 20 一些数学和统计函数一些数学和统计函数?最大值max(),最小值min(),均值mean(),标准差sd(),方差var(),相关系数cor(),求和sum(),积prod(),中位数median(),

12、分位数quantile(),对数log(),指数exp(),排列factorial(),组合choose(),四舍五入round(),向下 48 取整 floor(),向上取整 ceiling(),总结 summary(),?累加 cumsum(),秩 rank(),排序 sort(),倒序 rev(),矩阵转置 t(),逆矩阵 solve(),特征根 eigen(),?关于统计分布的四大金刚：pnorm(),dnorm(),qnorm(),rnorm(),(p,d,q,r+分布名称分别构成：分布函数值、密度函数值、分位数、随机数，如 pf()表示 F 分布函数值，runif()表示产生均匀分

13、布的随机数)；抽样 sample(),?线性模型 lm(),广义线性模型 glm(),t 检验 t.test()49 程序控制流：告诉计算机怎么走程序控制流：告诉计算机怎么走?选择控制?if(条件)怎样怎样?if(条件)怎样怎样 else 又怎样?ifelse()函数：ifelse(条件,满足则取我,不满足则取我)?switch()函数：多条件选择（至少我很少用它）50?循环?for(循环变量 in 一个序列中)怎样怎样?while(某条件满足则)怎样怎样 例：cond=T if(cond)winDialog(ok,Hello Kitty!)x=0 for(i in 10:1)51+x=x+i

14、+x 1 55 while(x=1&winDialog(yesno,Im a little pig.)!=YES)+x=x-1+#这个例子比较整人 52 Tips：括号怎么打？小括号一般用于函数，中括号获取元素，大括号把一段完整语句“保护起来”（如果只有一句话，通常大括号可用可不用）至此，我们有了数据结构，知道了怎么控制我们的程序，那么，编程的问题基本已经解决。剩下的工作就是适应 R 的语言环境（这是一个积累的过程，勿急躁）53 国有国法，家有家规。不要对你的计算机想当然。几个小例子：http:/ 不是逻辑运算；0 x|t|)(Intercept)0.009639 0.019638 0.491

15、 0.625 x 0.193976 0.032711 5.930 4.55e-08*-58 Signif.codes:0*0.001*0.01*0.05.0.1 1 Residual standard error:0.0979 on 98 degrees of freedom Multiple R-Squared:0.2641,Adjusted R-squared:0.2566 F-statistic:35.17 on 1 and 98 DF,p-value:4.548e-08 plot(x,y);abline(fit)0.00.20.40.60.81.0-0.2-0.10.00.10.20.

16、3xy 59 60 fit=lm(y x+I(x2)#高次的回归，注意 I()summary(fit)（2）方差分析）方差分析 x=gl(3,50);y=rnorm(150)summary(aov(y x)Df Sum Sq Mean Sq F value Pr(F)x 2 1.743 0.872 0.8053 0.4489 Residuals 147 159.092 1.082 boxplot(y x,notch=T)12 61 3-2-1012 62（3）位置参数的检验）位置参数的检验 参数方法 t 检验：t.test()非参数方法：wilcox.test(),kruskal.test()

17、,等等 条件检验：coin 包（Conditional Inference）（4）广义线性模型）广义线性模型 三个核心概念：?因变量：指数分布族(;,)exp()iiiiii/(,/)iif yA yyA =+63?自变量的线性组合1 1xppx=+?连接函数（link function）：用来建立因变量的分布参数与之间的连接 用得很广泛的所谓 logistic 回归只不过是 GLM 的特例：分布族为二项分布、连接函数取 logit 函数。普通的回归也是GLM 的特例：分布族为正态分布、连接函数为 f(x)=x R 函数：glm(formula,family=gaussian,data,.)6

18、4#Dobson(1990)Page 93:Randomized Controlled Trial:counts-c(18,17,15,20,10,20,25,13,12)outcome-gl(3,1,9)treatment-gl(3,3)print(d.AD|z|)(Intercept)3.045e+00 1.709e-01 17.815 table(rpois(100,5),rbinom(100,1,.5)0 1 0 0 2 1 1 2 2 4 4 3 10 5 72 4 5 10 5 7 8 6 8 7 7 3 7 8 7 0 9 6 0 10 2 1 11 0 1 73 Fisher

19、精确检验：fisher.test()卡方检验：chisq.test()x=matrix(c(12,5,7,7),nc=2)chisq.test(x)Pearsons Chi-squared test with Yates continuity correction data:x X-squared=0.6411,df=1,p-value=0.4233 74（7）混合效应模型）混合效应模型 nlme 包 lme4 包（8）分位数回归）分位数回归 quantreg 包 线性分位数回归函数 rq()75（9）等等等等）等等等等?根据自己的研究需要，可以自己写一点代码，或者用别人的包。?怎么知道该用什

20、么包呢？读书！如：A Handbook of Statistical Analyses Using R Mixed-Effects Models in S and S-Plus 76 4R 图形图形?R 是枯燥的吗？当然不是，你可以用它做很多千奇百怪、匪夷所思的东西看你的创造力如何。?图形由什么构成？点、线、形状、颜色、文本、坐标轴、图例。这一切，在 R 中都有相应的低层函数可以完成。?那么，还有什么作不出来的图？当然没有了。77 本学期我正极其头大地写一本关于 R 作图的书，困难在哪里？在于统计理论。自己回去花两天时间把 graphics 包的所有函数都看一遍，也就到了境界了。演示几个我自己

21、编的小例子：火海刀山？y=.1+runif(20,.2,1)78 xx=c(1,1:20,20)yy=c(0,y,0)plot(xx,yy,type=n,xlab=x,ylab=y)for(i in 255:0)yy=c(0,y-(1-i/255)*min(y),0)polygon(xx,yy,col=rgb(1,i/255,0),border=NA)Sys.sleep(.05)51015200.00.20.40.60.81.0 xy 79 80 时钟？x=seq(0,2*pi,.01)par(pty=s)plot(sin(x),cos(x),type=l)n=60 for(i in 1:n)

22、points(0,0,cex=2)81 arrows(0,0,cos(2*pi/n*(i-1),sin(2*pi/n*(i-1),col=white)arrows(0,0,cos(2*pi/n*i),sin(2*pi/n*i)Sys.sleep(1)arrows(0,0,cos(2*pi/n*(i-1),sin(2*pi/n*(i-1),col=white)lines(sin(x),cos(x)winDialog(,到点了，吃饭去啦，冲啊！nn（旁白：我晕，这个钟走反了吧？-_-/）)-1.0-0.50.00.51.0-1.0-0.50.00.51.0sicos(x)n(x)82 83 置信区

23、间：f=function(n=1000,alpha=0.95,rn=50)d=replicate(n,rnorm(rn)m=colMeans(d)#s=apply(d,2,sd)z=qnorm(1-(1-alpha)/2)y0=m-z*1/sqrt(50)y1=m+z*1/sqrt(50)84 plot(1,xlim=c(0.5,n+0.5),ylim=c(min(y0),max(y1),type=n,xlab=,ylab=)abline(h=0,lty=2)for(i in 1:n)arrows(i,y0i,i,y1i,length=0.05,angle=90,code=3,col=ifel

24、se(0 y0i&0 y1i,blue,red)points(i,mi)85 Sys.sleep(.1)Sys.sleep(3)replicate(2,f(100,.95)020406080100-0.6-0.20.20.6 86 87 重要的不在于图形的花哨，而是表达的信息，还有对于图形原理的深刻理解！（你确信你理解直方图吗？）一句话：在 R 面前，我时常感觉到我渺小得如同蚂蚁的左脚小指的指甲盖儿。88 四、写在最后（给师弟师妹们）四、写在最后（给师弟师妹们）总听师长说要打好基础，我们究竟需要什么基础？我是统计学方向的学生，在统计学院已经呆了近五年，到我这种“老龄”阶段，“打基础”只能幻想了

25、（“我想早恋，但是已经晚了”）那么，想早恋的都趁早吧！89 这是我的个人晚恋总结：这是我的个人晚恋总结：?高等代数和数学分析：矩阵的乘法、矩阵特征根/特征向量（优化问题中常遇到）、迹的性质 tr(AB)=tr(BA)、行 列 式 的 乘 积|AB|=|A|*|B|、正 定（positive definite）与正交（orthogonal）矩阵的定义、各种分 布/密 度 函 数 的 微 积 分、泰 勒 展 开（Talor expansion）、凸函数、偏导数（太重要了）、极值和条件极值（记好拉格朗日乘数法）?概率论：贝叶斯公式、分布的变换（一个随机变量/随机向量的函数的分布）、随机变量的定义、期

26、望和方差的定义（你一辈子都会和它们打交道）、Chebyshev不等式、柯西-施瓦兹不等式、熵和母函数、大数定律和中心极限定理（注意条件！）?数理统计：统计量及其分布（统计量从正态到 t、卡方、F 的一连串构造、经验分布函数）、赵选民那本数理统计书上 P58 公式 2.13：()/(XS/)(0,1)LnnN、90 点估计（矩估计与极大似然估计：发生的都是大概率事件）与区间估计、假设检验（小概率事件不太可能发生）P165 之后的非参数假设检验很重要！?回归：?(1)X XX y=、帽子矩阵 H 等幂且迹为 p+1、?2211niie=1np=、?21()()91 DX X=2a、系数检验的 t统

27、计量构造、F 统计量构造、回归诊断的数学本质、选择子集的三条传统原则：；AIC；RpC、岭回归、非 92 线性回归（迭代求解的方法）、广义线性模型?多元统计分析：玩弄矩阵?非参数统计：玩弄秩，到时候编程序就死命地用rank()函数（有秩走遍天下）?时间序列：平稳的概念、自相关、偏自相关、关于ARIMA：初学的时候盯着一件事AR(1)模型的系数怎样用最小二乘法求解出来，不然学完将一无所获。?抽样：实际上是统计的基础，但是这个地位很少被强 93 调。关于这门课程，思考一个问题就达到境界了：为什么样本均值近似服从正态分布？悟透了这一点，便打通了从概率论到数理统计到抽样的经脉。（我相信有很多人在这里有

28、着根深蒂固的误解，参见某些教科书）94 现在我期待着弄清楚什么：现在我期待着弄清楚什么：因子分析的计算过程以及数学表达式的现实意义、IV（工具变量）与 GMM（广义矩估计）、PLS（偏最小二乘）的现实意义、分位数回归的系数估计方差、对数线性模型、Zero-Inflated 系列模型（ZIP、ZIB 等）、机器学习算法模型的推广性能、SEM（结构方程模型）的收敛性能与估计的稳定性、线性混合效应模型（Linear Mixed Effects Models）的随机效应的现实意义、罚似然、伪似然 95 统计之都的正式建立时间是 2006 年 5 月 19 日，也就是说，过两天就是COS一周年“诞辰”了

29、。一年时间中，我们迅速积累了近 17000 名会员，发表近 3500 篇主题、26000篇帖子（比率说明什么？），每月网站独立访问IP近 30000，每天论坛访问IP数超过 1500。Last but not least,（排名不分先后）96 感谢勤劳的师弟王剑（无痕），通过和他的谈话，我才得以了解到师弟师妹们的学习情况，同时他作为 COS 论坛的另一位管理员付出了辛勤的劳动。感谢热情的 COS 管家杨俏（shashouzhang），她为本次活动以及 COS 论坛的发展都做了大量的工作。感谢李静萍老师、王燕老师、谢邦昌老师、薛薇老师、吴喜之老师，以及 amzon007、netcow、zouwu

30、 等各位老师一直以来学术上的鼎力支持，还有统计学院各位领导老师 97（如王晓军老师、宋大我老师等）长期的关怀！感谢 xifan、colinisstudent、ypchen、abel、jyma、涤尘、cran、nurseshark、shoeda、areg、蟋 蟀、woodcutter、lucky、gracefeng 等各位版主在各自岗位上的认真负责。感谢工作在 COS 网站幕后的几位功臣：鲁明、叶舟、刘秋艳、郝明月、陈喆、雷斯帏！当然，更要感谢广大会员的支持！我不敢列举出一个名单，98 因为对本论坛做出贡献的会员实在太多。若仅从我个人的印象来说，我们要感谢 rtist、micro、liutyy、hcq930、eshanzi、ilikemath、anning189、drewlee、bjt、sociology、anita_jiu、statax、感谢select username from COS 会员数据库 青春年少的人们，让我们再接再厉，用我们的努力打造一个统计学的精品平台，继而推广、普及统计学知识，谢谢！99 IN THE NAME OF STATISTICS,UNITE!WWW.COS.NAME