zhanxw

Feb 032015
 

R函数式编程
Functional Programming in R

R语言的循环很慢,同时语法有很复杂,陷阱很多。举例来说,逐行遍历一个矩阵有两种写法:

(1)

for (i in seq_len(nrow(m))){print(i)}

(2)

for (i in 1:nrow(m)){print(i)}

这两个写法那个对?或者两个都对?
熟悉R的用户自然会用(1),因为当m矩阵的行数为零时,只有(1)的结果是正确的:

> m <- matrix(, nr = 0, nc = 0)
> for (i in seq_len(nrow(m))){print(i)}
> for (i in 1:nrow(m)){print(i)}
[1] 1
[1] 0

为了避免手动写循环,一种折中的解决方案是用Filter, Reduce, Map等函数式编程的概念。
这些函数可以比较自然的处理列表(list)或者向量(vector)类型。举个例子:

l <- as.list(seq(9, 1, -1))
d <- data.frame(a = c(1, 2), b = c(3, 4))
m <- matrix(seq(9, 1, -1), 3, 3)
v <- seq(9, 1, -1)

isEven <- function(x) { x %% 2 == 0 }
isEven(1:3)

Filter(isEven, l)
## Filter(isEven, d) ## this does not work
Filter(isEven, m)
Filter(isEven, v)

Position(isEven, l)
# Position(isEven, d) ## this does not work
Position(isEven, m)
Position(isEven, v)

Find(isEven, l)
# Find(isEven, d) ## this does not work
Find(isEven, m)
Find(isEven, v)

可以看到,对于列表(list)或者向量(vector)类型,R内建的函数式函数(更多的见[1])还是很方便的。
如果想把函数式编程玩出更多花样,可以看看Hadley的新书Adv-R(参见[2])。

这里还有个问题,这种函数式的写法速度如何?下面的代码告诉我们函数式编程节省了写程序的时间,在某些情况下也许也能节省计算时间。

library(microbenchmark)
v <- seq(10000)
m1 <- function(v) {
  v[v%%2 == 0]
}
m2 <- function(v) {
  res <- vector()
  for(i in v) {
    if (i %% 2 == 0) {
      res <- c(res, i)
    }
  }
  res
}
m3 <- function(v) {
  Filter(isEven, v)
}

microbenchmark(res1 <- m1(v),
               res2 <- m2(v),
               res3 <- m3(v), times = 1000)

下面是结果:

Unit: microseconds
          expr       min        lq       mean    median         uq       max neval
 res1 <- m1(v)   341.038   348.526   529.9475   356.225   380.6965  58359.46  1000
 res2 <- m2(v) 33687.987 38433.148 48014.3656 39791.000 41751.3955 303604.96  1000
 res3 <- m3(v)  6178.547  6648.369  8493.9365  6999.976  7445.0725 124781.90  1000

在上面的例子中,
第一种是向量化的赋值,速度最快。
第二种是用循环,速度最慢,
第三种是函数式,速度居中。

不过据我的经验,当v的维数很小时(比如只有100个数),函数式相对于循坏的优势越来越小,甚至比循环要慢。

题外话:据说R的核心是类似于Scheme,也是一种函数式语言。但是在速度上对函数式编程的支持看来还有很大的提升空间。

[1]https://stat.ethz.ch/R-manual/R-devel/library/base/html/funprog.html
[2]Advanced R http://adv-r.had.co.nz/Functional-programming.html

Dec 262014
 

拉格朗日乘数检验
Lagrange Multiplier test (Score test)

拉格朗日乘数检验,英文是Lagrange multiplier test,或者叫做Score test是一种常用的统计检验。
拉格朗日乘数检验的名称来源于这个检验用的是拉格朗日乘数的分布,见[2]。
Score test的名称则来自于Score本身。
为了写起来方便,下面都用Score Test来代替拉格朗日乘数检验。

Score Test,Likelihood Ratio Test和Wald Test的图形表示

假设似然函数\( L(.) \) 只有一个参数, 这三种检验可以在一张图里表示出来:
FAQ__How_are_the_likelihood_ratio__Wald__and_Lagrange_multiplier__score__tests_different_and_or_similar_
(摘自:http://www.ats.ucla.edu/stat/mult_pkg/faq/general/nested_tests.htm)

Likelihood Ratio Test计算的是\( \xi^R = 2 (l(\hat{a}) – l(0)) \), 就是两倍红色竖线的长度,这个统计量近似有自由度为1的卡方分布(假设只有一个自由变量)
这里的\(\hat{a}\)是最大似然估计(Maximum Likelihood Estimator):

$$
\xi^R = 2 (l(\hat{a}) – l(0)) = 2 log(L(\hat{a})/L(0)) \sim \chi^2_1
$$

Score Test只需要考虑\( 0 \) 这一点上似然函数的性质,在图像上通过蓝色部分表示(通过考察这一点的斜率和曲率)。
写成公式是:
$$
\begin{align}
斜率 &= U = l'(0) \\
|曲率| &= V = |l”(0)| = -l”(0) \\
\xi^S &= U’ * V^{-1} * U = \frac{l'(0)^2}{|l”(0)|} \sim \chi^2_1
\end{align}
$$

而Wald Test只需要考虑\(\hat{a}\)这一点上似然函数的性质,在图像上用绿色的横线表示的(通过比较绿色横线的宽度和在\(\hat{a}\)的曲率)。
$$
\begin{align}
斜率 &= l'(a) \\
\xi^W &= \hat{a}^2 * |l”(a)| \sim \chi^2_1
\end{align}
$$

Score Test等价于Likelihood Ratio Test和Wald Test

从渐进性质(Asymptotic property)上讲,Score Test和Likelihood Ratio Test(LRT),Wald Test是等价的。
详细的严格的证明可以参见[2],文献里有更一般的结论(不限于自由度等于一)。
但是在不严格的意思下,我们可以用一个简化的例子来展示这种等价性质。
仍然假定\( H0: a = 0 \) 以及 \( H_a: a \neq 0 \)。

对于Likelihood Ratio Test (LRT),统计量(test statistics)是两倍红线的长度:
$$
\chi^{R} = 2 * ( l(\hat{a}) – l(0))
$$

首先证明Wald Test和LRT是等价的。
$$
l(0) = l(\hat{a}) – \hat{a} l'(\hat{a}) + \frac{1}{2} \hat{a}^2 l”(\hat{a})
$$
根据MLE的性质,\(l'(\hat{a}) = 0 \)。因此:
$$
l(0) = l(\hat{a}) + \frac{1}{2} \hat{a}^2 l”(\hat{a})
$$
整理一下,得到:
$$
\xi^R = 2 (l(\hat{a}) – l(0)) = \hat{a}^2 * (-l”(\hat{a})) = \hat{a}^2 * |l”(\hat{a})|= \xi^W
$$

接下来证明Score Test和LRT是等价的。
根据泰勒展开和MLE的性质:
$$
\begin{align}
l(\hat{a}) &= l(0) + \hat{a} * l'(0) + \frac{1}{2}\hat{a}^2 * l”(0) \\
l'(\hat{a}) &= l'(0) + \hat{a} * l”(0) = 0
\end{align}
$$
第二行的式子可以写成:
$$
\hat{a} = – \frac{l'(0)}{l”(0)}
$$
第一行可以写成:
$$
\begin{align}
\xi^R &= 2* (l(\hat{a}) – l(0)) = 2 * \hat{a} * l'(0) + \hat{a}^2 * l”(0) \\
&= – 2 * \frac{l'(0)^2}{l”(0)} + \frac{l'(0)^2}{l”(0)} \\
&= \frac{l'(0)^2}{-l”(0)} \\
&= \frac{l'(0)^2}{|l”(0)|} \\
&= \xi^S
\end{align}
$$

在上面的推导中,我们用到了多次\(l”()\),这个函数一般是负数。对这个函数取期望就是Fisher Information。
而Fisher Information的倒数或者是逆矩阵,就是Score statistics的方差。

Score Test的形式

上面的例子都是检验一个参数,如果参数个数多于一,那么有更一般的结论。
假设参数是\(\theta = (\theta_1, \theta_2) \),这里\(\theta_1,\theta_2\)分别是长度为\(p_1,p_2\)的向量,
假设检验的是: \(H_0: \theta_1 = 0, H_a: \theta_1 \neq 0\).

首先计算\(H_0\)下的MLE,假设是: \(\hat{\theta} = (0, \hat{\theta}_2) \)。然后计算U和V
$$
\begin{align}
U &= Score = l'(\hat{\theta}) = l'( (\theta_1 = 0, \hat{\theta}_2) ) \\
V &= I^{11}_{\hat{\theta}} = (I_{11} – I_{12} (I_{22})^{-1} I_{22} )_{\hat{\theta}}
\end{align}
$$
这里\(I\)是Fisher information:
$$
I(\theta)
= \left[ \begin{array}{cc}
I_{11} & I_{12} \\
I_{21} & I_{22}
\end{array} \right]_\theta
= \left[ \begin{array}{cc}
I_{11}(\theta) & I_{12}(\theta) \\
I_{21}(\theta) & I_{22}(\theta)
\end{array} \right]
$$
最后,统计量\( U * V^{-1} * U’ \sim \chi_{p_1} \)。

为什么上面的\(V\)不是\(I_{11}(\theta)\),而是\(I_{11}(\theta)\)减去额外的一项呢?
可以认为检验\(\theta_1\) 实际上是在控制\(\theta_2\)的同时检验\(\theta_1\),
因此\(\theta_1\)是在\(\theta_2\)的residual space里面,因此它本身的变化(Variance)就小了。

举例来说(下一节的例子也会提到),假设\(H_0: b = 0, H_a: b \neq 0\)。
如果模型是\( Y = X b + \epsilon, \epsilon_{ii} \sim N(0,1) \)。
那么\( V = \sigma^2 (X’X) \).
如果模型是\( Y = X b + Z r + \epsilon, \epsilon_{ii} \sim N(0,1) \),
那么按照公式,第一种方法的得到的\( V = X’X – X’Z (Z’Z)^{-1} Z X = V_{XX} \)

另一种方法是求出X的Residual adjusting for Z:
$$
X_Z = (I-H) X = (I – Z (Z’Z)^{-1} Z’ ) X
$$
可以得到:
$$
V_{XX} = X_Z’ X_Z = X’ (I-H)’ (I-H) X = X’ (I – H) X = X’X – X’Z(Z’Z)^{-1}ZX
$$
这和第一种方法得到的\(V\)是等价的。

上面的例子假设误差是独立同分布的,我们也可以推广到generalized linear regression。
那么相应的\(V\)项就是\(X’WX – X’W (X’W^{-1}X)^{-1} W X’\),
其中\(W\)可以是已知的,或者通过mean function求出来的。
比如Logistic regression里\(W = \mu * (1-\mu)\)。

更一般的结论可以参考[1]。

一些例子

下面给出一些具体例子,都假设 \(H_0: \theta_1 = 0, H_a: \theta_1 \neq 0\)。

1. 简单的线性回归(Simple Linear Regression) :\( Y = X b + \epsilon, \epsilon_{ii} \sim N(0,\sigma^2) \)

$$
l(b, \sigma^2) = – \frac{n}{2} \log(\sigma^2) – \frac{(Y-Xb)'(Y-Xb)}{2 \sigma^2}
$$
可以求出\(\hat{\sigma}^2 = \frac{1}{n} Y’Y \)
$$
U_b = \frac{\partial l}{\partial b} = Y’X / \hat{\sigma}^2 \\
V_{bb} = -\frac{\partial l^2}{\partial^2 b} = X’X / \hat{\sigma}^2\\
$$
因此\(\xi^S = \frac{(Y’X / \hat{\sigma}^2)^2}{ X’X / \hat{\sigma}^2 } = \frac{ (Y’X)^2 }{(X’X)(Y’Y)/n} \sim \chi^2_1\).
当\(X\)和\(Y\)的均值都是零的时候(centered),\(\xi^S = n r^2\).
这里\(r\)是皮尔森相关系数(Pearson correlation coefficient)。

2. 一般的线性回归:\( Y = X b + Z r + \epsilon, \epsilon_{ii} \sim N(0,1) \)
$$
\begin{align}
U_b &= \frac{\partial l}{\partial b} = (Y – Z \hat{r})’X / \hat{\sigma}^2 \\
V_{bb} &= -\frac{\partial l^2}{\partial^2 b} = (X’X – X’Z(Z’Z)^{-1}Z’X) / \hat{\sigma}^2\\
\hat{\sigma}^2 &= \frac{ (Y-Z \hat{r})'(Y-Z \hat{r})}{n} \\
\hat{r} &= (Z’Z)^{-1} Z’Y
\end{align}
$$

上面的形式也许不太好记,不过如果定义\(H_Z = Z(Z’Z)^{-1}Z’\),则有:
$$
\begin{align}
U_b &= Y'(I-H_z)X \\
V_{bb} &= X'(I-H_z)X /\hat{\sigma}^2 ]\\
\hat{\sigma}^2 &= \frac{ Y’ (I-H_z) Y }{n} \\
\end{align}
$$
可见\( I-H_z \)在这里有关键作用,如果让:
$$
X_z = (I-H_z)X \\
Y_z = (I-H_z)Y
$$
那么上面的Score test等价于第一种情况: \( Y_z = X_z b + \epsilon, \epsilon_{ii} \sim N(0, 1)\).

此外还可以证明\(\) \xi^S < \xi^R < \xi^W [/latex],见[2]。 3. Logistic回归: [latex]logit(E(Y)) = X b [/latex] $$ \begin{align} l &= \sum_i y_i \log(p_i) + (1-y_i) * \log(1-p_i) \\ & = \sum_i y_i \log(\frac{p}{1-p}) + log(1-p_i) \\ & = Y'X b - \sum_i \log(1+\exp(X_i b)) \end{align} $$ $$ \begin{align} U &= \frac{\partial l}{\partial b} = Y'X - (\frac{\exp(Xb)}{1+\exp(Xb)})' X = (Y - \hat{Y})' X \\ V &= - \frac{\partial^2 l}{\partial b^2} = - \frac{\partial -X' \hat{Y}}{\partial b} = X' \frac{\exp(Xb)}{(1+\exp(Xb))^2} X = X' W X \\ \end{align} $$ 上式中[latex]W = diag(\hat{y} * (1-\hat{y})) [/latex]. 在[latex]H_0[/latex]下,[latex]b = 0 [/latex],因此: $$ \begin{align} U_b &= (Y - \frac{1}{2})' X \\ V_{bb} &= \frac{1}{4} X' X \\ \end{align} $$ 4. Logistic回归: [latex]logit(E(Y)) = X b + Z r[/latex] $$ \begin{align} l &= \sum_i y_i \log(p_i) + (1-y_i) * \log(1-p_i) \\ & = \sum_i y_i \log(\frac{p}{1-p}) + log(1-p_i) \\ & = Y'(X b + Zr) - \sum_i \log(1+\exp(X_i b + Z_i r)) \end{align} $$ 经过一些跳步: $$ \begin{align} U_b &= \frac{\partial l}{\partial b} = (Y - \hat{Y})' X \\ V_{bb} &= - \frac{\partial^2 l}{\partial b^2} = X' W X - X'WZ (Z'WZ)^{-1} Z'W X \\ \hat{Y} &= \frac{1}{1 + \exp(- Z \hat{r} )} \end{align} $$ 类似的,上式中[latex]W = diag(\hat{y} * (1-\hat{y})) [/latex]. [latex]\hat{r}[/latex]是在[latex]b=0[/latex]时的MLE估计。 上面这个形式和[3]里的公式是等价的。 [1] Chen, C.-F. Score Tests for Regression Models. Journal of the American Statistical Association (1983).doi:10.1080/01621459.1983.10477945 [2] Lin, D.-Y. & Tang, Z.-Z. A General Framework for Detecting Disease Associations with Rare Variants in Sequencing Studies. Am. J. Hum. Genet. 89, 354–67 (2011). [3] Liu, D. et al. Meta-analysis of gene-level tests for rare variant association. Nature Genetics 46, 200204 (2013).

Dec 112014
 

加快C++编译速度
Speed up C++ compiling speed

在开发rvtests时,常常需要重新编译整个项目,但是这是个耗费时间的事。

比如,编译最耗费时间的regression库,用“make”,需要50秒。
在多CPU的机器上,打开并行编译选项后,“make -j”,需要11秒(系统时间是11秒,用户时间是64秒)。

这就是说重新编译regression库,最快也要11秒。

今天我发现有一个好用的工具,ccache
可以把这个时间缩小到5秒 (加速编译速度一倍)。
方法是用:

make CXX="ccache g++"

ccache首次使用会稍微慢一点。因为它会缓存(cache)源程序。但之后编译的速度就变得飞快。
比如上面的方法用5秒。如果还想更快,打开“-j”选项,只需要不到1秒。

Nov 172014
 

最近学习和研究有点松懈,转一篇古文激励一下自己。
小时候读这篇文章的时候没什么感受,现在读了就觉得应该立刻干活。

为学一首示子侄
彭端淑

天下事有难易乎?为之,则难者亦易矣;不为,则易者亦难矣。人之为学有难易乎?学之,则难者亦易矣;不学,则易者亦难矣。

吾资之昏,不逮人也,吾材之庸,不逮人也;旦旦而学之,久而不怠焉,迄乎成,而亦不知其昏与庸也。吾资之聪,倍人也,吾材之敏,倍人也;屏弃而不用,其与昏与庸无以异也。圣人之道,卒于鲁也传之。然则昏庸聪敏之用,岂有常哉?

蜀之鄙有二僧:其一贫,其一富。贫者语于富者曰:“吾欲之南海,何如?”富者曰:“子何恃而往?”曰:“吾一瓶一钵足矣。”富者曰:“吾数年来欲买舟而下,犹未能也。子何恃而往!”越明年,贫者自南海还,以告富者,富者有惭色。

西蜀之去南海,不知几千里也,僧富者不能至而贫者至焉。人之立志,顾不如蜀鄙之僧哉?是故聪与敏,可恃而不可恃也;自恃其聪与敏而不学者,自败者也。昏与庸,可限而不可限也;不自限其昏与庸,而力学不倦者,自力者也。

选自《白鹤堂文集》,作于乾隆九年(公元1744年)。

Sep 082014
 

拉普拉斯方法(Laplace’s Method)

拉普拉斯方法又称为拉普拉斯近似(Laplace Approximation)。它可以用来计算一元或多元积分[1]。

举例来说,假设 \(f(x)\) 是一维函数,我们要计算
\( \int_{-\infty}^{\infty} f(x) \mathrm{d} x \).

如果\(f(x)\)形式很复杂,我们往往找不到定积分的公式(Close form)。
如果用数值方法来计算积分,计算量又很大。
所以要想个办法,得到一个比较精确的积分结果。

拉普拉斯方法可以适用于这种情况。我们先用泰勒展开(Taylor Expansion):

$$ f(x) \approx f(x_0) + f'(x_0)(x-x_0) + \frac{1}{2}f”(x_0) (x-x_0)^2 $$
如果选取\(x_0\)使得\(f'(x_0) = 0\),则可以进一步简化为
$$ f(x) \approx f(x_0) + \frac{1}{2}f”(x_0) (x-x_0)^2 $$

再引入一个假设,\(\int_{-\infty}^{\infty} f(x) \mathrm{d} x \)
可以被写成 \(\int_{-\infty}^{\infty} e^{f(x)} \mathrm{d} x \) 的形式,那么就有:

$$ \begin{align}
\int_{-\infty}^{\infty} e^{f(x)} \mathrm{d} x
&\approx \int_{-\infty}^{\infty} e^{f(x_0) + \frac{1}{2}f”(x_0) (x-x_0)^2} \mathrm{d} x \\
&= e^{f(x_0)} \int_{-\infty}^{\infty} e^{- \frac{1}{2}|f”(x_0)| (x-x_0)^2} \mathrm{d} x \\
&= e^{f(x_0)} \sqrt{\frac{2\pi}{|f”(x_0)|}}
\end{align}
$$

注意上式中\(f”(x_0)\)要取绝对值。因为\(f'(x_0) = 0\),\(x_0\)是\(f(x)\)的极值点。
对于概率密度函数而言,一般它也是最大值点(mode), 因此\(f”(x_0) < 0[/latex]. 从几何上讲,拉普拉斯方法是要用一个[latex]e^{-x^2}[/latex]形式的函数近似[latex]f(x)[/latex]. 或者说,要用一个长的像正态函数的函数取近似,这样的函数服从[latex]N(x_0, 1/f''(x_0))[/latex]. 下面举个例子来计算 $$\int_{-\infty}^{\infty} x^2 e^{-x^2/2} \mathrm{d} x $$. 这个积分的结果是[latex]\sqrt{2\pi}[/latex],这可以这样计算出来。 假设[latex]Z[/latex]是一个标准正态分布,那么: $$\int_{-\infty}^{\infty} x^2 e^{-x^2/2} \mathrm{d} x = \sqrt{2\pi} E(Z^2) = \sqrt{2\pi} Var(Z) = \sqrt{2\pi} = 2.507$$. 用拉普拉斯方法: $$\int_{-\infty}^{\infty} x^2 e^{-x^2/2} \mathrm{d} x = \int_{-\infty}^{\infty} e^{-x^2/2 + 2 \log{x}} \mathrm{d} x $$. 对[latex]f(x)[/latex]求导: $$ \begin{align} f(x) &= -\frac{x^2}{2} + 2 \log{x} \\ f'(x) &= -x + \frac{2}{x} \\ f''(x) &= -1 - \frac{2}{x^2} \end{align} $$ 令[latex]f'(x) = 0[/latex],可以解出[latex]x_0 = \sqrt{2}[/latex]. 我们有[latex]f(x_0) = -1 + \log{2}[/latex] 和 [latex]f''(x_0) = -2[/latex]. 积分结果为: $$ e^{f(x_0)} \sqrt{\frac{2\pi}{|f''(x_0)|}} = e^{-1 + \log{2}} * \sqrt{\frac{2\pi}{|-2|}} = \frac{2}{e} \sqrt{\pi} = 1.304 $$ 比较2.507 和1.304,我们发现理论结果和拉普拉斯方法的结果几乎差了一倍。 为啥差别那么大呢? 看下面的图: [caption id="attachment_564" align="alignnone" width="300"]Demo Laplace Method Demo Laplace Method[/caption]

我们发现要求的积分值是黑色曲线下的面积。拉普拉斯方法想用灰色曲线([latex]e^{-\frac{x^2}{2}}\))通过拉伸得到红色曲线,然后用红色曲线下的面积来近似积分。但是由于黑色曲线有两个峰值,这个近似显然不算成功。
所以拉普拉斯方法是有局限性的:被积分的函数有一个峰值,并且和正态曲线长的像。这种情况下的近似才能比较精确。

对于概率密度函数,(无论一维还是多维),大部分都是一个峰值。或者因为中心极限定理,统计量均值的分布和正态函数很像,也是一个峰值。因此拉普拉斯方法的用处还是很多的。此外,拉普拉斯方法计算很快。比如在线性混合效果模型中,拉普拉斯方法是用的最广的方法,也可能是唯一能在实际中使用的方法[2-4]。

[1]维基百科 http://en.wikipedia.org/wiki/Laplace’s_method
[2]Approximate Inference in Generalized Linear Mixed Models. N. E. Breslow and D. G. Clayton. Journal of the American Statistical Association Vol. 88, No. 421 (Mar., 1993) , pp. 9-25
[3]Variance component testing in generalised linear models with random effects. XIHONG LIN. Biometrika (1997) 84 (2): 309-326.
[4]lme4 package Douglas Bates et al. http://cran.r-project.org/web/packages/lme4/lme4.pdf (see nAGQ parameter)

Jul 062014
 

卡方检验
Chi-square test

卡方检验的全称是Pearson’s Chi-square test,这是统计中最重要的检验。
从用途来讲, 卡方检验可以做两件事:
(1)适配度检定(Goodness of fit),检查样本是否符合某种随机分布;
(2)独立性检定(Independence test),检查多个变量之间是否独立。

这篇不讲怎么用卡方检验,而是讲讲卡方检验的公式,推导,以及和似然检验的联系,目的是温故知新。

从公式来讲卡方检验:

$$ \chi^2 = \sum_i \frac{(O_i – E_i)^2}{E_i} $$
这里\(O_i\)是观察到的次数,\(E_i\)是期望的次数,
统计量\(\chi^2\) 服从自由度\(df\)的卡方分布(\(\chi^2\) distribution)。

假设有一个有可能是多项分布的随机变量\(X \sim Multinomial(N, p_1, p_2, \ldots, p_k)\)。
观察到\(k\)类的次数分别是\(x_1, x_2, \ldots,x_k\).
为了统计检验\(X\)是否服从这个分布,列出零假设\(H_0\) 为\(X\)服从多项分布,对立假设为\(X\)不服从这个多项分布,可以得到:

$$ \chi^2 = \sum_i \frac{(O_i – E_i)^2}{E_i} = \sum_i \frac{(O_i – N \cdot p_i)^2}{N \cdot p_i}$$

注意这个每一项的分母是\( N \cdot p_i \), 但根据中心极限定理(Central Limit Theorem),
$$\frac{O_i – N \cdot p_i}{\sqrt{N \cdot p_i \cdot (1-p_i)}} \rightarrow Normal(0, 1)$$

同卡方检验的公式相比较,这里面有两个问题:
1. 为什么卡方检验不是\(N\)个\(\frac{O_i – N \cdot p_i}{\sqrt{N \cdot p_i \cdot (1-p_i)}}\)的平方和?
2. 为什么\( \chi^2 \)的自由度是\( n-1 \),而不是 \( n \)?

回答这两个问题的一个方法是推导一下 \( \chi^2 \) 的分布。
根据[1],推导的思路如下:

1. 令\( Z_i = \frac{O_i – N * p_i}{\sqrt{N*p_i}} \),因此\( \chi^2 = \sum_i Z_i^2 \).

2.
根据中心极限定理,
\( Z_i \rightarrow \sqrt{1-p_i} N(0, 1) = N(0, 1-p_i) \)
\( Cov(Z_i, Z_j) = \frac{Cov(O_i, O_j)}{n*\sqrt{p_i*p_j}} = -\sqrt{p_i * p_j} \)

3. 为了得到\( Z_1^2 + Z_2+ \ldots + Z_k^2 \) 的极限分布,
我们先假设\( \vec{g} = (g_1, g_2, \ldots, g_r) \) 是i.i.d.的正态分布,
再假设单位向量 \( \vec{p} = (\sqrt{p_1}, \sqrt{p_2}, \ldots, \sqrt{p_k}) \).
当\( \vec{g} \) 向量 减去 它在\( \vec{p} \)的投影时,
\( \vec{g^p} = \vec{g} – ( \vec{g} \cdot \vec{p} ) \cdot \vec{p} \) 的分布和 \( Z_1, Z_2, \ldots, Z_k \) 是 一样的。
为了证明这一点, 可以证明 \( E(Z_i^2) \rightarrow (g^p_i)^2 \) 和 \( E(Z_i Z_j) \rightarrow g^p_i \cdot g^p_j \) 。
有了\( \vec{g^p} \), 就可以把\( \chi^2 \) 表示为 \( \chi^2 = \sum_i (\vec{g^p}_i)^2 \) 即 \( \vec{g^p} \) 长度的平方。

4. 下面计算\( \vec{g^p} \) 的分布
\( \vec{g} \) 的每个分量是i.i.d的正态分布。
根据正态分布的性质,对\( \vec{g} \) 的坐标进行正交变换后,在新坐标下,每个坐标仍然是i.i.d的正态分布 (见[2])
因此我们构造出一个特殊的正交变换使得\( \vec{g} \) 变换之后有一个坐标系和\( \vec{p} \)同一个方向,
因此可以把\( \vec{g} \) 表示为\( (g’_1, g’_2, \ldots, g’_k) \),并且
\( \vec{g^p} = (g’_1, g’_2, \ldots, g’_{k-1}, 0) \)。
因此\( \chi^2 = \sum_{i=1}^{k-1} (g’_i)^2 = \chi^2_{k-1} \)。

以上四个步骤可以从理论上证明卡方检验。另外一种思路是把似然检验(G-test)和卡方检验联系起来:
$$ \begin{align} G
&= 2 \sum_i O_i \log(\frac{O_i}{E_i}) \\
&= 2 \sum_i O_i \log(1 + \frac{O_i – E_i}{E_i}) \\
&\sim 2 \sum_i O_i \frac{O_i – E_i}{E_i} \\
&= 2 \sum_i \frac{O_i^2 – E_i * O_i}{E_i}\\
&= 2 \sum_i \frac{O_i^2 – 2 * E_i * O_i + E_i^2}{E_i}\\
&= 2 \sum_i \frac{(O_i – E_i)^2} {E_i}
\end{align}
$$
上面的计算中用到了\( \log(1+x) \sim x \) 和 \( \sum_i O_i = \sum_i E_i \)。
可见两种检验基本上是等价的(上面推导来自[3],但更简洁)。

然而,似然检验和卡方检验相比,似然检验需要计算对数,在皮尔森的年代,这是很不方便的。
卡方检验则便于手算,这或许是得到广泛应用的原因之一吧。

[1]Pearson’s Theorem http://ocw.mit.edu/courses/mathematics/18-443-statistics-for-applications-fall-2003/lecture-notes/lec23.pdf
[2]Orthogonal Transformation of Standard Normal Sample http://ocw.mit.edu/courses/mathematics/18-443-statistics-for-applications-fall-2003/lecture-notes/lec15.pdf
[3]The Two-Way Likelihood Ratio (G) Test http://arxiv.org/pdf/1206.4881v2.pdf

Apr 132014
 

R package的UBSAN测试
How to perform UBSAN tests on R packages

最近需要提交一个R package (seqminer)到CRAN,我们需要保证这个包代码通过UBSAN 测试(Undefined bahavior sanitizer tests,详见这个链接).
这个测试可以有效的检测结果不确定的指令(例如打印INT_MAX+1)。

目前只有Clang(版本3.3以后)才能支持UBSAN。为了对R package进行UBSAN测试,我们首先需要一个特殊的R版本,即用打开UBSAN支持(-fsanitize=undefined)的Clang编译的R。参考了网上一篇Blog,写一下简单步骤:

1. 下载编译工具

sudo apt-get install valgrind subversion r-base-dev clang-3.4 texlive-fonts-extra texlive-latex-extra
sudo apt-get build-dep cran-base

2. 下载R源代码(devel版本)

svn co https://svn.r-project.org/R/trunk ~/R-devel

3. 替换R-devel/config.site

CC="clang -std=gnu99 -fsanitize=undefined"
CFLAGS="-fno-omit-frame-pointer -Wall -pedantic -mtune=native"
F77="gfortran"
LIBnn="lib64"
LDFLAGS="-L/usr/local/lib64 -L/usr/local/lib"
CXX="clang++ -std=c++11 -fsanitize=undefined"
CXXFLAGS="-fno-omit-frame-pointer -Wall -pedantic -mtune=native"
FC=${F77}

4. 编译

cd R-devel
./configure --with-x=no --without-recommended-packages
make

这里–without-recommended-packages会省略一些R的packages(例如boot,nlme,Matrix) 以提升编译速度。如果不用这个选项,我发现Matrix package总是没法编译

5. 建立~/.R/Makevars

CC = clang -std=gnu99 -fsanitize=undefined -fno-omit-frame-pointer
CXX = clang++ -fsanitize=undefined -fno-omit-frame-pointer
PKG_LIBS = /usr/lib/llvm-3.4/lib/clang/3.4/lib/linux/libclang_rt.ubsan_cxx-x86_64.a

这里CXX 一定要用clang++,如果用clang的话会出现链接错误,比如找不到这个符号: _ZTVN10__cxxabiv117__class_type_infoE

这里还需要指定PKG_LIBS包括LLVM的库文件,不然就找不到UBSAN的符号。

6. 测试package

R CMD build seqminer/
R CMD check seqminer_2.7.tar.gz

UBSAN测试的结果会显示在:seqminer.Rcheck/tests/tests.Rout
如果有错误,你会看很长的一行, runtime error….

其他技巧

检测内存还可以用valgrind:

R CMD check --use-valgrind seqminer_2.7.tar.gz
R -d valgrind --vanilla < mypkg-Ex.R
R -d "valgrind --tool=memcheck --leak-check=full" --vanilla < mypkg-Ex.R
Mar 092014
 

A pitfall in R

先看下面的代码:

freq <- 0.29
altAllele <- freq * 200 * 2
print(altAllele) ## this will output 116
print(length(rep(1, altAllele))) ## what is the output of this??

最后一行输出的是115, 不是116。为什么?
答案是:

altAllele是一个numeric类型的数,尽管打印出来是116,但在电脑里存储的是一个比整数116稍小的浮点数。因此rep(1, altAllele)是一个长度是115的vector。

R的程序有很多陷阱。这个陷阱是因为R不需要定义数据类型,但它的底层还是有类型的。因此写R程序一定要注意这点。解决方法如下:

freq <- 0.29
altAllele <- as.integer(freq * 200 * 2)
print(altAllele) ## this will output 116
print(length(rep(1, altAllele))) ## what is the output of this??

Python里有没有类似的陷阱?结论是一样有 (Python 2.7.5 NumPy 1.8.0):

freq = .29
altAllele = freq * 200 * 2
len(np.repeat(1, altAllele)) ## output 115 here.

但是Python里打印altAllele的值不是116,而是115.99999999999999。这总比R打印出16要好一点。

Mar 012014
 

自动更新IP
Create dynamics DNS

办公室电脑的IP经常变动。为了能远程登录,我需要知道它的IP是什么。
我的解决方法是写一个小程序,这个程序把办公室电脑的IP绑定到域名 zhanxw.x64.me
然后,我让办公室电脑每小时运行这个脚本一次,这样zhanxw.x64.me就能指向办公室电脑了。

这里x64.me域名是由DNSDynamic (www.dnsdynamic.org)提供的.

下面是脚本,注意使用的时候要替换username和password。

#!/usr/bin/env python
import socket
ip=str(socket.gethostbyname(socket.gethostname()))
print ip

import urllib2, base64
request = urllib2.Request("http://www.dnsdynamic.org/api/?hostname=zhanxw.x64.me&myip=%s" % ip)
username = "######"
password = "######"
base64string = base64.encodestring('%s:%s' % (username, password)).replace('\n', '')
request.add_header("Authorization", "Basic %s" % base64string)   
result = urllib2.urlopen(request)

data = result.read()
print data

对于Windows主机,可以用计划任务每小时运行脚本一次。对于Linux主机,可以用crontab。

Dec 262013
 

gcc参数顺序
Ordering of gcc parameters matters

2013年最后一个月最后一个星期谈一个奇怪的bug,我在尝试pinfo的时候,configure脚本会提示ncurses有问题:

checking location of curses.h file... /usr/include/ncurses.h
checking if curses is usable... no
configure: error: Curses not found. You need curses to compile pinfo

打开config.log文件,发现这个提示:

configure:12563: checking if curses is usable
configure:12587: gcc -o conftest -g -O2 -I/usr/include    -L/usr/lib -lncurses conftest.c  >&5
/tmp/ccBLkSqd.o: In function `main':
/net/fantasia/home/zhanxw/software/pinfo-0.6.10/conftest.c:38: undefined reference to `initscr'
/net/fantasia/home/zhanxw/software/pinfo-0.6.10/conftest.c:39: undefined reference to `printw'
/net/fantasia/home/zhanxw/software/pinfo-0.6.10/conftest.c:40: undefined reference to `stdscr'
/net/fantasia/home/zhanxw/software/pinfo-0.6.10/conftest.c:40: undefined reference to `wrefresh'
/net/fantasia/home/zhanxw/software/pinfo-0.6.10/conftest.c:41: undefined reference to `stdscr'
/net/fantasia/home/zhanxw/software/pinfo-0.6.10/conftest.c:41: undefined reference to `wgetch'
/net/fantasia/home/zhanxw/software/pinfo-0.6.10/conftest.c:42: undefined reference to `endwin'
collect2: ld returned 1 exit status

从这个命令行来看,curses库的头文件目录和库文件目录都对,但编译的时候(链接/tmp/ccBLkSqd.o)却找不到几个函数的定义(比如:initscr,printw等等)。为了重现这个看起来不该出现的问题,我把conftest.c创建好,然后用相同的gcc命令行编译代码,的确是可以重现错误。

更令我不解的是如果把conftest.c放到命令行的末尾,这个编译就能通过了:

# This does not work
gcc -o conftest -g -O2 -I/usr/include    -L/usr/lib -lncurses conftest.c 
# This works
gcc -o conftest conftest.c -g -O2 -I/usr/include    -L/usr/lib -lncurses  

为什么把conftest.c移到“-L/usr/lib -lncurses”前面就能解决问题?
看到StackOverFlow上的一个帖子

这个帖子讲到gcc编译和链接一个源程序时,会保留参数的顺序,也就是说在链接时,把conftest.c放到前面和后面,效果不同。具体来讲:

conftest.c 放到后面

gcc -o conftest -g -O2 -I/usr/include    -L/usr/lib -lncurses conftest.c 
<=> 等价于
gcc -o conftest -g -O2 -L/usr/lib -lncurses tmporary_object_file.o

conftest.c 放到前面

gcc -o conftest -g -O2 -I/usr/include conftest.c -L/usr/lib -lncurses 
<=> 等价于
gcc -o conftest -g -O2 tmporary_object_file.o -L/usr/lib -lncurses 

也就是说,gcc在链接目标文件时,会依照命令行的顺序,对于每个在前面目标文件中没有定义的函数,都试着在其后面的目标文件找找,如果最后一个目标文件也找不到,就报错。因此在我们的例子里,tmporary_object_file.o一定要放到-L/usr/lib -lncurses前面才能正确编译。

回到本文最开始,我们如何让configure把conftest.c放到“-L”前面?答案不是改动configure,而是使用LIBS环境变量:

LIBS="-L/usr/lib -lncurses" ./configure

与LIBS不同,LDFLAGS会放到conftest.c前面,这样configure还是没法正确的使用configure。

在本文最后,回想我们本来的问题,就是说gcc的参数有一定的顺序要求。也许gcc的手册里写了顺序的要求,但是很少有人能够看完gcc的手册。如果见到类似问题,应该怎么办呢?我觉得可以有两个思路,第一个就是想想写gcc软件的人是怎么想的。比如现在gcc的选择就是一个很方便的方法,如果当前没某个函数没有定义,就在他后面制定的目标文件里找找。如果这个要求反过来,或者任何顺序都支持,那么链接的函数就会复杂(因为链接程序需要存储所有定义过的函数,参见链接);第二个就是想想传统的(或者标准的)写法是什么,把已知的经验和这个错误联系起来。比如configure的错误看起来很奇怪,但是先按照传统的写法先编译后链接,找到一个能凑合工作的方案(把conftest.c放到前面),再逐步改变到出错时的命令行,这个逐步的过程可以帮助查错也能巩固一下已有的知识。