Entries Tagged '统计备忘录' ↓

更新(技术与工作,2011-10-04)

大半年没跟新了。一些个人技术方面的动向在英文博客里:

1.  CDISC Express

CDISC Express是一款开源的SDTM转化器,基于SAS和Excel(其实背后是XML)。我为这个产品写了一系列的入门教程,并参加了一个相应的比赛(赢了一台iPad)。益辉最近也抽中了一台Kindle,加上上半年在SAS Global Forum拎回来的一台Kindle,不禁感慨,知识(和概率)真是咱农民的第一生产力呀。对着手头这台上世纪的诺基亚6300,我坚定的目光瞄上了iPhone5。有一句话说,做自己喜欢的事,然后找一个人或机构买单,人生一乐也。

又,把临床数据自动导入到符合CDISC标准的数据集,现在有不少vendor在做,很多药厂内部也有开发。感兴趣的朋友不妨也鼓捣鼓捣。

2. Big Data

现在还提这个概念就较不时尚了,重要的是如何做,或者如何为以后做些准备。SAS也将加入对Hadoop的支持,不过现在Hadoop的框架还是基于Java,拣回来些C/C++、Java或者Python还是很有戏场的(甚至R也能够玩会Hapoop),对一个SAS程序员来说,是时候更新自个的工具箱了。

迄今为止,这个领域还是IT人士的天下。这一方面说明了我们这些传统的数据分析人员还需要夯实技术(去看看Hadoop的‘hello world’, 所有的操作都要分解为map()和reduce(),而这并不是所有人都习惯的),也说明目前这一块,发展还在起步阶段,因为这些框架,最终都是要为了让数据产生价值,而这就不可能有数据分析师的缺席。

3. SAS之外

这些日子,在SAS之外,试图拣回来些数据挖掘的东西。本来准备四月份在拉斯维加考下SAS数据挖掘的认证,后来还是被其他事情给挪掉了。最后所得就是,还是稍微看了下这块的东西,在概念上是过了一遍。这次在拉斯维加,碰到以前在北大一起做SAS俱乐部的朋友光辉,他现在加拿大工作。光辉说他趁着开会的缝隙,随手把SAS的数据挖掘认证(SAS Certified Predictive Modeler using SAS Enterprise Miner 5 or 6 Credential)给考了。这叫行动力啊。

机缘巧合,这上半年有朋友给我介绍Feature Selection(大白话就叫变量选择,在数据挖掘里面属于前期的工作),我就对着一个R的包看了一遍,也就是在脑袋里过一下,最后没有什么成果出来。另一个额外的好处是,这次我是对R在学习方面的优势有了切身体验,对R算是有了些个人感情。一边找论文看,然后几乎就能找到一个相应的R包跑一跑,增加些手感,这效果真是不错。

再说些R。我还真尝试过认真学一学R。做为一个SAS程序员,想,还是从input/output开始吧,就先研读官方文档R Data Import/Export,转了一圈,发现,这数据读取转存还是SAS更为方便有效,对我来说,鼓捣这个有点得不偿失,遂罢。

又想,R对SAS的一个优势在于它支持多种数据结构,像array、list之类,就从这入手了。不知道怎么搞的,在学习R的数据结构时,可能是我学习不够深入,修为不够,总觉得有些不够畅快——具体我说不上来。想,我还是把自己的Python拣回来吧,R跟Python有类似的地方。后来我就能说上来了,Python的List等数据结构,对我来说,的确是更为优雅(换一种说法,或者是因为我还没有读到一本“优雅地”介绍R数据结果的文档)。好吧,到此为止,我的R还停留在用它尝试算法的阶段。

然后就是捡回Python。我两年多的SAS程序员工作经历,一个直接的后果就是把之前学过的一点C++、Java和Python又打回“hello world”的水平,手都生疏了。我读了很多Jian Dai(他在加州做SAS程序员)的代码,包括Perl, C, JavaSAcript等等,想,做SAS程序员也不至于把其他的语言都丢掉,为了方便,又重新学习下Python。

为学而学毕竟有局限,在工作中使用当然来得最快。现在我还没有这机会,除了翻文档之外,就尝试用Python写些小东西,把它就当一回事一样。其实,除了工作中直接使用,还有一项就是参加开源运动——我看益辉这几年写R代码是风生水起,一个极大的外部刺激就是参与开源运动,整个积极性都给起来了。一个普通的统计学博士生,用不了那个多的代码量,那动力,就在于统计之外。

统计。今年的基础统计学习,没有去年那么饱满。每天跟统计师打交道,就想好好学习一下他们的语言。今年,到现在为止,算是在工作中(向统计师)学习,整了下CI calculation和euivalence and noninferiority test。下一步准备鼓捣些Sample Size之类。有个哥们写了一篇文章叫Programmers Need To Learn Statistics Or I Will Kill Them All,好吧,作为一个程序员,为了生命安全,每天跟统计师打交道,还是学下统计了。

4. 工作

不是说自己的工作,是说工作与实习机会。感兴趣的朋友,不妨跟我联系。中南大学统计系的本科生韩帅,不久前结束在我公司的实习,在Sxlion主持的SAS中文门户写了两篇非常有价值的总结文章,大伙可以去看看:

一个SAS菜鸟的故事 学习篇

一个SAS菜鸟的故事 实习篇

人是这样的,只有经历过了(我说的是第一份实习或其他),才能自信地自称“菜鸟”,否则就真是论坛里“弱弱地问”的菜鸟。我推荐我公司,当然是因为最熟悉。说实在的,在北京,在药厂,对SAS程序员来说,这的确是一份不错的的实习机会。去年我也推荐了汤耀华来我公司做实习生,他现在中科院读研。这两位都来自中南大学,怎么说?对我和我的同事来说,他们提升了中南大学在我们心目中的分量。在上海,我想,很多SAS程序员对中南大学也很有兴趣。

一些大学生朋友给我写邮件打电话,说些职业或技术选择的事。这不是选择两个到手的offer,大哥,随便扎进一个方向,都比在那犹豫强。

5. SAS

还是得提一下SAS。现在用上SAS9.2了(然后SAS9.3就华丽丽地发布了。。。),对我来说,直接地,有两个触动:

1)SAS9.2及以后,作图更漂亮,我终于有些动力研究下SAS作图了。在公司,去年与宾州的同事一道,把公司作图的macro(基于9.1.3)的更新以后,以后都直接调用,压根就没写过proc gplot之类,然后我就发现自己作图的那一点点手艺正在迅速凋零。。。

2)SAS9.2中,BASE中正式支持自定义函数(仿佛听到其他门类的程序员在冷笑,我弱弱地辩解一下,SAS/IML早就可以自定义函数,而且对SAS程序员来讲,一直使用macro写一种叫function-like macro,而且,macro也支持递归。。。)。这将给我们的编程生活带来不少乐趣,其中之一就是可以把以前的一些macro改写成函数(方便debug啊)。

在写函数之前,我把SAS BASE里面的函数都看了一遍,省得以后重复造轮子:

SAS_Func

 

 

 

 

 

 

 

把所有的函数分门别类导入到Excel,然后动用一点VBA,把每个函数的解释都搁进每个函数名的comments里。这么一梳理,又发现不少有趣的东西,再议。

理解统计学

像外行一样思考,像专家一样实践。

刚在“统计之都”(COS)贴了一篇读书笔记,关于假设检验的基本概念。作为一个非统计出身的SAS程序员,我特别在意如何用大白话来解释统计学的概念:数学公式的推导,我可能走不了多远;但用普通话给表达出来,我其实没什么劣势。写这类文章,我都把自己的认知界定在高中水平,然后一路路给出示例说明,最后达到理解某些概念的目的。比如这篇关于假设检验,基于这段时间来生物统计学基本教材的阅读(这里这里)而做的一个笔记,就是从日常的逻辑推理说开去。从常识说起,能推到多远,就写到多远。

这一年的统计学阅读,多集中在生物统计方面。生物统计基础,是统计学中最保守也是最经典最严格的那部分——这对我来说,一个好处就是,可以借这个机会,细致地补补统计学101。如果有时间,我倒是愿意就每个主题都写一篇。以后如果自己长时间不接触统计学了,通过自己的笔记还可以迅速捡回来(《程序员统计书》?)。时间,写这个东西最耗的就是时间:什么东西都得不厌其烦地写出来。最后自己读一读,都有些琐屑的味道。

———–

在那篇博文里,我借用了一个在线统计学习网站的一个t分布生成器,挺方便的东西:

http://onlinestatbook.com/calculators/

秋冬又是读书季:数学与统计

最近又开始读些基础数学与统计书,看规模,够我读上好几年了。老三样,微积分、线性代数和概率论与统计,大致想让自己的数学水平,随时保持在一个能进研究生院的水平。

工作中一个感受是,如果一个问题成功地转成了数学形式,那你就可以轻松地解出来(随便叫一个学习还算认真的大学生也能解出来)。困难的是如何把问题转换成数学形式。

1. Calculus by Michael Spivak, Publish or Perish Press, 1994

这本书的出版社很有意思,叫“不出版,就完蛋”(Publish or Perish),是作者Spivak拥有的一家小型出版社,而且出的大多是他自个的书。为什么选这本微积分教材?因为IBM的一位工程师Antonio Cangiano在他的博客里有推荐。先推荐一下他的博客

http://math-blog.com

一位软件工程师写的数学博客,内容生动适度。如果打不开网页,你可以直接在google reader里订阅。在一篇名为the most enlightening calculus books的博文里,Cangiano(当然)还提到Apostol以及Hardy的经典教材,但咱嫌它们太难。Spicak的这本书,是introduction to real analysis的水平,觉得对自个是一个挑战,但不妨一跳那种。

2. Introduction to Linear Algebra, 4th edition by Gilbert Strang, Wellesley-Cambridge Press, 2009

Strang这本,可能是世界上最好的线性代数入门书。他今年在MIT开的这门课,可以在MIT开放课程里找到:

http://ocw.mit.edu/courses/mathematics/18-06-linear-algebra-spring-2010/

前些年孟岩写了一些列关于线性代数的文章,现在读起来还是很让人大脑激荡。自个也要去线性空间去巡巡礼了。

又,前些年上海师范大学的数学老师eijuhz,托我在北大图书馆借过这本书。当时想就应该是好书。

3. Introduction to Probability, by Dimitri P. Bertsekas and John N. TsitsiklisAthena Scientific, 2002

这家雅典娜科学出版社(Athena Scientific )又是一家跟本书作者息息相关的小型出版社,这本书,也是MIT数学教授写给工程师的教科书,同样可以在MIT的开放课程里找到。我手头是第一版。

读这本书的缘由,大概是“MIT”这个关键字。无论怎么样,这是一本不错的书,可着劲儿读掉就是。我看各种教材就像我们古时候的字典,没有索引,你只有背下来了才能使用。我常翻阅的资料,很多都是以前用的熟透的教科书。遗憾的是,以前用的教科书,并不是每一部都是经典(所以,现在有些垃圾书也在翻着,路径依赖啊)……

4. Biostatistical Analysis (5th Edition) by Jerrold H. Zar, Prentice Hall, 2009

Zar教授引以为豪的,除了这部最受欢迎的生物统计入门书,还有一首关于spelling check的打油诗。我拿这本书捡回些统计的基础知识。这本书的数学保持在高中水平,这也是我手头要备着上面三个大部头数学书的原因之一。

5. 沈明来编著《生物检定统计法》,台北:九州图书文物有限公司,2007

台版书。一些表达比较有趣,比如提到T分布,说Gosset与1908年“得此式”。

又,内页介绍出版社,说它在台湾大学侧门麦当劳楼上。麦当劳楼上的出版社,你猜我想到哪家?哦,中国人民大学出版社。

生物检定(biological assay)所需之统计法,与一般统计指导书所述法不尽相同。

技术博客重新开张

把以前在space写的文字都导入到这个新博客里了。

这新得白花花扎眼的一年,还想多写些关于SAS程序员本身的文字,关于这个职业,它依托的行业环境等等。SAS程序员在国内还不是一个很兴盛的职业。

还会有关于SAS本身的文字,关于SAS语言,SAS公司,关于它的创始人等等。最近我对SAS的创始人Tony Barr比较感兴趣。

技术本身,这个跟饭碗相关,除了SAS技术,很多笔墨可能会停留在CDISC上面。当然还会有自个兴之所至的其他文字,才年初呢,啥都没定。作为跟“统计之都”的约定,所有跟统计相关的文字,我会首先发布到“统计之都”,然后在自个的博客做个备份:

http://cos.name/author/hujiangtang/

白话统计(4):P-value:一个注脚

/*本文发布在cos.name(统计之都)。*/

*************本书给你数理统计的直观****************************

资料来自美国G.H.维恩堡等著的《数理统计初级教程》(常学将等译,太原:山西人民出版社,1986)

《白话统计(1):平均数、中位数、众数》

《白话统计(2):中心极限定理》

《白话统计(3):决策与风险》

*************************************************************************

郑冰刚提到P值,说P值的定义(着重号是笔者加的,英文是从WikiPedia摘来的):

P值就是当原假设为真时,所得到的样本观察结果更极端的结果出现的概率。

The P-value is the probability of obtaining a result at least as extreme as the one that was actually observed, given that the null hypothesis is true.

以下延续白话系列,解释一下,“什么是P值,什么是极端”,算是郑文的一个长长的注脚。

回到上次的硬币试验,那是一次二项试验,每次试验投100次,记下出现正面的次数,比如,如果

每次出现的正面数都是50,你就有把握认为这是一枚均匀的硬币;

正面数等于45或者等于55,你就有一点点的怀疑它是均匀的;

正面数等于30或者等于70,比较怀疑;

正面数等于10或者等于90,非常怀疑。

如上,正面数和反面数的差异越大,你就越有把握认为硬币不是均匀的(拒绝原假设)。重复一下P值的定义,“P值就是当原假设为真时,所得到的样本观察结果更极端的结果出现的概率”,把这个定义套入上述硬币试验的场景中,比如你观察到“正面数是10或者90,正反面次数差异是80”:

如果原假设为真(硬币是均匀的),P值就是你投100次,所得的正反面数差异大于80的概率。

如果这个P值很大,表明,每次投100次均匀的硬币,经常有正反面差异大于80的情形出现。如果这个P值很小,表明,每次投100次均匀的硬币,你很难看到正反面的差异会超过80。

以前说过,10-90是A博士的接受区域。如果一枚硬币投出的正反面次数,差异大于80,——这真是一个“极端”的情形,连保守的A博士看了都摇摇头,不能接受原假设,只好认为原假设不对,硬币是有偏的。这里的逻辑是:

在假定原假设为真的情况下,出现所看到的偏差(正反面差异为80),是这么地不可能(P值很小),以至于我们不再继续相信原假设。

参考资料:

1. 维恩堡《数理统计初级教程》(常学将等译,太原:山西人民出版社,1986,Statistics: An Intuitive Approach By George H. Weinberg and John Abraham Schumaker)

2. Statistics I: Course Notes, 2008 SAS Institute Inc. Cary, NC, USA

Technorati Tags: ,,,,,,,

白话统计(3):决策与风险

/*本文发布在cos.name(统计之都)。*/

*************本书给你数理统计的直观****************************

资料来自美国G.H.维恩堡等著的《数理统计初级教程》(常学将等译,太原:山西人民出版社,1986)

《白话统计(1):平均数、中位数、众数》

《白话统计(2):中心极限定理》

*************************************************************************

/*读书笔记,白话统计系列,力图用普通话讲述统计学的基本概念。这里的题目是“决策与风险”,讲的就是两类错误(type I and type II errors)。以下改编至维恩堡《数理统计初级教程》(常学将等译,太原:山西人民出版社,1986),英文名叫Statistics: An Intuitive Approach By George H. Weinberg and John Abraham Schumaker。这书几近绝迹,当回文抄公,以期重见天日。*/

1、假设与决策:场景

原假设:硬币是均匀的。   备择假设:硬币是有偏的。

/*当我们难以拒绝原假设时,只能得到结论:原假设也许是真的,现在不能拒绝它。而当我们能够拒绝它时,结论是:它肯定不真。以下的口语表述不如这里明确(和拗口)的,以这里的表述为准。*/

试验:在平坦的地方,独立地投掷硬币100次,每次投掷的结果都做记录。最后,正反面出现的次数分别是:

正面:55  反面:45

提问:根据你所看到的结果,判断一下,你接受还是拒绝”硬币是均匀的“这一假设?

-R博士回答:“拒绝这个假设,因为所得到的正面数超过了反面数的允许界限,这表明硬币是有偏的。”
-A博士回答:“接受硬币是均匀的这一假设。我们不能非难硬币掷出55个正面,45个反面,一个均匀的硬币也能掷出这个比率。”
-R博士:“那什么样的结果才能使你拒绝那假设呢?我的意思是,正面数和反面数应该有多大的差异,才能使你认为硬币是有偏的?“
-A博士:“至少90个正面对10个反面,或者90个反面对10个正面。如果我们的决策是拒绝一个掷出55对45这个比率,或者更高一些比率的硬币,那么这个决策将使我们把许多由于偶然掷出上述比率的均匀硬币都宣判为有偏的。你的看法使得非难一个均匀的硬币太容易了。”
-R博士:“太过分了!至少要掷出90对10的比率你才说硬币是有偏的。你过度的轻信,将几乎不可能拒绝关于硬币是均匀的假设。诚然,你很少拒绝一个均匀的硬币,但对一个有偏的硬币,你也很难拒绝。”

上面的对话应该让大伙体会到了一些假设检验的意思。可以总结一下,对照下面的表格,思路会清晰一些:

判定 \         假设
拒绝 第I类错误α  没有错误1-β
接受 没有错误 第II类错误β

A博士(Accept,接受)的法则是,除非试验得到的比率超过90比10,否则就接受硬币是均匀的这一假设。A博士厌恶犯否定均匀硬币的错误(”弃真“,第I类错误),他的法则使得犯这种错误的概率最小。由于均匀的硬币几乎不会出现超过90比10的比率,他很少冒把一个均匀的硬币说成有偏的风险。然而,他付出的代价是,大大降低了试验的检测能力(power,见下),他的法则使得拒绝假设是极端困难的。大量有偏的硬币也不会出现如90对10这样大的差异,因此它们也会被当成均匀的硬币而没有被检测出来。可以说,A博士对接受假设有偏爱,当假设为真时,他很少犯拒绝它的错误;但当假设不真时,他会常犯接受它的错误。

R博士(Reject,拒绝)的法则是,除非比率低于55对45,否则就不能接受硬币是均匀的这一假设,也即,仅当硬币的正反面数差异在一个狭窄的界限之内,她才接受假设。她把试验看成类似9.11时美国进行的安全检查(”宁可错杀三千,不可错过一个“),重要的是检测出有偏的硬币。R博士的法则在接受错误的假设方面所冒的风险极小(”取伪“,第II类错误),代价是增加了把一个均匀硬币判成有偏的风险。可以说,R博士对拒绝假设有偏爱,当假设碰巧不真时,她很少犯接受它的错误;但当假设碰巧为真时,她常犯拒绝它的错误。

2-1、决策与风险(用均匀的硬币做试验,第I类错误)

一次试验,不足以判断两位博士谁的法则是正确的。现在,用一个均匀的硬币(我们知道,两位博士不知道,这里的原假设是硬币是均匀的),把上面提到的投硬币试验,重复100次(每个试验由100次投掷构成),那么,记录下的正面数X,将构成一个二项分布,X~B(n,p),其中,n=100,p=0.5。根据某个中心极限定理,正态分布是二项分布的极限分布,上面的二项分布可以由均值为np=50,方差为np(1-p)=25的正态分布来近似。又因为二项分布只取整数值,在近似它的正态曲线下会出现很多空隙,为了校正这种情况,可以把整数的两头各扩大0.5个单位,以这个区间表示正态曲线下的那个数。

对R博士来说,仅当掷出的正面数多于45,少于55时,她才接受假设。在正态曲线下,这两个端点可以写成45.5和54.5。

——|-/////-|———
      45.5    54.5

标准化,(45.5-50)/5=-0.9,(54.5-50)/5=0.9,根据标准正态表,可知45.5-54.5这个接受区域包括了总面积的63%。也即,投掷均匀硬币所产生的样本中,有63%的样本,其正面数落在接受区域,相应地,其正面数落在R博士提出的否定域的概率为37%。也就是说,当硬币是均匀的时,R博士犯第I类错误的概率为37%。对A博士来说,他的接受区域在10-90之间,他几乎不会犯第I类错误。

2-2、决策与风险(用有偏的硬币做试验,第II类错误,功效)

现在取一个有偏的硬币(我们知道,两位博士不知道,这里的原假设还是硬币是均匀的),即投出正面的概率不等于二分之一(注意,说硬币是有偏的,并不必对p的值作出指定,因为硬币有偏可以有无限多种方式)。为了评价两位博士的法则在拒绝假设方面有多大的成功,我们需要对硬币指定一个偏度,比如是掷出正面的概率是0.6,做上面同样的100次试验(每次试验有100次投掷),近似成一个正态分布,均值np=60,方差是np(1-p)=24。

对A博士来说,他的判定法则是,只要得到的正面数在10到90之间就接受假设。显然,即使一个有偏的硬币所得到的正面数,也位于A博士的接受区域里。即,当硬币出现正面的概率为0.6时,A博士还是经常

白话统计(2):中心极限定理

*************本书给你数理统计的直观****************************

资料来自美国G.H.维恩堡等著的《数理统计初级教程》(常学将等译,太原:山西人民出版社,1986)

《白话统计(1):平均数、中位数、众数》

*************************************************************************

定理1:(中心极限定理)假定大量的等容量随机样本都是从同一无限总体采样的,算出每个样本的和,并把不同样本的和放在一起以形成一个新的分布,于是这个新的分布就是渐近正态的(其中要假定产生这些和的随机样本每个容量是足够大)。

在传统概率论教科书上,一般会这么陈述这个定理:

定理1`:(独立同分布的中心极限定理)设随机变量X1,…Xn,…相互独立,服从同一分布,且具有相同的数学期望和方差,则随机变量之和ΣXi的标准化变量服从标准正态分布。

演示性例子

想像一个很大的箱子,装满了小纸条,可供我们无穷无尽地抽取,每张纸条上写有一个数字。为简单起见,假定只有0、1、2三个数字,且每个数字出现在每张纸条上的可能性都是1/3。记住,这个箱子里的纸条如此之多,以致我们可以抽取任一数目的任一种纸条,而不必担心会改变箱中剩下的各种纸条之间的比例。

箱子有一个小口,通过它,每次可以释放出一张纸条。箱子还有一个洗牌装置,这种装置会把纸条洗得这样得均匀,以至当我们决定抽取一张时,每张纸条有同样的被释放出来的机会。因此,我们的观察室独立的,而且我们的样本是随机的。

现在我们就来抽取等容量的随机样本,假设每个样本都包含200张纸条。

我们一张一张地抽取200张纸条。比如头一张纸条上的数字是2,第二张纸条的数字是0,第三张纸条是2,如此等等。假设构成这个第一份样本的200张纸条上的数字总和是210,这个和成为所产生的新的分布的第一项。

第二个样本的200张纸条上的数字之和比如是194.对大量的样本,每个样本都包含200张纸条,重复这个过程。定理1告诉我们,这种样本和数越来越多时,样本和的分布近似于正态分布。

如何实际运用定理1

关于定理1,对被抽取样本的那个总体没有要求任何限制。不管被抽取样本的那个总体,其分布的形状如何,样本和的分布都是正态的。

定理1说明,为什么正态分布出现在如此多的不同的问题之中。我们用于纸条取样的那种方法,看来是实际中特别喜欢使用的一种方法。在每次情况中出现的、构成一个正态分布的那些数,都可以看作独立观察资料的等容量样本的和

例子1。考察射击时围绕靶子构成正态分布的子弹。每一颗子弹击中的位置实际上是许多随机影响的和,比如姿势、风向、光线、心理等等。这些因素和诸如此类因素的影响,同时在一位特定射手的身上起作用;且对于不同的射手,它们是不同的。一个射手的得分,表明他的子弹最终射到何处去了,这个得分是那些随机影响的样本之和。具体地,比如每一个射手的分布式70项主要影响之和,因而每一发子弹的得分,都可以看作是70项的一个样本和(与70张纸条上的那些数字的和相对应)。这样一来,不同射手的得分,就可以看作是不同的等容量样本的和。根据定理1,子弹得分的分布式正态的。

例子2。考察每个人的智力水平,也可以当作出自不同根源的小影响的和来看待,包括营养、机会、性格、遗传等等等等。这么看来,大量的人的智力水平的分布式正态的。

定理2:(定义1的一个变形,平均数的中心极限定理)假定,大量的等容量随机样本是从同一无限总体中采集的,算出每一个样本的平均数,并把不同样本的平均数放到一起形成一个新的分布,于是这个新的分布就是渐近分布的(假定产生这些平均数的随机样本容量是足够大的)。

样本平均数的集合可以通过样本和集合直接得到,因此平均数的分布就是和的分布的一个小比例的变形。样本平均数的分布用两个有用的性质:

  1. 假定无穷多个等容量随机样本是从同一无限总体中抽取的,而且把这些样本的平均数放到一起,以构成一个新的分布,那么这个新分布(样本平均数构成的)的均值与原总体的均值相同。
  2. 假定无穷多个等容量随机样本是从同一无限总体中抽取的,以n表示每一个样本的容量,这些样本的平均数的分布有一个标准差,它等于原总体的标准差除以n的平方根。

定理2及其两个性质就是我们熟悉的mean(X)~N(μ,σ**2/n)。

白话统计(1):平均数、中位数、众数

********************示例数列为4、4、5、7、10********************

定义1. 各项累加之和除以项数,所得之数值,叫做平均数(6)。

定义2. 众数是出现次数最多的那一项的数值(4)。

定义3. 中位数是这样一个项的数值(或两项之平均数):它(或该两项)的数值大于或等于其余一半项的数值,而小于或等于另外一半项的数值(5)。

注1:平均数受每个改变的牵动,中位数和众数却只受某些改变的牵动。由于这个原因,平均数常常被认为是“敏感的”,是“反映整个分布的”。

注2:当分布含有少数在一端远离的极端项时,平均数的敏感性,对它的代表性反而可能是不利的。

注3:中位数不受少数极端值的影响。

注4:众数显示最大频率,而最大频率是最普遍化的同义语。

一、平均数的两个性质

物理模型:设想一块实验板,其上有4至10的等距刻度。假定5个一磅重的砝码置于板上,砝码的位置为4、4、5、7、10。现在假设有一支点,令该实验板连同其上的砝码在支点上保持平衡,并设该实验板本身无重量。几经失败之后,我们设想找到了平衡点。

 

image

 

 

 

 

结果,支点在各项的平均数之下。任何项的集合的平均数都是那些项的平衡点,而且平均数是唯一的平衡点。

考虑支点左边每个砝码与支点(平均数)的距离。最左端的两个砝码与平均数相距2个单位,第三个相距1个单位,总和是5个单位。同样,支点右边的砝码与支点的距离也是5个单位。正是这个等量平衡了实验板。左右两边砝码与平均数的距离之和总是相同的,也即平均数两边的距离“相抵”。

为了得到上面的距离,从各项减去平均数。有负号的是平均数左边砝码得到的结果,有正号的是右边砝码得到的结果。前面提到,把负号抹去,左边各项的和等于右边各项的和(距离“相抵”)。留住负号,则所有项的和等于0 ,Σ[x-mean(x)]=0,即,与平均数之差的和恒等于0。总结一下:

定理1:在任何分布中,各项对平均数的差之和等于0。反之,若分布中各项对某数的差数之和等于0时,该数即为此分布的平均数。

又考虑对平均数的差数平方:

各项 平均数6做参考点 差数 差数平方
4 6 -2 4
4 6 -2 4
5 6 -1 1
7 6 1 1
10 6 4 16
    总和 26

备忘录之主成分分析

//这个暑期小学期,跟管理技术系的同学上营销数据分析。这门选修,主要是为凑学分了,偶尔也在课程bbs发些贴,权当备忘录了。这个主成分分析,主要来自高惠璇老师的两本教材,《应用多元统计分析》(北京大学出版社,2005)和《实用统计方法与SAS系统》(北京大学出版社,2001)。这个贴就发在学院bbs上,其他ip登不上,贴过来。在记事本上写的,像Z1这样的标记就是Z加一个下标1。

在实际问题中,如果变量个数太多:
1.变量彼此存在一定的相关性,因而观测到的数据在一定程度上反映的信息有所重叠(多重共线性);
2.在高维空间中研究样本的分布规律比较复杂,变量太多增加分析问题的复杂性。

我们希望:
1.用较少的综合变量来代替原来较多的变量(降维,减少维度);
2.这几个综合变量能够尽可能多地反映原来变量的信息;
3.这几个综合变量彼此互不相关。

主成分分析、因子分析、典型相关分析和偏最小二乘回归等统计方法就是以上降维思想的体现。

一、主成分分析

1.定义

这里说的主成分分析,又称主分量分析、主轴分析,是将多个指标(变量)化为少数几个综合指标的一种统计方法。
把p个变量X1,X2,…,Xp,记为一个p维的随机向量X=(X1,X2,…,Xp),其协方差阵为D(X)。考虑X的线性变换:

Z1=A1*X
Z2=A2*X
……
Zp=Ap*X,

这里的X和A1、A2、…、Ap等都不妨看成向量形式。假如我们想用Z1来代替原来的p个变量,这就要求Z1尽可能多地反映原来p个变量的信息。这里“信息”可以用Z1的方差Var(Z1)来表示,方差Var(Z1)越大,表示Z1包含的信息越多。当然,这需要强加一些数学上的限制,否则Var(Z1)就可能是无限大了,这里的限制是向量A1和它自己的转置之积等于1,记为A1*Trans(A1)=1。就这样:

         若存在满足A1*Trans(A1)=1的A1,使得Var(Z1)最大,则称Z1为为第一主成分,或第一主分量,Z1=A1*X。

如果第一主成分不足以代表原来p个变量的绝大部分信息,我们就可以考虑X的第二个线性组合Z2=A2*X。此时,我们要求,已经体现在第一主成分Z1中的信息不要出现在Z2中,即Z1和Z2的协方差Cov(Z1,Z2)=0。就这样:

          在Cov(Z1,Z2)=0时,若存在满足A2*Trans(A2)=1的A2,使得Var(Z2)最大,则称Z2为为第二主成分,或第二主分量。

类似我们可以定义X的第三主成分,以致第p主成分(当然,对p维的随机向量X来说,第p主成分就没有必要了)。

2.直观解释

从代数上讲,主成分就是p个原始变量的一些特殊的线性组合。
从几何上讲,这些线性组合是把由X1,X2,…,Xp构成的坐标系通过旋转而产生的新坐标系。

3.求解

从上面的定义可知,求第一主成分Z1=A1*X的问题,就是求A1,使得在A1*Trans(A1)=1下,Var(Z1)达到最大,这是一个条件极值问题,可以用拉格朗日乘子法求解。

   Max Var(Z1)=Var(A1*X)
   s.t A1*Trans(A1)=1

运用拉格朗日乘子法,令L=Var(A1*X)-t(A1*Trans(A1)-1),接下来就是一阶条件了,按下不表。
说说结果。求解以上条件极值问题,就是求X的协方差阵D(X)的特征值(就是拉格朗日乘子t)和特征向量(就是我们感兴趣的A1、A2、…、Ap)。一个定理:

设X=(X1,X2,…,Xp)是p维的随机向量,其协方差阵为D(X),D(X)的特征值为t1>t2>…>tp,A1、A2、…、Ap为相应的单位正交特征向量,则X的第i主成分为Zi=Ai*X.

4.一些名词

原总体X的总惯量(总方差):就是上面我们求出来的特征值之和,sum(t1,t2,…,tp)
贡献率:单个特征值ti与总方差之比,称为主成分Zi的贡献率。同样可以定义累计的贡献率。
因子负荷量(因子载荷量):主成分Zi与原始变量Xi的相关系数。

5.计算机求解(SAS)

(在SAS系统中,主成分分析的程序步是princomp,因子分析的程序步是factor。所以可以知道主成分分析和因子分析毕竟不是一样的东西,具体见因子分析的部分。)

这里有一个数据,5个变量,id、x1身高、x2体重、x3胸围、x4坐高,30个观测值。我们想对x1身高、x2体重、x3胸围、x4坐高做主成分分析。

data princomp;
  input id x1 x2 x3 x4  ;
  cards;
1 148 41 72 78  
2 139 34 71 76
3 160 49 77 86
4 149 36 67 79
5 159 45 80 86
6 142 31 66 76
7 153 43 76 83  
8 150 43 77 79
9 151 42 77 80
10 139 31 68 74
11 140 29 64 74
12 161 47 78 84
13 158 49 78 83
14 140 33 67 77
15 137 31 66 73
16 152 35 73 79
17 149 47 82 79
18 145 35 70 77
19 160 47 74 87
20 156 44 78 85
21 151 42 73 82
22 147 38 73 78
23 157 39 68 80
24 147 30 65 75
25 157 48 80 88
26 151 36 74 80
27 144 36 68 76
28 141 30 67 76
29 139 32 68 73
30 148 38 70 78
;
proc princomp;
  var x1 x2 x3 x4  ;
run;

结果解释:

(1)结果首先有相关矩阵,可以看出这四个变量之间有较强的相关性。
                                          Correlation Matrix
                                      x1          x2          x3          x4
                          x1      1.0000      0.8632      0.7321      0.9205
                          x2      0.8632      1.0000      0.8965      0.8827
                          x3      0.7321      0.8965      1.0000      0.7829
                          x4      0.9205      0.8827      0.7829      1.0000

(2)然后是以上相关阵的特征值t(Eigenvalue,上面例子中的t)
                                         Eigenvalues of the Correlation Matrix
                                   Eigenvalue    Difference    Proportion    Cumulative
                        1    3.54109800    3.22771484        0.8853        0.8853
                        2    0.31338316    0.23397420        0.0783        0.9636
                        3    0.07940895    0.01329906        0.0199        0.9835
                        4    0.06610989                             0.0165        1.0000

从中可以看到,第一主成分的贡献率就达到 0.8853。单个贡献率就是某个特征值与总特征值的比,这里0.8853=3.54109800/(3.54109800 +0.31338316 +0.07940895+0.06610989) .还有,前两个主成分的累计是0.9636,所以用这两个主成分就可以很好地概括这些数据。

(3)还有特征向量( Eigenvectors,我们例子中的A)
                                                    Eigenvectors
                                      Prin1         Prin2         Prin3         Prin4
                      x1      0.496966      -.543213      -.449627      0.505747
                      x2      
0.514571      0.210246      -.462330      -.690844
                      x3      0.480901      0.724621      0.175177      0.461488
                      x4      0.506928      -.368294      0.743908      -.232343

我们可以由最大的两个特征值对应的特征向量,写出第一和第二主成分(结果中的Prin1和Prin2):
Z1=Prin1= 0.496966*X1+0.514571 *X2+0.480901*X3+0.506928*X4
Z2=Prin2= -0.543213 *X1+0.210246 *X2+0.724621*X3 -0.368294*X4

这样我们就写出了开头那个线性组合。接下来说些因子分析,因子分析是主成分分析的推广,稍复杂些。(待续)

Technorati Tags: , , ,