标签归档:直方图

从中心极限定理的模拟到正态分布

昨日翻看朱世武老师的《金融计算与建模》幻灯片(来源,幻灯片“13随机模拟基础”),其中提到了中心极限定理(Central Limit Theorem,下文简称CLT)及其SAS模拟实现。由于我一直觉得我们看到的大多数对CLT的模拟都有共同的误导性,因此在此撰文讲述我的观点,希望能说清楚CLT的真实面目,让读者对它有更深刻的理解。

一、中心极限定理说了什么

从广义的角度来讲,CLT说的是一些随机变量之和(在n趋于无穷的时候)趋于正态分布,条件是这些随机变量之间相关性不要太严重,每个随机变量自身的方差相对来说不要太大。我们平时看到的都是CLT的最简单的形式:随机变量独立同分布,方差有限。实际上独立同分布不是必要的条件。要理解这段话,我们就得去思考一下Lindeberg条件或Liapounov条件的含义,以及非独立情况下的CLT。而CLT为什么偏偏选中了正态分布呢?这主要是泰勒展开的功劳(注1),加上特征函数作为“催化剂”,捣鼓半天,最后一看,特征函数的极限形式正是正态分布。本文不打算走到这样偏僻的数学角落,只是先介绍一下CLT的前世今生,当我们在统计理论上犯迷糊的时候,回溯到事情的最原始状态,往往能让头脑更清醒(注2)。

即使是最简单的CLT,我觉得几乎所有做模拟的人都走入了一个误区,就是把钟形密度曲线解释为正态分布。给一个非正态的总体分布,重复计算若干次样本均值,把密度以直方图或密度曲线的形式画出来,给观众们说,看:钟形曲线!群众们一看,哇,原本是左偏/右偏/水平的密度,现在真的变成了对称的钟形曲线。于是乎不得不信服CLT的魔力。

二、中心极限定理的模拟

我们缺的不是钟形曲线,而是样本均值的分布和正态分布差多少的指标。这一点在下面第3小节中再谈,先看几个常见的模拟例子。

1、对称的钟形曲线不代表正态分布

很多人喜欢用掷骰子的例子来讲中心极限定理,大意是将两次独立掷骰子的结果加起来,看这个和的分布。我犹记当年人大一位老师给我们上抽样就举了这个例子说明CLT的魔力:看,即使样本量为2,得到的分布也是正态分布!恰好昨天在朱老师的幻灯片中又看到这个例子,不禁感叹这样一个糟糕的例子竟然经久不衰。这个例子的骗局在哪儿呢?其实很简单:在于掷骰子得到的结果的“对称性”,换句话说,结果是1、2、3、4、5、6,这6个数字围绕其均值3.5左右对称,因此两次骰子的结果加起来也围绕7对称,再画个图,两边低,中间高,看似正态分布好像出来了。朱老师的SAS代码如下:

data a;
do x1=1 to 6;
do x2=1 to 6;
output;
end;
end;  /*模拟掷骰子两次,生成36行数据*/
data a;
set a;
x=sum(x1,x2);
proc univariate data=a noprint;
var x;
histogram/normal (mu=est sigma=est);
run;

先跑题说一处细节:这段代码中的直方图不合适,因为我们知道结果x只可能是整数,而(SAS)直方图默认的分组也会以整数为边界,每个边界上的整数都被归到右边的条中去了,而这个所谓的直方图其实只是个条形图:展示每种结果的频数。所以这种情况下认为指定非整数的分组边界才能显示数据真实的分布。以下是R代码:

x = rowSums(expand.grid(1:6, 1:6))
hist(x, breaks = seq(min(x) - 0.5, max(x) + 0.5, 1), main = "")
两次掷骰子得到的结果之和的分布
两次掷骰子得到的结果之和的分布

样本量为2的时候真的得到正态分布了么?CLT确实有它的神奇之处,但还没那么神奇。以上结果如我前文所说,仅仅是由于样本空间中的元素的对称性,所以得到了一副对称的图形,看起来像是正态分布,如果我们再用另外6个数字试验一下,马上就看穿这种迷惑了,我把最后一个数字换成15再看直方图:

d = c(1, 2, 3, 4, 5, 15)
x = rowSums(expand.grid(d, d))
hist(x, breaks = seq(min(x) - 0.5, max(x) + 0.5, 1), main = "")
6个数字中2个样本之和的分布
6个数字中2个样本之和的分布

如何?“正态分布”哪里去了?客官是否能看清骰子中的中心极限定理谎言了?

2、直方图上添加正态密度曲线有误导性

在我还会用SPSS的年代(大约是SPSS 15.0),我就发现SPSS的一个荒唐之处,它的直方图有个选项,可以控制是否添加正态密度曲线。我们被“正态”毒害得多深,从这些软件设置就可以看出来。为什么只能加正态密度曲线,而不能加数据自身反映出来的核密度估计曲线?换句话说,数据的分布一定要是正态么?

我们看CLT的模拟,很多人也喜欢在直方图上加上正态分布密度曲线,这有一定道理:可以看直方图跟正态密度的差异有多大。然而,我们却很少见直方图上的核密度曲线(注3)。既然是要作比较,就拿最可比的去比,比如曲线对曲线,直方图对直方图,人眼本来就是不精确的测量工具,那么制图者就应该提供尽量准确而方便的参照系。

3、给出拟合好坏的指标

综上,我在两年多以前已经把这些CLT模拟的想法写成函数收录在R包animation中,参见函数clt.ani()。代码及输出如下:

if (!require(animation)) install.packages("animation")
library(animation)
ani.options(interval = 0.1, nmax = 100)
par(mar = c(4, 4, 1, 0.5))
clt.ani()
中心极限定理模拟:从指数分布到正态分布
中心极限定理模拟:从指数分布到正态分布

这个CLT模拟的过程很简单:给一个总体分布(默认为右偏的指数分布),在给定样本量n时不断重复抽样分别计算样本均值,一直这样计算obs个均值,并画出它们的直方图和相应的核密度估计曲线;然后随着n增大,看相应的样本均值分布如何。此外,我使用了Shapiro正态性检验来检验这些均值的正态性,并把P值取出来画在下半幅图中。这样我们就很清楚地知道,对于每一种样本量(n = 1, …, 100),我们的样本均值究竟离正态分布多远。此处P值就充当了一个拟合好坏的指标。可以看出,上面的动画中,当样本量n超过20之后,P值会普遍偏大,也就是样本均值的分布和正态分布比较接近(严格来说,是“不能拒绝正态分布”),但也不能保证样本量大就一定意味着正态分布,譬如上图中n=40的时候P值就很小。

三、小结

正态分布在统计中的地位如此之重要,以至于人们几乎认为正态是一种自然而然的分布,本文想说明的是,正态有它的自然性(参见注1),但我们不能逮着两边低中间高的东西都叫正态(参见小学课文《小蝌蚪找妈妈》);做模型或分析数据之前,先清空脑子里的这种“正态教义”,用事实说话。

另外,前面胡江堂挖了个大坑,大家跳进去争论了半天SAS和R的问题,我在这里也挖个小坑:就模拟而言,如果有SAS用户愿意研究一下,我想知道朱老师的SAS模拟代码是否有改进的余地。SAS自70年代创立,过了二十多年才引进“函数”这种杀人越货编程必备之工具,一直以来都是“宏”的天下,说它“恐龙”应该也不冤枉吧。朱老师幻灯片中的长篇SAS代码,如果用R改写,应该都不会超过三五行甚至一行(注4),而反过来可能就麻烦了,比如上面的动画(用SAS怎么写?),它用R之所以方便,是因为R大量的统计相关的函数可任意调用,而图形也很灵活,可以把P值动态写在图例的位置,也可以愿意把样本量n作为下标动态写在$\bar{X}$旁边。

脚注

注1:以我的浅见,统计学中绝大多数极限正态分布的来源都是泰勒展开,主要是因为泰勒展开中有一个平方项,这个东西和正态分布的密度函数(或特征函数)看起来形式相似,再加上一些对高次项的假设(使它们在极限中消去),正态分布也就来了。另一个典型例子是极大似然估计的渐近正态性:泰勒展开中一阶导数为零,剩下二阶像正态。此外,正态分布有个独特的特点,就是它的密度函数和特征函数长得很像,冥冥之中也主导了很多统计量的分布(这一点我还没太想清楚)。

注2:我很少见统计学家强调这种想法的重要性,仿佛大家都埋在公式堆里都推导得倍儿高兴。有两个例外:一个是Brian Ripley教授,我在看他的一份幻灯片Selecting amongst large classes of models的时候发现他竟然回顾了AIC的假设条件——这是我见过的唯一一个讲模型选择时会回到数学根源的人;另一个是Ripley的老搭档,Bill Venables,他在他那著名的手稿Exegeses on Linear Models中竟然从泰勒展开来说线性模型的来历,这也是我在所有我读过的线性模型相关的文章和书中看到的唯一一份从泰勒展开角度谈线性模型的材料。我个人并不喜欢数学,但我喜欢看思路。

注3:我的观点是,画直方图尽量加上核密度估计曲线,因为它能刻画数据在任意位置上的密度大小,而直方图的形状则完全受制于分组边界的位置。当然,核密度估计曲线也取决于核函数,但大多数情况下,核函数不会对曲线的形状有太大的影响。不过邱怡轩在前面“有边界区间上的核密度估计”一文中提到的问题很值得注意。

注4:伯努利随机变量之和的分布。SAS代码:

symbol;
goptions ftext= ctext= htext=;

%macro a(n);
data rv;
retain _seed_ 0;
do _i_ = 1 to &n;
binom1 = ranbin(_seed_,&n,0.8); /*ranbin函数返回以seed、n、0.8为参数的二项分布随机变量*/
output;
end;
drop _seed_ _i_;
run;
%mend;

%a(100);

proc univariate data=rv noprint;
var binom1;
histogram/normal( mu=est sigma=est);
inset normal ;
run;

R代码:

a = function(n) hist(rbinom(n, n, 0.8))

什么叫向量化编程?这就是R它爹S语言诞生的原因之一:统计学的编程不应该涉及到那么多繁琐的循环。

漫谈相关与回归

老师不断提醒我要对统计学的基本概念、定义及背景反复思考,这样才不会本末倒置,迷失方向。但是这个做起来很难,因为那些概念定义等看起来实在”太简单”、”没什么东西”,可能还是不能够平心静气吧!

最近静下来看了David Freedman等著的《统计学》的”相关与回归”部分,以及一篇关于直方图的文章,不免有些感慨!其实统计学中的很多概念、工具、方法等的实际意义或作用可能要比我们认为的要大很多,同时,当我们从一些概念定义等中发现出一些新东西时我们总会欣喜若狂。世界上的很多事物又何尝不是如此,人们对事物的了解总易受到传统或他人的影响仅仅停留在表面,很少达到全面而深刻,而一旦我们获得了那种深刻的洞察力,才发现真实世界是何等的精彩!一直以为直方图很简单,无非是一些代表频数的柱状图的组合而已,感觉没什么作用,但是看了一篇关于直方图制作方面的论文时,才认识到直方图的威力。直方图其实是非参数统计中估计总体分布特征的一项重要工具,选择好适当的组距和边界点(组距和最小边界点是关键),随着样本量的增大,它可以非常接近地反映数据的真实分布情况。其实,在统计中使用一种工具方法的目的也应该是使现有的数据尽可能多地反映出真实的信息,而这项工作往往是一个无底洞(这时又要考虑到效率问题了)。

散点图亦是如此。散点图给出了所有数据点的信息,但是如何从这些数据中获得结论或拟合模型,甚至用来预测?面对一张散点图,相关和回归应该是最容易想到的吧!这里主要谈谈两个变量间的相关和回归。

在研究两个变量的关系时,一般会先看看它们的散点图,在图中两变量的关系还是比较直观的,大致可以判断是否线性相关及相关性大小如何,是否是非线性相关等。而到底什么是相关呢?相关其实就是知道一件事对了解另一件事的帮助的大小。实际中,如果对某一事物不太了解,但是对与其有一定联系的另一事物有所了解,如果这种联系很强,那我们对于那件不了解的事物就有了更多的信息,或者说对这个不了解的事物有了更大的自信去预测。其实这也是研究中的一种常用的方法。

关于两个变量间的相关系数的计算。我们都知道两个变量X与Y的相关系数的计算公式为$Cov(x,y)/(SD(x)*SD(y))$,然而这已经是一个结果性东西了,我更推崇David Freedman等著的《统计学》中计算方法:先分别对两个变量做标准化,比如对变量X做标准化$(x_i-\bar{x})/SD(x)$,然后对应的标准量相乘,最后加总再求平均. 这种求法反映到散点图中,相当于对散点图的坐标刻度标准化,从而使两个坐标轴具有了相同的刻度,同时在直观两个变量之间的相关性大小时不会受到各自的标准差大小的影响。这个新的坐标系把所有的点(数据对)分到了不同的象限,通过观察各个象限的点的个数和大致分布情况便可以对相关性的大小与正负有直观的了解,比如更多的点都分布在一、三象限且群集于一条直线周围,那么这两个变量的一般具有较强的正线性相关。

我们都知道相关系数是-1到1之间的一个实数,那么相关系数为0.8是不是表示百分之八十的点群集在一条直线的周围吗?当然不是,相关系数是基于全体数据的一个综合信息,它反映的是所有点与某一条直线的群集程度,而不是一部分的点。由此也不免想到,我们在用到一些概念或定义时,也必须清楚这个概念或定义是基于怎样的对象,或有哪些局限条件或假定,比如概率论中的”事件”, “事件”是基于特定条件的,在具体使用过程中大家对这些特定条件都太”熟悉”以至于很少关注到它们,所以一旦条件改变了,大部分的人认为还是指同一个东西。

相关与因果也是这样,我们都知道相关不能同因果划等号,但实际过程中人们总”自然而然”地得出一些结论。一般来说,体重和升高相关,那体重大是不是就因为身高高呢?除非有一个固定的身材标准,而所有人都是这个标准。(随便提个问题:如果两个变量严格线性相关,即相关系数为1,那是否可以说这两个变量中一个为因一个为果呢?)。其实两个变量的相关更经常的情况是它们同时受到另外的一个或多个因素的影响,在这里可以通过对照试验或观察研究来进一步研究。另外,相关是可逆的,而因果则不可以。所以我们分析相关时总是如此谨慎地说,某某变化,与此相关的某某”相应地”如何变化。研究发现,个人收入与教育水平相关,高教育水平是不是高收入的原因呢?实际情况是它们相互影响:教育水平高的人收入一般较高,收入高的一般也更有能力获得继续教育的机会。虽然相关不是因果,然而有时我们并不需要弄清所有的因果关系,盯住输入和输出,只要存在相关,即使不是因果关系也不妨碍人们利用这种关系来进行推断。比如利用公鸡打鸣来预报太阳升起,虽然公鸡打鸣绝对不是日出的原因(虽然打鸣发生在先)。

在对两变量的相关关系有一定了解后,接下来的自然想法便是拟合回归模型。”回归”这一词来自于高尔顿的父子两代身高的研究,身高较高的父亲其儿子的平均身高要比父亲矮些,身高较矮的父亲其儿子的平均身高要比父亲高些,用高尔顿的话说就是”回归到平常”。虽然现在统计学上的”回归”这一概念已经远远超出的当时的定义,但是回归的原始思想依然有着非常重要的作用。”回归”,个人认为其实就是向中心的回归。在知道某地区18-24岁男子的身高的大致情况时,如果没有其他信息,让我们估计该地区中某一特定区域18-24岁男子的平均身高时(当然不是侏儒或篮球运动员之类的人),自然是用平均数(包括中位数)去估计了,这便是回归,没有其它的辅助信息时我们总倾向于平均值,这当然是符合统计思想的。两个变量的相关系数绝对值为1时,那么知道一变量的值就立即知道了另一变量的值;相关系数为0时,那么知道一变量的值对预测另一变量没有任何意义,那么我们就估计其值为平均值;相关系数绝对值介于0与1之间时,相关程度越大,我们越不倾向于取平均值。其实回归模型也是基于平均意义的,让我们来看看回归的本质(暂以两个变量x和y为例),回归是对每一个x值的y的平均值的估计,所以用回归模型来预测或估计总是平均意义的(这也是回归的思想),而针对某个特别的个体的预测则就需要非常的慎重了。

有这样一个例子,某学前班在儿童入学和结业时均要做智商测验,结果发现前后两次测验的分数平均都接近于100分,标准差为15分。但是仔细观察发现入学分数低于平均值的儿童结业时分数平均提高了5分,相反入学分数高于平均值的儿童结业时分数平均降低了5分,难道学前班会使儿童的智商平均化?其实没那么夸张,这只是回归效应的一个表现,只要两次测验分数的散点图中所有点不在同一条直线(这条直线的斜率为1)上,那么就会存在回归效应。观察得到的数据并不是真实值,都有或大或小的、或正或负误差,在大多数对称的概率分布中,观察值大于平均值的往往是其真实值加上了一个正的机会误差,观察值小于平均值的往往是其真实值加上了一个负的机会误差。所以在那个学前班中,入学分数较平均分低的儿童其真实分数一般是大于观察值的,因此在结业时的分数一般是要比入学时高,因为在向观察值的平均值,即真实值回归。

相关与回归是一定范围内的相关与回归,超出范围没有任何意义(经常实践的人应该会很少犯此类毛病的吧)。回归其实并不能增加信息量,它是一种结论(结论的准确性还有待评价),或对数据以某一种方式的总结,超出范围的估计预测是没有任何意义的。收入与教育水平有关,无休止的教育显然不会带来收入的持续的增加,所以人为地改变一个变量,希望通过回归模型的”魔力”来改变另一个变量是很荒谬的。另外,变量也是有范围或区域限制的,因此在使用回归模型做预测时是要非常谨慎的。

现在研究的回归往往都是多元回归,往往比较复杂,其实这是符合实际情况的,因此往往要用多个变量作为因子来拟合,但是这些变量是不是考察某一方面的较好指标呢,比如收入与教育水平有关,还可能与父母的社会地位有关,那这个”父母的社会地位”这一因子又该如何度量呢?这又是一个问题,尽管多元回归是一种非常有用的技术,但是永远代替不了对数据间内在关系的了解。由此可见实践经验的重要性!