所有由邱怡轩发布的文章

关于邱怡轩

中国人民大学统计学院硕士,普渡(众生)大学博士研究僧

共轭梯度法计算回归

共轭梯度示意图(图片来源:维基百科)
轮回眼 共轭梯度示意图(图片来源:维基百科

引子

之所以写这篇文章,是因为前几天统计之都的微信群里有同学提了一个问题,想要对一个很大的数据集做回归。然后大家纷纷给出了自己的建议,而我觉得共轭梯度算回归的方法跟这个背景比较契合,所以就正好写成一篇小文,与大家分享一下。

说到算回归,或许大家都会觉得这个问题太过简单了,如果用 $X$ 表示自变量矩阵,$y$ 表示因变量向量,那么回归系数的最小二乘解就是 $\hat{\beta}=(X’X)^{-1}X’y$。(本文完)



哎等等,别真走啊,我们的主角共轭梯度还没出场呢。前面的这个算系数的公式确实非常简洁、优雅、纯天然、不做作,但要往里面深究的话,还是有很多问题值得挖掘的。

最简单暴力的方法,就是从左向右,依次计算矩阵乘法,矩阵求逆,又一个矩阵乘法,最后是矩阵和向量的乘法。如果你就是这么算的,那么可以先默默地去面壁两分钟了。

更合理的方法,要么是对 $X’X$ 进行 Cholesky 分解,要么是对 $X$ 进行 QR 分解,它们基本上是现在算回归的软件中最常见的方法。关于暴力方法和矩阵分解方法的介绍和对比,可以参见这个B站上的视频。(什么?你问我这么严肃的话题为什么要放B站上?因为大部分时间都是在吐槽啊)

好,刚才去面壁的同学现在应该已经回来了,我们继续。前面这些通过矩阵运算求回归系数的方法,我们可以统称为直接法。叫这个名字,是因为它们都可以在确定数目的步骤内得到最终的结果。而与之相对的,则叫做迭代法,意思是通过不断更新已经得到的结果,来逐渐逼近真实的取值。打个比方,你想要知道一瓶82年的拉菲值多少钱,直接法就是去做调研,原料值多少,品牌值多少,加工费多少,运输费多少……然后加总起来得到最终的定价;而迭代法就是去问酒庄老板,你先随便蒙一个数,然后老板告诉你高了还是低了,反复循环,总能猜个八九不离十。

说到这里,你自然要问了,既然算回归的软件大都是用直接法,为什么还要考虑迭代法?莫非直接法有什么不好的地方?这就说到问题的点子上了。

继续阅读共轭梯度法计算回归

标题党统计学

如果你是被这个标题骗进来的,那么说明标题党的存在的确是有原因的。在网络高度发达(以及“大数据”泛滥)的今天,数据动不动就是以 GB 和 TB 的级别存储,然而相比之下,人类接受信息的速度却慢得可怕(参见大刘《乡村教师》)。试想一下,你一分钟能阅读多少文字?一千?五千?总之是在 KB 的量级。所以可以说,人们对文字的“下载速度”基本上就是 1~10KB/min。如果拿这个速度去上网的话你还能忍?

既然如此,每天网上有成千上万的新闻、报告、文章和八卦,怎么看得过来呢?没办法,只能先对正文进行一次粗略的筛选——看标题。俗话说得好,这是一个看脸的世界。于是乎,文章的作者为了吸引读者,就要取个足够博眼球的标题,而所谓标题党便是充分利用这种心理,用各种颇具创意的标题来吸引读者的注意。

好了,既然看官已经看到了这里,我就可以承认本文其实也是标题党了。这篇小文并不是要讨论标题党的前世今生,而是研究一个与此有关的统计问题:怎样的标题会更加吸引读者的关注?

这个问题有点太大了,所以我们缩小一下范围。既然是统计问题,就拿自家的一个例子下手吧:做统计学研究的,都得读各种各样的统计论文,那么论文的标题是否会对这篇文章的阅读量产生影响呢?巧的是,美国统计协会期刊(JASA)的网站上正好提供了该期刊旗下文章的下载访问量,所以我们可以以此做一个小分析,来研究一下标题与文章阅读量之间的关系。

可能有读者要问,为什么要使用文章的访问量,而不是引用率呢?这是因为 JASA 在其网站上说明,访问量数值是指从 JASA 官网下载的统计量,不包括从其他途径(比如购买的论文数据库)的来源。在 JASA 网站上,下载文章之前读者能获取到的主要是文章的标题和作者信息,所以访问量的主要驱动因素就是读者在阅读标题和作者之后产生的好奇感,从而减少了数据中的噪音。相反,引用一篇文章,通常是对文章有了充分理解之后产生的行为,这时候标题的作用可能就非常微弱了。总而言之,JASA 文章的下载量可以较好地代表读者在获取了文章的基本信息后对它感兴趣的程度。

jasa

继续阅读标题党统计学

美国统计协会开始正式吐槽(错用)P值啦

(图片来源:https://xkcd.com/1478,一幅讽刺滥用P值的漫画)

今天美国统计协会(ASA)正式发布了一条关于P值的声(吐)明(槽),算起来可以说是近期统计学界的一件大事了。为什么这么说呢?首先,P值的应用太广,所以对P值进行一些解释和声明非常有必要。其次,对P值的吐槽历来有之,但今天是第一次被一个大型的专业协会以非常正式的形式进行澄清,多少带有一些官方的意思。声明的全文可以在这个页面中下载。

那么这则声明里面都说了什么呢?小编整体读了一遍,把我认为重要的信息概括在这篇文章之中。

首先,ASA介绍了一下这则声明诞生的背景。2014年,ASA论坛上出现了一段如下的讨论:

问:为什么那么多学校都在教 p = 0.05?

答:因为那是科学团体和期刊编辑仍然在用的标准。

问:为什么那么多人还在用 p = 0.05?

答:因为学校里还在这么教。

看上去多少有点讽刺的味道,但事实却也摆在眼前。从舆论上看,许许多多的文章都在讨论P值的弊端,小编摘录了几条言辞比较激烈的:

这是科学中最肮脏的秘密:使用统计假设检验的“科学方法”建立在一个脆弱的基础之上。——ScienceNews(Siegfried, 2010)

假设检验中用到的统计方法……比Facebook隐私条款的缺陷还多。——ScienceNews(Siegfried, 2014)

针对这些对P值的批评,ASA于是决定起草一份声明,一方面是对这些批评和讨论作一个回应,另一方面是唤起大家对科学结论可重复性问题的重视,力图改变长久以来一些已经过时的关于统计推断的科学实践。经过长时间众多统计学家的研讨和整理,这篇声明今天终于出现在了我们面前。

P值是什么

这份声明首先给出了P值一般的解释:P值指的是在一个特定的统计模型下,数据的某个汇总指标(例如两样本的均值之差)等于观测值或比观测值更为极端的概率。

这段描述是我们通常能从教科书中找到的P值定义,但在实际问题中,它却经常要么被神话,要么被妖魔化。鉴于此,声明中提出了六条关于P值的准则,作为ASA对P值的“官方”态度。这六条准则算是这条声明中最重要的部分了。

继续阅读美国统计协会开始正式吐槽(错用)P值啦

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

标准正态分布的分布函数 $\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

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

极简 Spark 入门笔记——安装和第一个回归程序

现在的各种数据处理技术更新换代太快,新的名词和工具层出不穷,像是 Hadoop 和 Spark 这些,最近几年着实火了一把。事实上听说 Spark 也有一段时间了,但一直是只闻其名不见其实,今天就来简单记录一下初学 Spark 的若干点滴。

Spark 是什么

按照 Spark 官方的说法,Spark 是一个快速的集群运算平台,以及一系列处理大型数据集的工具包。用通俗的话说,Spark 与 R 一样是一套用于数据处理的软件和平台,但它最显著的特点就是处理大型数据(我就是不说大数据 ( ̄^ ̄))的能力。

极简安装

Spark 本身面向的是大规模的分布式计算,但对学习和测试来说,利用单机的多核 CPU 就已经足够了,所以作为入门,我并没有打算去涉及多台计算机相连的情形。在这个基础上,第一件出乎我意料的事情就是,Spark 的安装和配置其实可以是异常简单的。

在网上出现的各种资料中,Spark 经常与 Hadoop 和 Scala 这两个名词一起出现。前者也是一个大型分布式计算的框架,诞生得比 Spark 更早;后者是 Spark 主要使用的一种编程语言。这就给不明真相的群众造成了一种印象,好像要使用 Spark 的话就得先安装配置好 Hadoop 和 Scala,而要安装它们又得有更多的软件依赖。但实际上,要在单机上使用 Spark,真正需要的只有下面几样:

  1. 一台金光闪闪的电脑
  2. 在上面这台电脑里面装一个金光闪闪的 Linux 操作系统
  3. 在上面这个系统里面装一个金光闪闪的 Java 开发环境(JDK)

这三样可以说是大部分计算环境的标配,如果系统还没有安装 JDK,那么一般都可以用系统的包管理工具,比如 Fedora 下是

sudo yum install java-1.8.0-openjdk

Ubuntu 下是

sudo apt-get install openjdk-7-jdk

有了上面的开发环境,安装 Spark 就非常容易了,基本上只要下载预编译包,解压缩,然后添加系统路径即可。首先,到 https://spark.apache.org/downloads.html 选择最新的 Spark 版本和 Hadoop 版本(实际上我们暂时用不上 Hadoop,所以任何版本都行),然后下载压缩包。 继续阅读极简 Spark 入门笔记——安装和第一个回归程序