标签归档:概率

需要相亲几次才能找到靠谱的对象?

谈到相亲就不得不提到著名的麦穗问题。说有一天,苏格拉底带领几个弟子来到一块成熟的麦地边。他对弟子们说:“你们去麦地里摘一个最大的麦穗,但要求只能摘一次,只许进不许退,我在麦地的尽头等你们。”可以看得出,相亲这种活动就有点类似于摘麦穗,在等待和决断之间达成平衡是解决问题的重点。

将上述的麦穗问题进一步抽象就是一个经典的概率问题。若一个袋子里有100个不同的球。每个球上标明了其尺寸大小。我们每次随机无放回的从袋中取一个球出来,观察其大小属性之后需决定要或是不要。如果要,取球就此停止。如果不要, 再继续取球,但不准再回头要原先的球。这样下去,直到100个球取完为止。目的就是取到那个最大的球。

对于这个概率问题,一种思路就是取1到100之间的某个数字n,以它作为分割点将整袋球划分为两组,第一组即从第1个到第n个球,第二组即从第n+1个到第100个球。我们以第一组为观察对象,找到第一组中最大的球M,记录其大小但并不行动。然后从第二组中寻找大于M的第一个球,取该球为最终选择。那么n应该设为多少,才能使取到最大球的概率尽可能的高呢?Ross在其《概率论基础教程》中已经给出了精确的解析(英文版第8版P345;或者陈木法的《随机过程导论》 P105),最优的n公式表达为 $n^*=\inf\{r:\sum_{n=r}^{N-1}\frac{1}{n}\le1\}$;在N充分大时,n应该取在1/e的比例处,也就是所有球数目的37%处。在解的过程中先运用条件概率得到全概率表达式,再用连续化的积分来近似离散化的概率,将积分求出后进行求导得到最终答案。看到这里各位是否有些惊讶,这个e可是无所不在啊。

继续阅读需要相亲几次才能找到靠谱的对象?

真理在缩水,还是上帝在掷骰子?

最近在Google Reader中看见科学松鼠会有两篇文章被频繁分享,名为《真理在缩水——现代科学研究方法并不尽善尽美?》()与(),下文简称《缩水》。文章很有意思,而实际上说的是我们的老本行——统计学,因此我在这里也发表一些我的想法和理解,包括这两年我在美帝学习的一些思考,部分内容受益于两位老师Kaiser和Nettleton教授,先向他们致谢(尽管他们永远都不会看到这篇文章)。同时我也要先说明一下,读这篇文章可能会很花时间(至少我花了大约二十小时写这篇文章),即使我的观点没有价值,我相信里面的引用文献是有价值的。

初读文章,我脑子里冒出的一句话是“上帝在跟我们掷骰子”,文中给出了大量的不可重复的试验,仿佛就像那些号称“具有统计学意义”(下文我再说这个所谓的“意义”)的试验结果只是若干次骰子中的一次运气好的结果而已。读完文章,我们可能不禁要问,到底是真理在缩水,还是它根本就不曾存在?下面我从四个方面来展开,分别说明人对随机性的认识、统计推断的基石、让无数英雄折腰的P值、以及可重复的统计研究。

一、感知随机

随机变量在统计分析中占据中心地位,数学上关于随机变量的定义只是一个“干巴巴的函数”,从样本空间映射到实数集,保证从实数集上的Borel域逆回去的集合仍然在原来的sigma域中即可。随机变量的性质由其分布函数刻画。写这段话的目的不是为了吓唬你,也不是为了作八股文,而是来说明我为什么不喜欢数学的理由,对我而言,我觉得有些数学工具只是为了让自己不要太心虚,相信某时某刻某个角落有个理论在支撑你,但后果是弱化了人的感知,当然,也有很多数学工具有很强的直觉性(如果可能,我想在未来下一篇文章里面总结这些问题)。我一直认为很多人对随机性的感知是有偏差的,对概率的解释也容易掉进陷阱(参见Casella & Berger的Statistical Inference第一章,例如条件概率的三囚徒问题)。

《缩水》一文发表了很多不可重复的试验案例,我们应该吃惊吗?我的回答是,未必。举两个简单的例子:

第一个例子:很多数据分析人员都很在意所谓的“离群点”,论坛上也隔三差五有人问到这样的问题(如何判断和处理离群点),而且也有很多人的做法就是粗暴地删掉,我从来都反对这种做法。除了基于“数据是宝贵的”这样简单的想法之外,另一点原因是,离群点也许并非“异类”。离群点是否真的不容易出现?请打开R或其它统计软件,生成30个标准正态分布N(0, 1)随机数看看结果,比如R中输入rnorm(30),这是我运行一次的结果:

> rnorm(30)
 [1]  1.19062761 -0.85917341  2.90110515  0.59532402 -0.05081508 -0.06814796
 [7]  2.08899701  0.76423007  0.92587075 -1.16232929 -0.68074378 -1.40437532
[13] -0.17932604 -0.72980545 -0.53850923  0.21685537 -0.35650714 -1.32591808
[19] -0.88071526 -1.25832441  0.24001498 -0.41682799 -0.09576492 -0.17059052
[25] -0.99947485  0.25108253 -0.47566842 -0.28028786  0.79856649 -0.13250974

30在现实中是一个比较小的样本量,你看到了什么?2.901?它接近3倍标准差的位置了。还有2.089?……如果你不知道这批数据真的是从标准正态分布中生成出来的,现在你会有什么反应?把2.9删掉?标准正态分布是一个在我们眼中很“正常”的分布,而一个不太大的样本量一次试验足以生成几个“离群点”,那么要是成千上万的试验中没能产生几项震惊世界的结果,你会怎样想?(上帝的骰子坏掉了)

另一个例子和统计学结合紧密一点,我们谈谈残差的QQ图。QQ图是用来检查数据正态性的一种统计图形,与腾讯无关,细节此处略去,大意是图中的点若呈直线状(大致分布在对角线上),那么可以说明数据的正态性比较好,因此QQ图经常被用在对回归模型残差的正态性诊断上。我的问题是,即使数据真的是正态分布,你是否真的会看见一些分布在直线上的点?若答案是否定的,那么我们就得重新审视我们对分布和随机的认识了。下图是一幅教科书式的QQ图(仍然基于30个正态分布随机数):

“正常的”QQ图
“正常的”QQ图(来自R语言qqnorm(rnorm(30)))

随机性并没有这么美好,即使数据真的来自正态分布,你也有可能很容易观察到歪歪扭扭的QQ图,尤其是小样本的情况下。比如下图是50次重复抽样的正态数据QQ图,它和你想象的QQ图本来的样子差多远?

library(animation)
set.seed(710)
ani.options(interval = 0.1, nmax = 50)
par(mar = c(3, 3, 2, 0.5), mgp = c(1.5, 0.5, 0), tcl = -0.3)
sim.qqnorm(n = 30, pch = 19, col = "red", last.plot = expression(abline(0, 1)))
真实的正态分布QQ图
真实的正态分布QQ图(图中直线为y = x)

正态分布是统计学中比较“正常”的一类分布(台湾学者甚至译为“常态分布”),我们尚不能很好感知它的随机性,就不必说其它“怪异”的分布了。

这是我想说的第一点,作为人类,我们对上帝的骰子至少在感知上就可能有问题(别误会,我不信教),接下来从理性角度分析一下。

二、统计推断

《缩水》一文中提到的基本上都是统计推断方法带来的结果,为了理解这些结果,我们必须三思统计推断本身的过程。一般说来,统计推断有两种途径:随机试验和(概率)统计模型,或者说基于试验的推断和基于模型的推断。前者的代表性方法为置换检验(Permutation test),不过它似乎被大多数人遗忘了,更多的人拿到数据首先想的就是用哪个统计模型和分布;实际上,置换检验是一种极具代表性的统计推理方法,可以用典型的“三段论”来说明它(参见去年江堂的文章):

  1. 要么A,要么B
  2. 若有A,则有C
  3. 若非C,则非A,于是B

置换检验的场景是,一次试验中,我们为试验单元随机分配不同的处理(treatment),为了简单起见,假设这里只有两种处理水平A和B,我们想知道两种处理在试验单元的某项指标上是否有显著差异。逻辑是这样:假设处理毫无效果,那么某一个试验对象的指标将不受处理影响,不管我们给老鼠嗑的是A药还是B药,结果都一样,那么我们可以把处理的标签随机打乱(某些A、B随机互换),打乱后A组和B组的差异不应该会和原来的差异很不一样(因为药不管用嘛),否则,我们恐怕得说,药还是管用的。就这样,我们随机打乱标签很多次,形成一个“人工生成”的AB差异分布(因为我们生成了很多个差异值),看原来的差异在这个分布的什么位置上,如果在很靠近尾巴的位置上,那么就认为P值很小。当了个当,当了个当,P值出场了。对置换检验熟悉的人是否有想过,好像我们一直没谈到什么分布的假设,那这个概率值(P值)是从哪里生出来的?答案是最初的“随机分配处理到试验单元上”。这就涉及到试验设计的一大原则:随机化。为什么要随机化?因为试验单元之间可能本来就有差异,换句话说,如果你不随机化,那么你观察到的差异有可能来自试验单元本身。比如,你从笼子里抓前10只老鼠给嗑A药,后10只老鼠给B药,这就没有随机化,前10只老鼠可能很笨或是老弱病残,容易被你抓住,而后10只老鼠跑得贼快。随机化将这些个体差异转变为了随机的误差,例如两只老鼠之间的确本身有差异,但我要是把它们随机分配给处理,那么这种个体差异就会随机进入处理组,保证统计推断有根基。如果这一点没有理解透彻,试验人员很容易在数据收集阶段就已经收集了错误的数据。《缩水》一文中的试验都是怎么做的,我没空去细究。

基于模型的推断的一大特征就是开始对随机变量做一些分布上的假设,例如两样本t检验,假设样本来自独立同方差的正态分布。仔细想想这里面的问题,对建模和理解模型结果有至关重要的作用。一个最直接的问题就是,假设条件是否可靠?George Box大人很久很久以前就说了一句被引用了无数遍的话:所有的模型都是错的,但有些是有用的。统计学方法很“滑”(用麦兜的话说),它的科学与艺术双重身份,总让人感觉拿捏不准。学数学的人可能觉得艺术太不靠谱,其它外行人士可能觉得科学太神秘。这个问题我不想作答,也无法作答,搁在一边,先说一些我曾经考虑过的问题,它们与《缩水》一文可能并没有直接联系,但应该能或多或少启发我们从源头考虑统计模型,直接上手统计模型的隐患就在于你认为一切都理所当然,随着时间推移,假设渐渐变成了“公理”和“常识”,我觉得这是很可怕的。

第一个问题是似然函数(likelihood function),它是频率学派的命脉,而我们大多数人仍然都是频率学派的“教徒”。对于离散变量来说,基于似然函数的方法如极大似然估计很容易理解:我们要找一个参数值,使得观察到的数据发生的概率最大。这里的“概率”二字应该被重重划上记号!接下来我要提一个让我曾经觉得后背发凉的问题:

为什么对连续变量来说,似然函数是密度函数的乘积?

你是否想过这个问题?我们知道连续变量取任何一个值的概率都是0,也就是说,如果按照概率的方式解释似然函数,那么连续变量的似然函数应该是0才对,换句话说,你观察到的数据发生的概率是0。现在你觉得似然函数还是一个理所当然的统计工具吗?

有一位统计学家叫J. K. Lindsey,1998年在(英国)皇家统计学会宣读了一篇论文,叫Some statistical heresies(一些统计异端邪说),如果你没见过统计学家打仗,那么这篇论文将会让你看到一场超大规模的战争,后面的讨论者中包括Murray Aitkin、D. R. Cox和J. A. Nelder等老江湖。Lindsey的文章几乎是被大家围攻了(期待我这篇文章也能被围攻),不过他对于似然函数的解释倒是让我有点茅塞顿开。细节我不展开,大意是,似然函数也是一种近似(用积分中值定理去想)。

第二个问题是渐近统计(asymptotic statistics),同样,这也是统计学家的日常工具之一,因为太常见,所以也有一种理所当然的味道。比如我们看到列联表就想到卡方检验(检验行列变量的独立性),殊不知卡方检验只是一种大样本下的近似方法。渐近统计的基石是什么?如果你要挖,我想挖到头你一定会挖到泰勒展开。至少目前我认为渐近统计“基本上就是泰勒展开的第二项(或少数情况下第三项)”。理论上泰勒展开有无穷多项,我们往往根据一些假设条件,把后面的高阶项都消灭掉,剩下一次项或二次项。比如你展开似然函数,就得到了似然比检验的卡方分布;你展开极大似然估计的一个连续函数,你就有了Delta方法(当然,需要依分布收敛的前提);就这样左展右展,展出了中心极限定理(对特征函数或母函数展开),展出了拉普拉斯近似(对对数密度函数展开)。如果你能看到这一点,就不会奇怪为什么正态分布在统计推断中有如此中心地位(谁叫正态分布的对数密度函数是个二次函数呢)。

第三个问题是,贝叶斯方法又如何?既然频率学派中几乎处处是近似,贝叶斯学派是否会好一些?我的回答是好不到哪儿去。贝叶斯的逻辑很简单,但由于灵活性太强,应用时非常摧残人的脑力,导致争议不断(先验分布怎么取、MCMC是否收敛等)。在《缩水》一文中,恐怕是没有基于贝叶斯方法的试验,因为在所谓的科学试验中,人们往往排斥“先验”这种想法,就好像先验一定代表主观、而客观一定正确一样。逻辑上,这是荒谬的。关于这些问题的重磅讨论,可参考Efron去年发表的The Future of Indirect Evidence以及文章后面Gelman他们三个人的讨论。我总感觉这不是我这个年龄应该看的东西,太哲学了。

我提到这些问题,本意是给统计学家看的,如果你是一个合格的统计学家,你自己首先应该意识到统计学的局限性,而不是拿到数据就分析。

三、万能的P值?

早些年当我还是个无知轻狂小子的时候,我曾戏称“统计软件就是为了放个P”,这里的P指的是P值,不是粗话。这话好像也不全然轻狂无知。使用统计方法的人,难道不是在追逐一个小于0.05的P值吗?如果你的结果不显著,那么肯定发表不了。换句话说,发表的结果很有可能是我们在自欺欺人。下面的漫画生动刻画了人们寻找P值的过程(来自xkcd):

Significant
Significant

重大科学发现!吃绿色的软糖会让你长痘痘!95%置信度!

当你看到95%的时候,你看不到红色的、灰色的、褐色的、橙色的、黄色的软糖……这便是《缩水》一文中说的“发表偏见”(publication bias,“偏见”翻译似乎不妥),即发表的结果是经过人工选择的,你已经不能用正常的概率意义去解读它,或者说,概率已经变了样。

插一句“统计学意义”:这个概念本来的英文是statistical significance,但是被很多专业的人翻译为“统计学意义”,我一直认为这很不妥,给人一种错觉,仿佛统计学保证了什么东西有意义一样,我提倡的译法是“统计显著性”。尤其是“由于P值小于0.05,所以具有统计学意义”这种话,我觉得应该见一句删一句。

上面的软糖问题涉及到传统的多重比较(multiple comparison)与P值调整,这里“传统”指的是要控制族错误率(Familywise Error Rate,下文简称FWER)。所谓控制FWER,也就是要使得一族(多个)统计检验中,一个第一类错误都不犯的概率控制在一定水平之下,比如0.05。让多个检验全都不犯错和单个检验不犯错(指第一类错误)显然是有区别的,比如假设所有的检验都是独立的,一个检验不犯错的概率是95%,两个都不犯错的概率就变成了95% * 95% = 90.25%,检验越多,不犯错的概率就越小。把整体的第一类错误率控制在某个alpha值之下,就意味着单个检验必须更“严格”,因此我们不能再以0.05去衡量每个检验是否显著,而要以更小的值去衡量,这就是P值的调整,老办法有Bonferroni等方法。

控制FWER显得有些苛刻,比如有10000个基因都需要检验在不同处理下的表达差异,那么要是按照传统的P值调整方法,恐怕是很难得出某个基因显著的结论(因为alpha值会被调得极其小)。FWER的目标是一个错误都不能犯,但实际上我们也许可以容忍在那些我们宣称显著的结果中,有少数其实是犯错的,就看你是不是“宁愿错杀三千,也不放过一个”了。

Efron在前面我提到的文章中把Benjamini和Hochberg在1995年的论文称为“二战以来统计界第二伟大的成果”(他把第一名给了James & Stein的有偏估计),那么B&H做了什么?答案就是虚假发现率(False Discovery Rate,下文简称FDR)。FDR要控制的是在宣称显著的结论中错误的结论的比例,比如10000个基因中有100个基因显著,但实际上有5个是虚假的发现(即本来这5个基因分别在两种处理下的表达并没有差异)。尽管有错误存在,但我们认为可以承受。

初学者到这里应该可以意识到了,通过FDR选出来的结果在理论上就已经存在错误了,当然这只是小问题,更大的问题在于,FDR的定义实际上是一个期望的形式:真实的零假设个数除以被拒绝的零假设个数的期望(零假设是没有差异)。凡是涉及到期望的东西,我们都应该冷静三秒想一下,期望意味着什么?

假设有一个游戏,你获胜的概率是70%,要是赢了,你得到一百万,要是输了,你付出一百万;获利的期望是40万。现在我请你去玩,一把定输赢,你玩不玩?我相信大多数人不会玩(除非你实在太有钱了),为什么期望是赚40万你也不玩?因为期望往往是“样本量无穷大”你才能隐约看到的东西,玩一次游戏,输掉一百万的概率是30%,恐怕你不敢。

FDR是个期望,也就是你做很多次试验,平均来看,FDR在某个数值附近。一次试验中它究竟在哪儿,谁都不知道。就像(频率学派的)置信区间一样,我给你一个区间,其实你并不知道参数在不在区间里,你也无法用概率的方式去描述单个区间,比如你不能说“参数落在这个区间内的概率是95%”(只有无数次抽样重新计算区间,这无数个区间覆盖真实参数的概率才是95%)。

所以,某种意义上,概率论源于赌博,而统计学在骨子里一直都是赌博。旁观者看赌徒,总觉得他在赚钱。当然,统计学家是“高级赌徒”,他们不是随机乱赌。

四、可重复的统计研究

看到这里,你大概也脑子有点抽筋了(如果你把我提到的Lindsey和Efron的文章都看过了的话),我们换个轻松点的话题:可重复的统计研究。这不是我第一次提它,我们一直都在号召可重复的统计研究(如《Sweave:打造一个可重复的统计研究流程》)。还是老话一句,不谈道德问题,因为这是不可控的因素,我们可控的只有制度和工具。源代码不会撒谎。

我们期待学术研究过程的透明化,至少统计之都在努力。

假设检验初步

准备尝试一下,用大白话叙述一遍统计推断中最基础的东西(假设检验、P值、……),算是把这段时间的阅读和思考做个梳理(东西不难,思考侧重在如何表述和展示)。这次打算用一种“迂回的”表达方式,比如,本文从我们的日常逻辑推理开始说起。

0.普通逻辑

复习一下普通逻辑的基本思路。假设以下陈述为真:

你打了某种疫苗P,就不会得某种流行病Q。

我们把这个先决条件表述如下:

如果P 则非Q

其中,

P表示打了疫苗P,

Q表示得流行病Q

或者,更形式化一点:

if P then NOT Q

然后,如果观察到你得了流行病Q,那么就可以推出你没有打疫苗P——这个推断只不过是上述前提条件的逆反命题而已。我们把以上推理过程表述如下:

if P then NOT Q (先决条件)

Q (前提)

———————–

then NOT P (结论)

还有,如果你没有得流行病Q,就能推断出你打了疫苗P吗?显然不能。打疫苗P是不得流行病Q的充分条件,但非必要条件:你没有得流行病Q,可能是因为打了疫苗P,也可以是因为其他任何原因。即,if P then NOT Q,不能够推出if NOT Q then P。

到此为止没有任何令人惊奇的地方。下面将表明,假设检验背后的统计推断规则也只不过是我们以上日常逻辑推理的一个衍生而已。这只需要思维的一次小小的“跳跃”。

1.假设检验

在统计推断中,我们不说“你打了疫苗P,就不会得流行病Q”,而是说,比如,“你打了疫苗P,就有95%的把握不会得流行病Q”,即if P then probably NOT Q。把上面的逻辑推理规则改写成统计推断规则:

if P then probably NOT Q    (先决条件)

Q                                                     (前提)

———————–

then probably NOT P         (结论)

回到以前“万能”的硬币实验,我们做实验来考察一枚硬币是不是均匀的。改写成现在我们熟悉的形式:

P:硬币是均匀的。

Q:在100次投掷中,得到90次正面,10次反面。

我们说,如果是一个均匀的硬币,就不太可能发生这样的情形:投100次,出现90次正面,10次反面(if P then probably NOT Q)。现在如果在100次投掷实验中,观察到出现90次正面,10次反面(Q),那就可以有把握地说,这个硬币不是均匀的(NOT P)。这个推理可以写成与上面一致的统计推断的形式,其中,P是原假设H0,NOT P是备择假设Ha:

H0:硬币是均匀的  (P

Ha:硬币是有偏的 (NOT P

如果原假设为真,即硬币是均匀的,就不太可能发生这样极端的事情,比如:在100次投掷实验中,观察到出现90次正面,10次反面(Q)。如果真的观察到这样极端的事情,你就有把握认为硬币不是均匀的,即拒绝原假设(P),接受备择假设(NOT P)。

另外,如果在100次投掷实验中,观察到60个正面,40个反面(NOT Q)。这时你就不好下结论了,因为一个均匀的硬币可能投出这样的结果,一个有偏的硬币也可能投出这样的结果。最后,你只能说,如果实验结果是这样的,那就没有把握拒绝原假设。这枚硬币是否有偏,需要更多的证据来证明(这通常意味着更多的实验,比如,再投1000次)。

总结一下。在搜集数据之前,我们把想证明的结论写成备择假设,把想拒绝的结论写成原假设。之所以写成这个形式,因为从上面不厌其烦的讨论中得知,这是方便逻辑/统计推断的形式:当我们难以拒绝原假设时,只能得到结论,原假设也许是真的,现在还不能拒绝它;而当我们能够拒绝原假设时,结论是:它就很有把握是不真的。注意,在看到数据之前,我们不知道自己想证明的结论是否能够被证据所支持。

在确定假设检验的形式的同时,我们对之前一直随意说的“把握”、“可能”也做一个限定,即指定一个显著性水平α(significance level),也叫犯第一类错误的概率(type I error,在上面的硬币实验中,就是否定一个均匀硬币的错误,也叫“弃真”错误)。

根据某些保守或稳健的原则(比如,我们认为,把一个无辜的人判决为有罪,比放掉一个有罪的人,后果更为严重),我们要尽量把犯“弃真”错误的概率控制在一个很小的水平里。通常α=0.05,这时候就是说,如果拒绝了原假设,你就有95%的把握说原假设是不真的。这里,95%(=1-α)就是置信水平(confidence level)。

又,放掉一个有罪的人,即把一个有罪的人判为无罪,这犯的是第二类错误β(type II error,在硬币实验中,就是把一个有偏的硬币当成均匀硬币的错误,也叫“取伪”错误)。关于第一类和第二类错误之间的权衡取舍(trade off),详见《决策与风险》。在我们的假设检验里,我们认为犯一类错误的后果比犯第二类错误的后果更为严重。

需要注意的是,在这里,我强调的是先提出需要检验的假设,然后再搜集收据。这是统计推断的原则之一。如果看到了数据之后再提出假设,你几乎可以得到所有你想要的结果,这是不好的机会主义的倾向。强调这些,是因为在学校里,我们大多是看了别人搜集好的数据之后再做统计练习。

事先确定好你想拒绝/证明的假设,在看到数据之前,你不知道结果如何。

2.P值(P Value)

上面提到“极端”事件,比如,在100次硬币投掷实验中,观察到出现90次正面,10次反面(Q)。怎么样的事件才是“极端的”?简单地说,一个事件很极端,那么少比它本身“更极端”的事件就非常少(比如,只有“91次正面,9次反面”、“91次反面,9次正面”等情况才比它更极端)。

但这个Q只是从一次实验中得出的。我们可以重复做这个实验,比如100次,每次都投掷100次,记录下的正面数X,它构成一个二项分布,X~B(n,p),其中,n=100,p=0.5。根据某个中心极限定理,正态分布是二项分布的极限分布,上面的二项分布可以由均值为np=50,方差为np(1-p)=25的正态分布来近似。我们在这个近似的正态分布的两端来考察所谓“更极端”的事件,那就是正面数大于90或者小于10。

重复一遍,“P值就是当原假设为真时,所得到的样本观察结果更极端的结果出现的概率”。如果P值很小,就表明,在原假设为真的情况下出现的那个分布里面,只有很小的部分,比出现的这个事件(比如,Q)更为极端。没多少事件比Q更极端,那就很有把握说原假设不对了。

在上述近似的正态分布中,P值就等于X>90 或 X<10的概率值(记做,P{X>90 or X<10})。根据对称性,这个概率值等于2*P{X<10}=1.2442E-15。

上面我们的确求出了一个非常小的P值,但如何不含糊地确定它就是很“极端”呢? 事先确定的显著性水平α,本身就是一个判定法则。只要P值小于显著性水平α,我们就认为,在认为原假设为真的情况下出现的事件Q,是如此地极端,以至于我们不再相信原假设本身。一句话,我们的判定法则是:

P值小于显著性水平α,拒绝原假设。

3.一个手算示例

用一个双侧的单样本T检验做例子。假设我们想知道,螃蟹的平均温度,跟空气的温度(24.3)有没有统计差别(α=0.05)。事先确定的假设检验的形式表达如下:

零假设(H0):   μ=24.3°C

备择假设(Ha):  μ≠24.3°C

以下是25只螃蟹在温度为24.3°C下的体温(单位:°C):

25.8    24.6    26.1    22.9    25.1
27.3    24        24.5    23.9    26.2
24.3    24.6    23.3    25.5    28.1
24.8    23.5    26.3    25.4    25.5
23.9    27        24.8    22.9    25.4

一些基本的算术结果:

样本均值:$\bar{X}=25.3$

样本量:n=25

样本方差:$s^2$=1.8

样本均值的标准误差:$s(\bar{X})=\sqrt{s^2/n}=0.27$

这里T检验的思路如下:

  1. 我们先假设H0为真,即认为螃蟹的平均温度跟空气温度没有差异(P),  μ=24.3°C。有一个极端事件Q,如果原假设H0成立,Q就不成立(if H0 then probably NOT Q);但如果在原假设为真的情况下,出现了这么一个Q,那我们就有把握拒绝原假设。
  2. 样本均值:$\bar{X}$是总体均值μ的最好的估计,在本例中,$\bar{X}=25.03$。
  3. 这个样本均值只是一个估计值。它只是从总体的一个随机样本中得到的(样本是上述25只螃蟹)。我们不知道这次实验结果是不是“极端”事件。而判断一个事件是不是极端事件,根据第二节的讨论,我们可以重复做上述实验,比如100次,每次都抓25只螃蟹,都在空气温度24.3的状态下测量其体温,然后也各自求出一个样本均值来。
  4. 容易得出,这种实验出来样本均值,辅以适当的数学形式,就服从一个自由度为24(=25-1)的t分布,即$(\bar{X}-\mu)/s(\bar{X})\sim t(24)$。
  5. 样本均值$\bar{X}=25.03$,在这个自由度为24的t分布下,有一个对应的t值,t=25.03-24.3/0.27=2.704。现在我们可以在整个分布里考察这个t值。在这个自由度为24的t分布里,我们看 t=2.704是不是一个“极端”事件Q。根据对称性,比Q更极端的是那些大于2.704或者小于-2.704的点。

t

从上图可以看到,在这个t分布里,比t=2.704更“极端”的点占整个分布的0.0124。这个0.0124就是我们要求的P值。这个P值小于我们事先选定的显著性水平α=0.05,因此我们可以拒绝原假设,认为这批螃蟹的平均体温不等于空气温度。

这个双侧P值可以手算如下:

在SAS里,P=2*(1-probt(t,df))=2*(1-probt(2.704,24))=0.012392

在R里,     P=2*(1-pt(t,df))=2*(1-pt(2.704,24))=0.012392

———-

以上是用P值作为判定条件。一个等价的做法是用临界值来判断。我们事先给定的显著性水平α=0.05,在这个自由度为24的t分布里,就对应着一个临界t值2.064。下图的阴影部分,也称作拒绝区域。上面求出的跟样本均值$\bar{X}=25.03$对应的t值=2.704,处在这个拒绝区域内(2.704>2.064),于是我们一样拒绝原假设。

t2

又,上述临界值可以手算(或查表)如下:

在SAS里,tCritic=tinv(1-alpha/tail,df)=2.06390

其中,alpha=0.05,tail=2表示双侧检验,df=24.

在R里,tCritic=qt(1-alpha/tail,df)=2.063899

4.注

本文是对近期阅读做的一个笔记。作为一个非统计科班出身的程序员,我一直在思考,如何来理解统计概念,以及如何把自己的理解向同行传达。关于用日常逻辑推理来理解假设检验的思路,来自

Common Statistical Methods for Clinical Research with SAS Examples(2nd edition, SAS Inc., 2002, by Glenn A. Walker)

关于决策与风险的讨论,参考了

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

第三节示例的数据,来自

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

第三节的t分布图,来自一个在线的t分布生成器(很好用):

http://onlinestatbook.com/analysis_lab/t_dist.html

附录: 用SAS来计算

上面的文字尽量做到“平台无关”。这里附出SAS例子,是想把以上的手算结果跟机器结果做个对照,让读者更有信心一些。 欢迎读者贴出自己趁手的工具得出的结果。

/*data*/
data body;
input temp @@;
h0=24.3;
diff=temp-h0;
datalines;
25.8    24.6    26.1    22.9    25.1
27.3    24      24.5    23.9    26.2
24.3    24.6    23.3    25.5    28.1
24.8    23.5    26.3    25.4    25.5
23.9    27      24.8    22.9    25.4
;

/*method 1: use proc means*/
proc means data=body T PRT;
var diff ;
run;

结果是:

t Value    Pr > |t|
——————-
2.71      0.0121
——————-

上面的t Value 就是计算出来的t值,Pr > |t| 就是P值(这里的|t|就是上面计算出来的t值2.704,Pr > |t|求的是比t值更极端的概率,即P值)。proc means没有提供临界t值(即通常说的查表得出的t值),下同。

/*method 2 (prefered): use proc ttest*/
proc ttest data=body h0=24.3 alpha=0.05;
var temp;
run;

proc ttest的结果更为丰富:

N      Mean     Std Dev  Std Err    Minimum   Maximum

25     25.0280      1.3418      0.2684 22.9000        28.1000

Mean     95% CL     Mean       Std Dev     95% CL   Std Dev

25.0280 24.4741  25.5819            1.3418             1.0477   1.8667

DF    t Value    Pr > |t|

24       2.71           0.0121

统计学博文导读:内贾德大选作弊?流星撞飞机的概率?买双色球?

“统计之都”站的“网站导读”栏目的设立是为了以简短的形式向大家介绍一些有意思而且有水平的统计学文章,不求理论之复杂,但求统计学之生活化,让大家看到一些统计学的“另类”面目。若这个目的达不到,那么我希望大家读完这些导读文章之后能说一句“哇,原来统计不是会计啊/不是做报表的啊/不是数学啊”也足够了。另外,现在网上很多文章都是抄来抄去(更恶劣的是不加出处的抄袭),我们觉得这种做法极其无聊,是对原始作者的极大不尊重,也容易造成以讹传讹误导不明真相的围观群众,本站这个栏目的建立,也是基于这一点考虑之上提供一种“引用他人文章”的示例,很傻很天真地希望互联网的抄袭现象能够有所收敛。言归正传:

一、用数字作弊要小心:关于内贾德的选票

内贾德日前伊朗的大选可谓轰轰隆隆,颇引人注目,内贾德在胜出之后却引来一片质疑,今日又有传闻说穆萨维才是真正的胜出者。总之疑云重重,那么让政治家玩政治家的游戏吧,我们从另一个视角来关注一下这次选举。密歇根大学的教授Walter R. Mebane, Jr.这几天一直在分析选票数据,今天的文章参见:Note on the presidential election in Iran, June 2009(注意论文和数据以及R代码都在更新中,如果不能访问,请到他的主页上找)

文章主要基于两点理论去检验选票数据:

  1. 本福特定律(Benford’s Law):生活中的数据里,1~9这9个数字在首位的出现并非均匀分布,例如1出现在一个数字的首位的概率约为1/3,而不是想象中的1/9,越往后的数字出现概率越低。不过,作者得到的选票数据是按地区汇总的,而对于汇总数据,我们往往难以发现它作弊的嫌疑,因为汇总数据倾向于符合本福特定律,从这一点上,作者没有找到足够的证据证明选票数据作弊;
  2. 检验离群点:手头的数据只有选票数,没有其它变量,这种情况下统计建模的局限性很大,似乎各种工具都施展不开,不过作者还是在“艰苦”的条件下建立了过度散布的二项回归模型(overdispersed binomial regression),因变量是二分变量(内贾德 vs 穆萨维),自变量是选票数,看看哪些地区的选票数有离群点出现,熟悉伊朗政局情况的人从这些离群点可能发现违背常理的现象(选票过高?过低?)

看来统计学家总是有办法拷问数据,不知咱国内是否有这样的统计学家呢?

想起本科上课的时候我们有一位老师提到陈毅元帅的一句话:“莫伸手,伸手必被捉。”在数据上,还是莫造假吧,只要有人较真,一定有办法看出来的。

二、飞机被流星撞到的概率多大:法航失事之后的计算

前几天,本站作者、COS论坛元老刘思喆用R计算了一下飞机被流星击中而失事的概率,发现一架飞机在11小时飞行过程中被流星击中的概率是双色球中一等奖概率的1/100,看样子可以舒口气,不过坏消息是,20年中被击中的概率就陡然上升到5%了。看来还是双腿走路可靠……

然而我想知道的是:

  • 过去20年有飞机被流星击中过么?
  • Poisson分布合适么?或者,流星砸向地球的分布是什么?有没有天文学相关的实证呢?

三、让R帮你查看是否中了500万大奖:“懒惰”的彩民

这一篇依旧是把统计和R融入生活的刘思喆:500万?去买双色球!(被评论为标题党)彩民们每天眼巴巴对着彩票网以及自己的彩票的若干位数字看,实在是很辛苦,因此思喆老大体谅大家,挑灯夜战写了一段R代码,按照该文的描述,彩民们每天只要打开R,就知道自己有没有中奖。我提议,要是哪位读者因此中奖了一定要给思喆老大提成以及给统计之都捐赠

我后来的想法是,干脆写个批处理文件(R CMD BATCH),开机自动运行好了,连R都不用打开。懒惰是创新之源,这话没错。

R还有很多“奇怪”的应用,在作者的主页贝吉塔行星中能找到更多。

四、其它文章

  • 画函数曲线本来是一件简单的事情,为什么我们总是忘记初中的“描点法”呢:画曲线的通用办法:描点法画图(谢益辉)
  • 标签云将词语文本的大小与其某种属性(例如重要性、出现频率)关联起来,因此标签云图可以直观展示一些词语的属性。例如文本越大表示出现频率越高,那么一眼看去,最大的词就是最热频词了,这里介绍了一种用R生成标签云的方法:Creating Tag Cloud Using R and Flash / JavaScript (SWFObject)(谢益辉)
  • 最后强烈建议对统计软件感兴趣(尤其是对R)的同志们订阅Journal of Statistical Software的RSS,你会经常发现一些稀奇古怪的论文:http://www.jstatsoft.org/rss(啥?你到现在还不知道RSS是神马东西?快来人把这个火星观众撵出去)

若您平时看到觉得写得好的文章,请勿忘推荐给我们,联系邮箱:[email protected]