标签归档:计算

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

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

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

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案例学习