分类目录归档:优化与模拟

传统统计计算,包括优化方法和统计模拟

Persi Diaconis (1)

作为统计之美的开篇,我一直想找一篇我非常愿意写的统计故事,尽管有很多,但都不能让我觉得可以发泄笔头之愤。最近在听贝叶斯统计课,刘军老师(哈佛大学统计系教授)提起了叫Persi Diaconis的人,他的故事和他的工作,这让我找到了写这篇文章的灵感。

你能想象,一个人在14岁离家出走,学习魔术,浪迹江湖,24岁后潜心学术,之后成为斯坦福大学的教授?

Persi Diaconis(维基)确实如此,他在搞魔术的时候,为了想研究如何防止被其他魔术师骗,买了本 William Feller 的 An Introduction to Probability Theory and Its Applications,但是里面涉及到了微积分等知识,看不懂,那年他18岁。他发誓要回学校学习,以此能够看得懂这本书。24岁重返校园(City College of New York)。他向《科学美国人》投稿介绍他两个有意思的洗牌方法。这样被一个马丁·葛登能的人看重,给他写了推荐信去哈佛大学,当时哈佛的统计学家 Fred Mosteller 正在研究魔术,于是就要了他(http://blog.sciencenet.cn/home.php?mod=space&uid=1557&do=blog&id=418859)。

Persi Diaconis 做了几个很有意思的工作,如洗牌多少次能够洗得比较彻底(我希望在统计之美里面,能够有一篇来单独介绍洗牌问题)等。他还有个绝活,据刘军老师说,他总可以抛硬币时,抛出他想要的那一面。而他每次的学术报告之前,都会表演一番,很多人实际上不是去听他的报告的,而是看看他的绝活。

这次就要介绍的,是他在一篇 MCMC 算法(Markov Chain Monte Carlo,马氏链蒙特卡洛方法,有文章将其评为20世纪最有名的10个算法之一)综述的文章(The Markov Chain Monte Carlo Revolution)中,给出的破译犯人密码的例子。

有一天,一个来自了解关押囚犯心理的心理医生,来到斯坦福统计系,给出了如下一个囚犯写的密码信息:

上图是囚犯写的密码信息的一部分,你可以看到很多出现频繁的字符。

问题来了,该心理医生想知道,这个密码信息的内容是什么?

我们可以想到,上图看起来怪怪的字符,每个都应该对应一个字母,只要我们找到字符和字母的对应关系,我们就可以解码了。但是怎么找到这个对应关系呢?(我想到了福尔摩斯探案集里面有这样一个例子,不过那个例子中的字符代表体系和这个不同,但是福尔摩斯的推断相当惊人!)

在下一篇我给出他的想法,从而写完MCMC算法的引子。

win your ex girlfriend backfree advice on how to make her want you back How To Get Your Ex Girlfriend Back how to get back with your ex girlfriendhow to get your girlfriend back

MCMC案例学习

本文是R中mcmc包的一篇帮助文档,作者为Charles J.Geyer。经过knitr编译后的pdf文档可见此处,提供中文译稿的作者:
闫超,天津财经大学统计系2011级研究生,方向:非寿险准备金评估。
高磊,天津财经大学统计系2011级研究生,方向:非寿险准备金评估。

这个案例,我们不关心题目的具体意义,重点放在利用贝叶斯的观点来解决问题时,MCMC在后续的计算中所发挥的巨大作用。我们知道,贝叶斯的结果往往是一个后验分布。这个后验分布往往很复杂,我们难以用经典的方法求解其期望与方差等一系列的数据特征,这时MCMC来了,将这一系列问题通过模拟来解决。从这个意义上说,MCMC是一种计算手段。依频率学派看来,题目利用广义线性模型可以解决,在贝叶斯看来同样以解决,但是遇到了一个问题,就是我们得到的非标准后验分布很复杂。我们正是利用MCMC来解决了这个分布的处理问题。本文的重点也在于此。

在使用MCMC时作者遵循了这样的思路,首先依照贝叶斯解决问题的套路,构建了非标准后验分布函数。然后初步运行MCMC,确定合适的scale。继而,确定适当的模拟批次和每批长度(以克服模拟取样的相关性)。最后,估计参数并利用delta方法估计标准误。

1. 问题的提出

这是一个关于R软件中mcmc包的应用案例。问题出自明尼苏达大学统计系博士入学考试试题。这个问题所需要的数据存放在logit数据集中。在这个数据集中有五个变量,其中四个自变量x1、x2、x3、x4,一个响应变量y

对于这个问题,频率学派的处理方法是利用广义线性模型进行参数估计,下面是相应的R代码以及结果:

library(mcmc)
data(logit)
out <- glm(y ~ x1 + x2 + x3 + x4, data = logit, family = binomial(), x = T)
summary(out)

Call:
glm(formula = y ~ x1 + x2 + x3 + x4, family = binomial(), data = logit, 
    x = T)

Deviance Residuals: 
   Min      1Q  Median      3Q     Max  
-1.746  -0.691   0.154   0.704   2.194  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)   
(Intercept)    0.633      0.301    2.10   0.0354 * 
x1             0.739      0.362    2.04   0.0410 * 
x2             1.114      0.363    3.07   0.0021 **
x3             0.478      0.354    1.35   0.1766   
x4             0.694      0.399    1.74   0.0817 . 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 137.628  on 99  degrees of freedom
Residual deviance:  87.668  on 95  degrees of freedom
AIC: 97.67

Number of Fisher Scoring iterations: 6

但是,使用频率学派的分析方法解决这个问题并不是我们想要的,我们希望使用Bayesian分析方法。对于Bayesian分析而言,我们假定数据模型(即广义线性模型)与频率学派一致。同时假定,五个参数(回归系数)相互独立,并服从均值为0、标准差为2的先验正态分布。

定义下面的R函数来计算非标准的对数后验分布概率密度(先验密度与似然函数相乘)。我们为什么要定义成密度函数的对数形式?因为虽然我们是从分布中采样,但是MCMC算法的执行函数metrop()需要的一个参数正是这个分布的密度函数的对数形式。

x <- out$x
y <- out$y
lupost <- function(beta, x, y) {
    eta <- as.numeric(x %*% beta)
    logp <- ifelse(eta < 0, eta - log1p(exp(eta)), -log1p(exp(-eta)))
    logq <- ifelse(eta < 0, -log1p(exp(eta)), -eta - log1p(exp(-eta)))
    logl <- sum(logp[y == 1]) + sum(logq[y == 0])
    return(logl - sum(beta^2)/8)
}

继续阅读MCMC案例学习

议员是如何投票的?

[latexpage]

一、议员投票

这个数据在近几年的图模型文章中常能见到,并且已有很多深入的讨论——包括图结构随时间变化、多图联合估计等情况。本文只涉及单个图结构的估计,此外笔者对政治不了解,因此文中摘录wiki的相关评论。

从 http://www.senate.gov 可以看到senators每次投票的结果。那么,你关心的议员在每次投票中起到怎样的作用?有怎样的政治立场?本文数据选取美国第111届议会(2009年、2010年)roll-call vote数据,涉及110位参议员696个议案的投票结果(点击此处下载数据)。利用图模型,我们可以对议员投票行为作图。下图中,每个点代表一个议员,每条边表示两个议员的投票行为“很相似”——他们常常同时投赞成票或反对票。

图中,两派阵营的状态很明显。绿色的是Democrat,红色的是Republican,蓝色的是Independent。一般认为,民主党在政治上偏左,主张社会自由与进步;而共和党偏右,“reflects American Conservatism”。可以从图中获得很多信息,下面我们对此图做一点深究。

abc

首先,红绿交界处的Lincoln、Nelson、Collins等人是否代表了某种“中间力量”?以下从维基摘录的一段话(http://en.wikipedia.org/wiki/Bill_Nelson)或许回答这个问题。

Nelson’s votes have tended to be more liberal than conservative. He has received high ratings from left-of-center groups such as Americans for Democratic Action, and low ones from right-of-center groups such as the Eagle Forum and the Club for Growth. According to ratings by the National Journal, Nelson’s votes have been liberal on economic matters, moderate on social issues, and liberal but close to the center on foreign policy.

其次,我们再来看另一个极端。选取图中远离红点的一个绿点(图中左下角)——Cantwell;以及远离绿色区域的红点(图中右下角)——Enzi。同样是摘自wiki的两段话,验证了他们的“极左”和“极右”。

While she(Cantwell) scores high on a progressive chart from ProgressivePunch.org,[25] Cantwell has made several controversial votes during her time in the Senate that have created friction between her and members of the Democratic Party.

Enzi was ranked by National Journal as the sixth-most conservative U.S. Senator in its March 2007 conservative/liberal rankings.[7] Despite his strong support of the War in Iraq, he was one of 14 U.S. Senators to vote against the Iraq War funding bill in May 2007 because he opposes the clauses of the bill which increase domestic spending.

此外,被绿点包围的蓝点Lieberman——“A former member of the Democratic Party, he was the party’s nominee for Vice President in the 2000 election. Currently an independent, he remains closely affiliated with the party.”

除了以上分析,还可以获得很多信息。比如可以分析图中的重要节点、用局部估计的方法了解投票关联是如何随时间变化的,等等。

二、图模型——罚极大似然

如何才能获得这幅“议员投票图”呢?它的理论依据是什么?估计图结构的方法也有很多,这里只介绍罚极大似然的方法。

图结构估计得常用假设是数据来自正态分布——注意上例的数据为二项分布,banerjee(2007)在论文中探讨了二项分布数据的问题,最终所求解的最优化表达式和正态类似;本文中暂且忽略banerjee(2007)中的细节,把投票数据用正态分布处理。正态分布的概率密度函数为
\[
f(x)=\frac{1}{(2\pi)^{p/2}|\Sigma|^{1/2}}\exp\Big\{-\frac{1}{2}(X-\mu)^T\Sigma^{-1}(X-\mu)\Big\}
\]
假设数据$x_1, x_2, \cdots, x_n$来自正态分布,每个$x_i$(p维向量)记录p位议员的一次投票结果。如果用协方差矩阵描述议员投票的相关关系,那么它的极大似然估计为,
\[
S=\frac{1}{n}\sum_{i=1}^N(x_i-\bar{x})(x_i-\bar{x})^T
\]
协方差矩阵的逆$\Omega=\Sigma^{-1}=\big(\omega_{ij}\big)_{p\times p}$在图模型中被称为Concentration Matrix或者Precision Matrix。该矩阵与偏相关系数有如下关系:
\[
\rho_{ij|\\\{i,j\}}=-\frac{\omega_{ij}}{\sqrt{\omega_{ii}\omega_{jj}}}
\]

由此可以看出,$\Omega$矩阵中的零元素代表了对应议员投票行为的条件独立关系。

为了在高维问题中更精确地估计$\Omega$矩阵,可以采用罚极大似然估计。首先回顾对数似然函数的表达式
\[
\begin{split}
\sum_{i=1}^n\log & f(x_i)=\frac{np}{2}\log(2\pi)-\frac{n}{2}\log|\Sigma|
-\frac{1}{2}\sum_{i=1}^n(x_i-\bar{x})^T\Sigma^{-1}(x_i-\bar{x})\\
&=\frac{np}{2}\log(2\pi)-\frac{n}{2}\log|\Sigma|
-\frac{1}{2}tr\Big(\sum_{i=1}^n(x_i-\bar{x})(x_i-\bar{x})^T\Sigma^{-1}\Big)\\
&\propto p\log(2\pi)-\log|\Sigma|-tr\big(S\Sigma^{-1}\big)
\end{split}
\]
此外,由于$\Omega=\Sigma^{-1}$,令上式最大等价于求解下式的最小值
\[
tr(S\Omega)-\log|\Omega|
\]
在$\Omega$矩阵非零元素稀疏的假定下,通过加罚,我们可以得到以下的最优化问题(解空间为正定矩阵):
\begin{equation}
\min_{\Omega\succ 0}\Big\{tr(S\Omega)-\log|\Omega|+\lambda\sum_{i\neq j}|\omega_{ij}|\Big\} \label{obj:main}
\end{equation}

针对此最优化问题,Yuan and Lin(2007)通过maxdet最优化算法求解,banerjee(2007)则借用内点搜索方法(interior point method),此外Friedman(2007)则将它转化为Lasso 的表达形式再通过shooting 算法求解。还有一类算法是从贝叶斯角度考虑,Bayesian Graphical Lasso,思想很像 Bayesian Lasso.

 求解(1)式,得到$\Omega$矩阵,它的非零元素便对应“投票图”中的一条边。

 

三、数据与代码

R数据下载——数据中包含三个矩阵(数据来自http://www.senate.gov)。

  • idList矩阵为参议员信息。每行对应一个参议员。5列分别对应:议员id,last name, first name, state, party(Democrat, Republican or Independent)
  • bilInfo为投票议案信息(bil info)。每行对应一项议案信息,4列分别为:议案编号、投票时间、议案描述、投票结果。
  • voteMatrix矩阵为投票信息,每行对应bilInfo的议案,每列对于idList的一个议员。若voteMatrix的第i行第j列数值为1,表示议员j对议案i投了赞成票;-1表示反对票;NA表示未投票。

代码很简单,只有以下几行:

library(spaceExt)
library(igraph)
load("senate.RData")
#移除投票缺失较多的议员
sel=which(!(colSums(is.na(voteMatrix))>100))
partyD=as.numeric(idList[sel,5]=="D")
partyI=as.numeric(idList[sel,5]=="ID" | idList[sel,5]=="I")
#用spaceExt做计算
res=glasso.miss(voteMatrix[,sel],20,rho=0.1,penalize.diagonal=FALSE)
#计算偏相关系数矩阵
p=-res\$wi
d=1/sqrt(diag(res\$wi))
p=diag(d)%*%p%*%diag(d)
diag(p)=0
#igraph包生成图模型、作图
g=graph.adjacency((p>(0.055)),mode="undirected",diag=F)
V(g)\$color=(partyD+2+partyI*2)
V(g)\$label=idList[sel,3]
par(mar=c(0,0,0,0))
plot(g,layout=layout.fruchterman.reingold, edge.color=grey(0.5),vertex.size=10)

相关文章补充:

  1. Yuan, M., and Lin, Y. (2007) Model selection and estimation in the Gaussian graphical model. Biometrika , 94, 1, pp. 19–35
  2. Friedman, J. H., Hastie T, Tibshirani R. “Sparse inverse covariance estimation with the graphical lasso.” Biostat (2008) 9 (3): 432-441.
  3. Banerjee, O., Ghaoui, L. E. and d’Aspremont, A. (2007), Model selection through sparse maximum likelihood Estimation, J. Machine Learning Research 101.
  4. Wang, H. The Bayesian Graphical Lasso
  5. Mladen Kolar, Le Song, Amr Ahmed, Eric P. Xing. Estimating time-varying networks. Annals of Applied Stat. arXiv:0812.5087v2

Sweave后传:统计报告中的大规模计算与缓存

学无止境。我曾以为我明白了如何在Sweave中使用缓存加快计算和图形,但后来发现我并没有真的理解,直到读了另外一些手册才明白,因此本文作为前文“Sweave:打造一个可重复的统计研究流程”之续集,向大家介绍一下如何在Sweave的计算和图形中使用缓存,以节省不必要的重复计算和作图,让那些涉及到密集型计算的用户不再对Sweave感到难堪。

如果你还没读前文,建议先从那里开始读,了解Sweave与“可重复的统计研究”的意义。简言之,Sweave是一种从代码(R代码和LaTeX)一步生成报告的工具,我们可以把整个统计分析流程融入这个工具,让我们的报告具有可重复性。然而,就普通的Sweave而言,这样做的一个明显问题就是,所有计算和作图都被融入一个文档之后,每次运行这个文档都要重复所有的计算和作图,这在很多情况下纯粹是浪费时间;比如,我只想对新添加的部分内容运行计算,而文档中的旧内容希望保持不变。这都是很合理的需求,我们需要的实际上就是一种缓存机制,将不想重复计算的对象缓存起来,需要它的时候再从缓存库中直接调出来用。

Sweave本身做不到缓存(或者很难做到),但由于R语言的扩展性以及一些疯子和天才的存在,我们可以通过R的一些附加包来完成缓存。这里的缓存有两种(不关心细节的读者可以略过):

  • 一种是对R对象的缓存,比如数据、模型(拟合好的模型)和其它计算结果等,这个缓存容易理解,我们可以把它简单想象为把计算过程中的对象存入某些特殊的缓存库(文件),下次计算的时候先查看缓存库中是否存在我们需要的对象,如果存在,那么就直接加载进来,否则就重新计算;如(伪代码):
    ## 如果存在,就加载
    if (file.exists('x.RData')) {
        load('x.RData')
    }
    ## 如果当前工作环境中不存在x,那么就重新创建它
    if (!exists('x')) {
        x = rnorm(100000)
    }

    R对象的缓存可以通过cacheSweave包(作者Roger D. Peng)实现,它真正采用的方法比上面的伪代码当然要高级(比如用代码的MD5值作数据库名等),不过用户可以不必关心这其中的细节。这种“当需要时才加载”的方式在术语上叫“延迟加载”,即Lazy Load,这在R里面也很常见,比如很多R包的数据都是采用延迟加载的方式(加载包的时候并不立刻加载数据,而是用data(name)在需要的时候加载)。

  • 另一种是图形的缓存,这是一个不可思议的魔法。cacheSweave的帮助文档和vignette中特别提到了它只能缓存R对象,无法缓存图形和其它附属输出,而开源世界的魅力就是你总能找到一些奇妙的解决方案。R包pgfSweave借力于pgf实现了图形的缓存。pgfSweave包的两位作者Cameron Bracken和Charlie Sharpsteen都是水文学相关专业出身,却在R的世界里做了这些奇妙的工作,让我觉得有点意外,这是题外话。pgf是LaTeX世界的又一项重大发明,它提供了一套全新的绘图方式(语言),而且它的图形可以通过LaTeX直接生成PDF输出(这是pgfSweave能实现缓存的关键)。关于pgf,我们可以翻一翻它726页的手册,如果你能坚持看上10分钟,感叹超过50次,那么你一定是一位超级排版爱好者。它的图形质量之高、输出之漂亮,真的是让人(至少让我)叹为观止。同样,这里我也不详细介绍细节。pgf图形的“缓存”是通过“外置化”(externalization)来实现的,简言之,pgf图形有一种输出形式如下:
    \beginpgfgraphicnamed{graph-output-1}
    \input{graph-output-1.tikz}
    \endpgfgraphicnamed
    

    如果LaTeX在运行过程中遇到这样的代码段,那么它会把tikz图形(tikz是高层实现,pgf是基础设施)先编译为PDF图形,等下次重复运行LaTeX的时候,这段代码就不必重新运行了,LaTeX会跳过它直接引入PDF图形。这样就省去了编译tikz图形的时间,同时也能得到高质量的图形输出。

    pgfSweave从R的层面上可以输出这样的代码段,并生成一个shell脚本来编译图形为PDF;而pgfSweave也会尝试根据R代码的MD5来判断一段画图的代码是否需要重复运行——如果R代码没有变化,那么图形就不必重制。这是自然而然的“缓存”。对了,pgfSweave对Sweave的图形输出作了扩展,新增了tikz输出,这是通过另一个R包tikzDevice实现的(同样两个作者),这个包可以将R图形录制为tikz绘图代码。pgfSweave依赖于cacheSweave包,把它缓存R对象的机制也引入进来,所以我们可以不必管cacheSweave,把精力放在pgfSweave上就可以了。

介绍完血腥的细节之后,我们就可以坐享高科技的便利了。下面我以一个统计计算中的经典算法“Gibbs抽样”为例来展示pgfSweave在密集型计算中的缓存功能。

Gibbs抽样的最终目的是从一个多维分布$f(X_1,X_2,\ldots,X_n)$中生成随机数,而已知的是所有的条件分布$f(X_1|X_2,\ldots,X_n)$、$f(X_2|X_1,\ldots,X_n)$、……、$f(X_n|X_1,\ldots,X_{n-1})$,并且假设我们可以很方便从这些条件分布中抽样。剩下的事情就很简单了:任意给定初始值$x_1^{(0)},x_2^{(0)},\ldots,x_n^{(0)}$,括号中的数字表示迭代次数,接下来依次迭代,随机生成$x_1^{(k)} \sim f(X_1|x_2^{(k-1)},\ldots,X_n^{(k-1)})$、$x_2^{(k)} \sim f(X_2|x_1^{(k)},\ldots,x_n^{(k-1)})$、……、$x_n^{(k)} \sim f(X_n|x_1^{(k)},\ldots,x_{n-1}^{(k)})$(迭代次数$k=1,\cdots,N$),也就是,根据上一次的迭代结果并已知条件分布的情况下,从条件分布中生成下一次的随机数(这个算法很简单,你什么时候看明白上标下标了就明白了)。这是个串行的过程,所以无法避免循环。

接下来的计算没有那么复杂,我们仅仅用一个二维正态分布为例。当我们知道一个二维正态分布的分布函数时,我们也知道它的两个边际分布都是一维正态分布,而且分布的表达式也很明确,后文我会详细提到,读者也可以参考维基百科页面。生成一维正态分布的随机数对我们来说当然是小菜一碟——直接用rnorm()就可以了。这样,我们只需要交替从$f(X|Y)$和$f(Y|X)$生成随机数,最终我们就能得到服从二维联合分布的随机数(对)。假设我们要从下面的二维正态分布生成随机数:

$
\left[\begin{array}{c}X\\ Y\end{array}\right]\sim\mathcal{N}\left(\left[\begin{array}{c}0\\ 1\end{array}\right],\left[\begin{array}{cc} 4 & 4.2\\ 4.2 & 9\end{array}\right]\right)
$

详细过程和结果参见下面这份PDF文档(点击下载):

A Simple Demo on Caching R Objects and Graphics with pgfSweave (PDF)

二维正态分布随机数及其等高线图
二维正态分布随机数及其等高线图

我们生成了50万行随机数,并画了X与Y的散点图。由于我们设定了相关系数为0.7,所以图中自然而然显现出正相关;而等高线也体现出多维正态分布的“椭球形”特征。均值在(0, 1)附近,都和理论分布吻合。所以这个Gibbs抽样还不太糟糕。

生成上面的PDF文档和图形的LyX/Sweave源文档在这里:

A Simple Demo on Caching R Objects and Graphics with pgfSweave (LyX)

如果你已经按照我前面的文章配置好你的工具(即使当时配置过,现在也需要重新配置,因为我最近作了重大修改),这个文档应该可以让你重新生成我的结果。文档中有两处关键选项:

  • cache=TRUE 使得文档中的R对象可以被缓存;
  • external=TRUE 使得图形可以被缓存;

第一次运行文档时,也许需要等待一分钟左右的时间,因为R需要长时间的Gibbs抽样计算,并编译生成的图形为PDF,但如果下一次再运行这个文档,则速度会快很多,因为不需要再计算任何东西,只需从缓存加载就可以了。我们前一次的抽样结果已经在缓存库中了。

本文至此可以结束,但真正想享用这个结果的读者可能还要继续看几个问题:

  1. pgfSweave对图形的缓存经常有点过度,即使你改动了代码,缓存仍然不会被替换,这种情况下你可以把生成的tikz和pdf图形删掉;
  2. 中文问题:tikzDevice包好像还没有完全解决中文问题,作者正在尝试我提供的中文图形,也许在未来不久的新版本中会有完全的中文支持,但目前实际上我已经在我的自动配置脚本中解决了这个问题。我们只需要保证我们的中文LyX文档用UTF-8编码就可以了——在文档选项中使用CTeX类(以及UTF8选项),并设置语言为中文、编码为Unicode (XeLaTeX) (utf8)。这种设置至少我在Windows和MikTeX下可以成功编译中文Sweave文档。目前高亮代码功能在中文文档中还不能用,但我已经向相关作者发了补丁,等CRAN管理员回来发布新版parser包之后这个问题也会解决。所以,中文用户也不必担心,一切都(将)可以使用。
  3. 依旧是样式问题:细心的读者也许能发现上面的PDF文档中的图形和这里网页中的图形略有区别,PDF文档中图形的字体和正文字体是一致的,这也是pgf的优势,前文第5节曾提到过。
  4. 如果使用缓存,那么代码段中的对象都无法打印出来,即使用print(dat)也无法打印dat,这是因为缓存的代码段只管是否加载对象,而无法执行进一步的操作(包括打印),对于这种情况,我们应该设计好我们的Sweave文档结构,让某些代码段作计算(并缓存),某些代码段作输出(不要缓存)。

这两篇文章都有点长,但最终的操作都极其简单。如果能理解这套工具,我想无论是写作业还是写报告,都将事半功倍(比如我就一直用它写作业)。

祝大家在新的一年能写出漂亮、利索的统计文档,告别复制粘贴的体力劳动。

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

昨日翻看朱世武老师的《金融计算与建模》幻灯片(来源,幻灯片“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语言诞生的原因之一:统计学的编程不应该涉及到那么多繁琐的循环。