解名缰 鸟倦飞

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

2016 年 01 月 30 日 | 分类于 学习中

标准正态分布的分布函数 $\Phi(x)$ 可以说是统计计算中非常重要的一个函数,基本上有正态分布的地方都或多或少会用上它。在一些特定的问题中,我们需要大量多次地计算这个函数的取值,比如我经常需要算正态分布与另一个随机变量之和的分布,这时候就需要用到数值积分,而被积函数就包含 $\Phi(x)$。如果 $Z\sim N(0,1), X\sim f(x)$,$f$ 是 $X$ 的密度函数,那么 $Z+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 <- pnorm(x))
## user  system elapsed
## 1.049   0.000   1.048
system.time(asy <- r_as26217ncdf(x))
## user  system elapsed
## 0.293   0.019   0.311
max(abs(y - asy))
## [1] 6.968772e-08

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

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

阅读全文→

周记(四)

2015 年 12 月 19 日 | 分类于 生活学习中

上周的周记缺了一次,原因是期末最后的这两周基本上是在各种死线中度过的,于是没有精力也没有心情去写博客了。

其实上上周有好几条突如其来的消息。谢老大之前身体不适,前两周一直在医院检查,直到这周一做了一晚上手术,向师姐了解情况后,得知最后问题不是太大,应该能很快恢复。也是两周之前,瑶总滑雪摔伤,经过前前后后的诊断,需要一个来月的时间恢复静养,其中滋味必然不好受,希望伤痛能尽快退去,新年无忧无虞。还有一条挺沉重的消息,是两周前陈博士在微信群里告诉我们的。当年和我们一起入学的有位瑞典小哥,他在13年暑假的时候去南美旅游,结果被报失踪,再也没有回来。陈博士偶然在浏览相关报道的时候发现今年5月份有更新的消息,说是他的遗骸被找到了,与当地的黑帮有关。这种在电影里面出现的情节就发生在自己身边的人,一时还难以接受。

负能量过后,生活依然还在继续,忙碌了两周,把剩下的期末报告都写完之后,今天就算是进入假期了。

周二和导师见面,闲聊的时候不经意提到了主动性的问题,他说读书期间要学会主动一点,如果有机会就应该多把自己做过的东西展示出来,之前十月份的时候他让我参加泛华统计协会的Poster比赛就是这个原因。他说这种类似的比赛有很多,如果有一些成果都可以投出去试试,拿不拿奖是其次,关键是要学会怎么让别人更多地了解你。

在这一点上,我确实得反思一下。很多时候我会觉得宣传自己是一件没有必要而且有些难以启齿的事,所以很少会主动跟别人介绍自己的能力或是成果,为此身边的朋友都或多或少觉得我过于低调和谦卑了。任何事情推到极致都是有问题的,关键是如何达到一个平衡,导师和朋友的劝诫我都记下了。

说起这事是因为和另一个话题有关。之前谢老大跟我提起过好几回参加John Chambers奖比赛的事,而我一直觉得没有什么值得参赛的作品,所以并没有太上心。当天听了导师一席话后,觉得确实得有所行动,所以就决定拿暑假写的一个软件包参加比赛。虽然我对结果并不抱什么希望,但我发现给自己人为地加一个截止日期,不管有意无意,都能极大提高工作的效率。事实上比赛当天就要截止,于是我一整晚都在写各种文档,完善一些细节,最终在截止前一个小时把申请投了出去。

参加这个比赛给我带来的另一个收获是它促使我做了一些之前一直想做但就是没有动力去做的事,比如给写的这个软件做一个网站。因为申请材料的篇幅有限,所以很大程度上得依靠材料之外的媒介对参赛作品进行描述,而网站一般是首选。由于前一阵为了做R会议的网站搭过一个简单的框架,所以这次上手也很快,我对做出来的效果也比较满意。让我有些吃惊的是,谢老大在第二篇周记里留过言,说“做一个好看的网站,看起来不务正业,却是有实在产出并对将来能产生影响的事情,所以是看起来没用但重要的事情”,不想这么快预言就成真了。

忙完这许多事情,算是可以稍微休息一阵了。今年村里的冬天很反常,气象温和,阳光常在,下午坐在书桌前,发现窗台前的花草在卧室墙上投出了一道映像。素影参差,年光催度,希望亲朋好友都能在年底有个完满的收官,把不如意之事放在旧年,给马上到来的2016一个崭新的开端。

影

阅读全文→

周记(三)

2015 年 12 月 05 日 | 分类于 生活学习中

周日在朋友家吃火锅,其间Z同学喊了句“藕下了啊”,我便调侃了句“这话在好多年前的 BBS 上就是我要下线了的意思”,然后W同学评论说暴露年龄了,于是接下来又互相聊出了更多暴露年龄的典故,像是多图杀猫,九城社区,东北人都是活雷锋,小小,Showgood 之类。这么列举一下,发现2000年初的网络文化其中很重要的一支就是 Flash 动画,而这其中 Showgood 毫无疑问是当时中国 Flash 界的主力军之一。

第一次听说 Showgood 我记得是初中的时候去P姑娘家玩,她推荐了大话三国系列中的某一部作品给我看,当时我就被那部动画的画风和搞笑台词惊艳到了,后来家里能拨号上网,就自己去 Showgood 的网站上看,现在依然还有印象的就是三集左右的《千里走单骑》,还有后来的《小兵的故事》。下午回到公寓后特意又搜了下,发现 Showgood 这个网站居然还在,只是已经成了一家做互联网儿童教育的创业公司。尽管已经改头换面,但进网站后迎面那句“这些年,你们还好吗?”还是让人唏嘘不已。

但其实在我心目中真正的 Flash 神作还是小小系列,也就是当时风靡整个 Flash 界的火柴人动画。我至今还记得小小的域名是 xiaoxiaomovie.com,用 whois 查了下,并没有被注销掉,只是网站已经无法访问了。小小的影响力之大,我清楚地记得曾经玩过一个英文的 Flash 格斗游戏,里面可选的角色是各类美国的超级英雄和日本的动漫角色,但列表中赫然有着 Xiao Xiao 的名字。可惜后来再想找到这个游戏,就如同大海捞针了。少年时的记忆就像是擦肩而过的路人,偶尔留下一抹难忘的影子,又很快地融入到黑压压的人群中,再不见踪迹。

阅读全文→

周记(二)

2015 年 11 月 27 日 | 分类于 生活

上周周六,入冬以来的第一场雪,气温骤降,让人猝不及防。突然从清爽中还带着一丝明媚的“深秋”时光变到压抑晦暗的冬天,一时还有些不习惯。

不过冬天依然是个好季节,很多事情都是在冬天的时候想通的。

初雪

周日做了一次大胆的尝试,从早上八点到晚上八点,把手机一直保持在关机状态,看能不能存活下来。后来事实证明不但成功地活了下来,还大大提高了做事的效率,先是把该干的活干完了,还抽空读掉了一本小说,就是上周因为写周记而没有读的那本,2014年诺奖得主帕特里克•莫迪亚诺的《暗店街》。

书是室友的,当时看封面内页的介绍,有些侦探小说的感觉,于是好奇就借来看了。书很薄,不到200页,总的故事走向也很简单,就是在二战后的巴黎,主角因为遗失记忆而四处寻找线索,希望能揭开自己的身世之谜。然而在读完最后一章后发现根本就没有真相大白的一幕,直到最后,主角还在寻找、怀疑和挣扎,依靠一些微弱的线索去拼凑曾经的自己。从理性的角度讲,主角是可以彻底与过去割裂开的,然而寻找本源似乎是所有人与生俱来的一种天性,以至于作者借另一个角色之口说出了“在生活中重要的不是未来,而是过去”。我非常喜欢书里面的一个比喻,他说我们都是一种 海滩人,我们无数次地出现在海滩度假照片中的一角或背景中,出现在欢快的人群中,但谁也叫不出我们的名字,说也说不清我们为何在那儿。“沙子只把我们的脚印保留几秒钟”,当我们离开沙滩以后,与这个世界的联系就仅剩许许多多互不相关的碎片而已。

豆瓣上一条评论说这本书得在一个阳光明媚的午后拿来读,否则容易想多,我觉得是这么回事。

阅读全文→

周记

2015 年 11 月 20 日 | 分类于 生活

今天晚上突然有个疯狂的想法:想开始写周记。

之所以我觉得疯狂,是因为别说周记了,月记都难得写一篇,上一篇博客更新还是半年多以前。但真的不知道是什么原因,让我突然想有规律地记录一些东西下来。

回想了一下之前写的博客,基本上是按主题成文的,比如有了些新的想法,自然就会写篇博客出来。但当你接触的东西越来越多,新鲜的想法反而会越来越少,所以更新的进度也越来越慢。如果要求自己按一定的时间间隔写一些东西出来,或许还能把思路放开一些,少一些拘束。这大概是潜意识里准备写周记的一个原因。

另一件事情我发现让我尤其苦恼,就是很多时候当我确实想写博客的时候,我会在脑子里快速地构思一遍要写的内容,但这个过程结束后,我也就失去了把想法写下来的动力。就在二十分钟前,我还准备用一晚上的时间读掉一本小说,但有了写周记的冲动后,我便立马开了机,生怕又像以往那样,计划构思完了就又把它们抛在脑后了。

回到写周记这事,最早应该是在初中的时候语文老师让我们写的,目的大概是锻炼写作能力之类的,而且内容通常都是写一件小事,然后在文末点题,抒发一下人生感悟什么的。用的作业本是那种半张A4纸那么大,画了方格的,每周老师要收上去改,偶尔还会写些评语。现在想想,当真是幼稚又温馨。如果不出意外的话,家里某个角落里应该还能翻出以前用过的课本,写过的作业,考过的试卷一类,说不定还能找到些十几岁时的周记本。记得某次暑假回家,在杂物堆里赫然发现了<10岁时画过的画,当时激动得真的快要哭出来。当时的心情,是不是就是十几年后再看这篇博客时的心情呢?时间真是一件让人着迷的东西。

这周的事情说多不多,说少不少,挑三件来写一写。

周二的时候,系里举办感恩节聚餐(当然真正过节是下周),每个人带一道菜,拼在一起吃。我做了烤排骨,早上把菜送到教务秘书那儿就听课去了。等下课到就餐的地方,排骨已经没了一半,于是拖着餐盘小心夹走三块不声张,深藏功与名。我座位对面是我现在待的统计咨询中心的首席指导教授,寒暄了几句发现不知道聊啥,于是便默默吃饭(他边上还坐着系主任)。几分钟后边上的美国同学开始聊体育,各种橄榄球棒球,教授便滔滔不绝和他以及另外一个美国教授攀谈起来。我只好在心里默默摊一摊手,深感文化隔阂一直有,融也融不进去。

周三听了威斯康星大学Ming Yuan教授的报告,并注册了和他吃午饭。他的名字可以经常在各大统计期刊的理论文章中见到,而我对纯理论的研究一向有些成见,认为他们不关注实际问题。但接触了袁教授后我的这种观点有了180度的转变,他的演讲在开头用了足足十多分钟的时间来介绍一个具体的问题,然后一步一步进行归纳总结,最后提炼抽象出核心的理论部分,招式精妙,内功深厚,俨然一派大师风范。由此觉得艰深的内容不一定就出不了精彩的呈现,关键是看演讲者本身的功力是否足够。

今天上午上刘传海教授的大数据计算课,收获颇丰。跟刘教授有关的故事太多,一篇周记放不下,以后得慢慢写。总而言之一句话,他是我接触过的最具黑客精神的统计学教授。据他说他最近几个月几乎把 Email 给关了,他的一些同事朋友因为联系不上他还特意打电话询问情况,实际上他只在干一件事:写!程!序!要知道他可是正教授啊,还是统计学的正教授,就算是统计系的学生,据我所知也很少有这么投入的。那么他在写些什么呢?他在写R!的!源!代!码!几个月前他就提到他在研究R的核心,当时我以为他只是说笑,万万没想到居然是真的。那么他研究R的源代码是要干嘛呢?答案居然是要写一个多!线!程!的!R!平时大家最恨R的哪点?不就是它单线程吗?那么为什么没有人把R改成多线程的呢?那是因为这要触碰到R的源代码,而这些底层的东西过于艰深,很少有人有能力也有精力去做这方面的研究。目前来说,R里面实现并行的方式基本都是要通过创建新的 R 进程 来实现,而他今天给我们展示的结果,是让R可以通过直接创建 线程 来并行地执行R代码。关键是,今天的演示还成!功!了!当这个功能稳定以后,R在并行计算中的可扩展性将大大提高。上完课以后,有种多年未能攻破的难题终被化解了一样,此事以后详谈。

写到这儿发现已经两个小时过去了,小说是看不成了,以后也不知道能不能坚持写下去,先就此收尾吧。

阅读全文→

showtext 新版本发布

2015 年 04 月 29 日 | 分类于 r学习中
quick_example

前两天把 showtext 包更新到了 0.4-1 版本。因为这次更新的内容比较多,所以简单写篇博客介绍下。

第一个比较重大的改变是现在可以在图形设备中自动使用 showtext 了,用法比以前更为简单。在这个版本之前,需要每次都打开图形设备,调用 showtext.begin()showtext.end(),然后关闭图形设备。这样反复调用函数显得非常麻烦。在新版本中,只要在最开始调用一下 showtext.auto(),之后所有的图形设备就会自动使用 showtext 来进行渲染。这对窗口设备尤其重要,因为我们平常用 R 画图时,很少手动去调用 windows() 或是 x11(),而是直接使用 plot() 函数让 R 自动打开一个窗口。调用了 showtext.auto() 之后,窗口设备也可以自动使用载入的字体了,这对于预览图形效果会很有帮助。

自动调用的另一个好处是能在 ggsave() 中使用 showtext 了。ggsave() 本身会自动打开和关闭图形设备,所以在以前的版本中我们无法在 ggsave() 中插入 showtext.begin()。现在自动化以后,ggsave() 也可以使用 showtext 的字体库了。

要关闭自动调用功能,只需执行 showtext.auto(FALSE)

第二个比较显著的改进是更好地支持了位图图形。在之前的版本中,showtext 主要用于矢量图,比如 pdf() 或者 svg()。如果在 png() 或是 jpeg() 中使用 showtext,你会发现画出的字体非常难看。这是因为那些位图设备没有较好的抗锯齿支持。而现在,showtext 可以为那些设备绘制好平滑后的字形,所以图片质量也会得到提升。唯一需要注意的是要保持图形设备的 DPI 与 showtext 的一致,比如用 png() 设备时,应该手动设置一下想要的分辨率(此处 DPI 为120)

library(showtext)
showtext.opts(dpi = 120)
png(..., res = 120)

支持位图图形也就意味着支持了窗口设备,比如 windows()x11()

最后一个改进其实是之前 sysfonts 包的更新内容,就是目前在国内无法直接连接 Google 的字体库,所以在使用 font.add.google() 时,可以使用 360 提供的代理。使用方法是加一个参数,比如

library(showtext)
font.add.google("Gochi Hand", "gochi", repo = "useso")

其他的例子和说明都可以参见 Github 上的介绍

阅读全文→

所有文章列表→