标签归档:正态分布

标准正态分布函数的快速计算方法

标准正态分布的分布函数 $\Phi(x)$ 可以说是统计计算中非常重要的一个函数,基本上有正态分布的地方都或多或少会用上它。在一些特定的问题中,我们需要大量多次地计算这个函数的取值,比如我经常需要算正态分布与另一个随机变量之和的分布,这时候就需要用到数值积分,而被积函数就包含 $\Phi(x)$。如果 $Z\sim N(0,1), X\sim f(x)$,$f$ 是 $X$ 的密度函数,那么 $Z+X$ 的分布函数就是

$$P(Z+X\le t)=\int_{-\infty}^{+\infty} \Phi(t-x)f(x)\mathrm{d}x$$

我们知道,$\Phi(x)$ 没有简单的显式表达式,所以它需要用一定的数值方法进行计算。在大部分的科学计算软件中,计算的精度往往是第一位的,因此其算法一般会比较复杂。当这个函数需要被计算成千上万次的时候,速度可能就成为了一个瓶颈。

当然有问题就会有对策,一种常见的做法是略微放弃一些精度,以换取更简单的计算。在大部分实际应用中,一个合理的误差大小,例如 $10^{-7}$,一般就足够了。在这篇文章中,给大家介绍两种简单的方法,它们都比R中自带的 pnorm() 更快,且误差都控制在 $10^{-7}$ 的级别。

第一种办法来自于经典参考书 Abramowitz and Stegun: Handbook of Mathematical Functions
公式 26.2.17。其基本思想是把 $\Phi(x)$ 表达成正态密度函数 $\phi(x)$ 和一个有理函数的乘积。这种办法可以保证误差小于 $7.5\times 10^{-8}$,一段C++实现可以在这里找到。(代码中的常数与书中的略有区别,是因为代码是针对误差函数 $\mathrm{erf}(x)$ 编写的,它与 $\Phi(x)$ 相差一些常数)

我们来对比一下这种方法与R中 pnorm() 的速度,并验证其精度。

library(Rcpp)
sourceCpp("test_as26217.cpp")

x = seq(-6, 6, by = 1e-6)
system.time(y 

可以看出,A&S 26.2.17 的速度大约是 pnorm() 的三倍,且误差也在预定的范围里,是对计算效率的一次巨大提升。

那么还有没有可能更快呢?答案是肯定的,而且你其实已经多次使用过这种方法了。怎么,不相信?看看下面这张图,你就明白了。

normal_table

继续阅读标准正态分布函数的快速计算方法

漫谈正态分布的生成

本文作者简介:王夜笙,就读于郑州大学信息工程学院,感兴趣的方向为逆向工程和机器学习,长期从事数据抓取工作(长期与反爬虫技术作斗争~),涉猎较广(技艺不精……),详情请见我的个人博客~
个人博客地址:http://bindog.github.io/blog/

感谢怡轩同学的悉心指导~

之前拜读了靳志辉(@rickjin)老师写的《正态分布的前世今生》,一直对正态分布怀着一颗敬畏之心,刚好最近偶然看到python标准库中如何生成服从正态分布随机数的源码,觉得非常有趣,于是又去查找其他一些生成正态分布的方法,与大家分享一下。

利用中心极限定理生成正态分布

设$X_1,X_2,\cdots ,X_n$为独立同分布的随机变量序列,均值为$\mu$,方差为$\sigma^2$,则

$$Z_n=\frac{X_1+X_2+\cdots+X_n-n\mu}{\sigma \sqrt n}$$

具有渐近分布$N(0,1)$,也就是说当$n \rightarrow \infty$时,

$$P\left \{ \frac{X_1+X_2+\cdots+X_n-n\mu}{\sigma \sqrt n} \leq x \right \} \rightarrow \frac{1}{\sqrt{2\pi} } \int_{-\infty }^{x} e^{ -\frac{t^2}{2} } \, dt$$

换句话说,$n$个相互独立同分布的随机变量之和的分布近似于正态分布,$n$越大,近似程度越好。当然也有例外,比如$n$个独立同分布的服从柯西分布随机变量的算术平均数仍是柯西分布,这里就不扩展讲了。

根据中心极限定理,生成正态分布就非常简单粗暴了,直接生成n个独立同分布的均匀分布即可,看代码

继续阅读漫谈正态分布的生成

正态分布的前世今生(下)

6. 开疆拓土,正态分布的进一步发展

19世纪初,随着拉普拉斯中心极限定理的建立与高斯正态误差理论的问世,正态分布开始崭露头角,逐步在近代概率论和数理统计学中大放异彩。在概率论中,由于拉普拉斯的推动,中心极限定理发展成为现代概率论的一块基石。而在数理统计学中,在高斯的大力提倡之下,正态分布开始逐步畅行于天下。

6.1 论剑中心极限定理

先来说说正态分布在概率论中的地位,这个主要是由于中心极限定理的影响。 1776 年,拉普拉斯开始考虑一个天文学中的彗星轨道的倾角的计算问题,最终的问题涉及独立随机变量求和的概率计算,也就是计算如下的概率值
$$ S_n = X_1 + X_2 + \cdots + X_n $$
$$P(a < S_n < b) = ? $$

在这个问题的处理上,拉普拉斯充分展示了其深厚的数学分析功底和高超的概率计算技巧,他首次引入了特征函数(也就是对概率密度函数做傅立叶变换)来处理概率分布的神妙方法,而这一方法经过几代概率学家的发展,在现代概率论里面占有极其重要的位置。基于这一分析方法,拉普拉斯通过近似计算,在他的1812年发表的名著《概率分析理论》中给出了中心极限定理的一般描述:

定理:[拉普拉斯, 1812]  $ e_i (i=1, \cdots n)$ 为独立同分布的测量误差,具有均值$\mu$ 和方差 $\sigma^2$。如果 $\lambda_1, \cdots, \lambda_2$ 为常数, $a>0$, 则有
$$ \displaystyle P\left(\left|\sum_{i=1}^n \lambda_i(e_i – \mu)\right|
\le a \sqrt{\sum_{i=1}^n \lambda_i^2}\right)
\approx \frac{2}{\sqrt{2\pi}\sigma} \int_0^a e^{-\frac{x^2}{2\sigma^2}} dx . $$

这已经是比棣莫弗-拉普拉斯中心极限定理更加深刻的一个结论了,理科专业的本科生学习《概率论与数理统计》这门课程的时候,通常学习的中心极限定理的一般形式如下:

[林德伯格-列维 中心极限定理] 设$X_1,\cdots, X_n$ 独立同分布,且具有有限的均值 $\mu$ 和方差 $\sigma^2$ ,则在 $n \rightarrow \infty$ 时,有
$$ \displaystyle \frac{\sqrt{n}(\overline{X} – \mu)}{\sigma} \rightarrow N(0,1) .$$

多么奇妙的性质,随意的一个概率分布中生成的随机变量,在序列和(或者等价的求算术平均)的操作之下,表现出如此一致的行为,统一的规约到正态分布。

central_limit_theorem中心极限定理

概率学家们进一步的研究结果更加令人惊讶,序列求和最终要导出正态分布的条件并不需要这么苛刻,即便 $X_1,\cdots, X_n$ 并不独立,也不具有相同的概率分布形式,很多时候他们求和的最终的归宿仍然是正态分布。一切的纷繁芜杂都在神秘的正态曲线下被消解,这不禁令人浮想联翩。中心极限定理恐怕是概率论中最具有宗教神秘色彩的定理,如果有一位牧师拿着一本圣经向我证明上帝的存在,我是丝毫不会买账;可是如果他向我展示中心极限定理并且声称那是神迹,我可能会有点犹豫,从而乐意倾听他的布道。如果我能坐着时光机穿越到一个原始部落中,我也一定带上中心极限定理,并劝说部落的酋长把正态分布作为他们的图腾。

继续阅读正态分布的前世今生(下)

正态分布的前世今生(上)

神说,要有正态分布,就有了正态分布。
神看正态分布是好的,就让随机误差服从了正态分布。
创世纪—数理统计

1. 正态分布,熟悉的陌生人

学过基础统计学的同学大都对正态分布非常熟悉。这个钟形的分布曲线不但形状优雅,它对应的密度函数写成数学表达式
$$ \displaystyle f(x)=\frac{1}{\sqrt{2\pi}\sigma}e^{-\frac{{(x-\mu})^2}{2\sigma^2}} $$
也非常具有数学的美感。其标准化后的概率密度函数
$$ \displaystyle f(x)=\frac{1}{\sqrt{2\pi}}e^{-\frac{x^2}{2}} $$
更加的简洁漂亮,两个最重要的数学常量 $\pi$、$e$ 都出现在这公式之中。在我个人的审美之中,它也属于 top-N 的最美丽的数学公式之一,如果有人问我数理统计领域哪个公式最能让人感觉到上帝的存在,那我一定投正态分布的票。因为这个分布戴着神秘的面纱,在自然界中无处不在,让你在纷繁芜杂的数据背后看到隐隐的秩序。

normal_curve正态分布曲线

正态分布又通常被称为高斯分布,在科学领域,冠名权那是一个很高的荣誉。2002年以前去过德国的兄弟们还会发现,德国1991年至2001年间发行的的一款10马克的纸币上印着高斯(Carl Friedrich Gauss, 1777-1855)的头像和正态密度曲线,而1977年东德发行的20马克的可流通纪念钢镚上,也印着正态分布曲线和高斯的名字。正态分布被冠名高斯分布,我们也容易认为是高斯发现了正态分布,其实不然,不过高斯对于正态分布的历史地位的确立是起到了决定性的作用。

10dm_with_gauss_curve 10dm_with_gauss_curve_detail   20-mark-gauss
德国马克和纪念币上的高斯头像和正态分布曲线

正态曲线虽然看上去很美,却不是一拍脑袋就能想到的。我们在本科学习数理统计的时候,课本一上来介绍正态分布就给出分布密度函数,却从来不说明这个密度函数是通过什么原理推导出来的。所以我一直搞不明白数学家当年是怎么找到这个概率分布曲线的,又是怎么发现随机误差服从这个奇妙的分布的。我们在实践中大量的使用正态分布,却对这个分布的来龙去脉知之甚少,正态分布真是让人感觉既熟悉又陌生。直到我读研究生的时候,我的导师给我介绍了陈希儒院士的《数理统计学简史》这本书,看了之后才了解了正态分布曲线从发现到被人们重视进而广泛应用,也是经过了几百年的历史。

正态分布的这段历史是很精彩的,我们通过讲一系列的故事来揭开她的神秘面纱。

继续阅读正态分布的前世今生(上)

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

最近在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:打造一个可重复的统计研究流程》)。还是老话一句,不谈道德问题,因为这是不可控的因素,我们可控的只有制度和工具。源代码不会撒谎。

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