分类目录归档:统计计算

贝叶斯统计学与传统优化方法

共轭梯度法计算回归

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

引子

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

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



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

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

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

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

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

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

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

标准正态分布的分布函数 $\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个独立同分布的均匀分布即可,看代码

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

嘿,朋友,抢红包了吗?

如果你有一台智能手机,如果你装了一个名叫微信的软件,那么你今年的春节很可能是在下面这样的场景中度过的(图片来自微信群):

lucky_money_3 lucky_money_2 lucky_money_1

这也使得众多的网络大V发出了下面的感慨:

lucky_money_4

而最近几天不少微信群里面又流行起来一种“红包接力”的玩法,大概的规则是:群里面先由一人发一个红包,然后大家开始抢,其中金额最大的那个人继续发新一轮的红包,之后不断往复循环。

这时候大家或许就会问了,一直这么玩下去会有什么结果呢?是“闷声赚大钱”了,还是“错过几个亿”了?是最终实现“共同富裕”了,还是变成“寡头垄断”了?要回答这些问题,我们不妨用统计模拟的方法来做一些随机实验,得到的结果或许会让你大跌眼镜呢。

继续阅读嘿,朋友,抢红包了吗?

COS每周精选:算法学习知哪些?

本期投稿:谢益辉 王威廉   冷静   王小宁

编辑:王小宁

算法

K-means是最常用的聚类算法之一:容易理解,实现不难,虽然会有local optimum,但通常结果也不差。但k-means也不是万金油,比如在一些比较复杂的问题和非线性数据分布上,k-means也会失败。普林斯顿博士David Robinson写了一篇不错的分析文章,介绍了几种k-means会失效的情形。

基于遗传算法的小车模拟, 还有遗传算法的行者,看着有一大拨僵尸来袭的感觉、遗传算法的猫

继续阅读COS每周精选:算法学习知哪些?