Sep 302016
 

卡方检验的公式(2×2)
A short-cut formula for Chi-square test a 2 by 2 table

最近发现一个快速计算卡方检验的公式。
对于2×2的表格

\(
\begin{array}{|c|c|c|}
\hline
a & b & (a+b)\\ \hline
c & d & (c+d)\\ \hline
(a+c) & (b+d) & n\\ \hline
\end{array}
\)
上式中\(n=a+b+c+d\)

卡方检验的统计量是:
\(
\chi^2_1 = \frac{n(ad-bc)^2}{(a+b)(c+d)(a+c)(b+d)}
\)

这个方法和R计算的结果(不要用continuity correction, e.g. chisq(…, correct = FALSE))是一致的。

-- February 17, 2017更新 --

刚刚发现Matthew’s Correlation Coefficient (MCC)的定义和2×2 table的卡方检验有如下联系
\(
|MCC| = \sqrt{\frac{\chi^2}{n}}
\)
因为MCC的范围是-1到1,因此\(\chi^2 \)的范围是0到\(n\)

Aug 122015
 

R的诡异错误 (3)
A pitfall in R (3)

下面的代码有两行,第一行计算Genotype的Minor Allele Frequency,第二行计算Weight matrix (用SKAT paper建议的Beta分布的weight)。
下面的代码有问题吗?

# G is a sample by marker matrix
maf <- colMeans(G) / 2
w <- diag(dbeta(maf, 1, 25))

粗粗一看,上面的代码的思路是正确的,不过这里有两个坑。
第一个坑是缺失数据(Missing data),如果Genotype矩阵G有缺失数据,那么maf和w都会有缺失的数据;
第二个坑比较难看出来:当G只有一列的时候,maf也会是一个标量(Scalar),同样dbeta函数的返回值是标量。
你可能想diag函数作用在一个标量上,其返回值是一个1×1的矩阵。
但是,diag()的返回的是一个….不确定的结果,看个例子吧:

diag(3)            => 3x3 单位矩阵
diag(3.9)          => 还是3x3 单位矩阵
diag(0)            => 是0x0 矩阵
diag(-1)           => 报错
diag(c(1, 2, 3)    => 3x3 矩阵, 对角线是1,2,3

总之,第二个坑告诉我们,如果dbeta返回的是一个长度为1的数值,别指望w一定是1×1的矩阵。

继续吐槽:
diag()这个函数设计的不好用。在R内部,标量就是被实现为长度为1的向量。
在现在的设计下,diag要求使用者同时考虑变量的类型(矩阵,向量,缺失)以及变量的长度和值(向量长度为1且大于等于0,向量长度为1且小于0,向量长度大于1)。
对于一个适应了使用有类型系统的编程语言的程序员,这种要求会增加使用者的脑力负担,同时阅读者在语义上也很容易出现混淆。
我建议的解决方案是:diag(向量)永远输出一个对角线是该向量的矩阵,同时定义一个新函数(可以叫diag_len,类比seq_len;或者叫eye,类比Matlab里的eye函数)专门用来生成单位矩阵。

最后总结一下,我觉得这样写比较

# G is a sample by marker matrix
maf <- colMeans(G, na.rm = TRUE) / 2
m <- length(maf)
w <- matrix(0, nr = m, nc = m) 
diag(w) <- dbeta(maf, 1, 25) 
# 或者:w <- diag(dbeta(maf, 1, 25), nrow = m, ncol = m)
Jul 112015
 

深度学习MNIST数据

Deep learning on MNIST data

最近试了一下Deep Learning算法,在Kaggle的Digit Recognizer排行榜上目前排19名。

简单说说经验。

从工具上讲,建议先用Keras来做原型(代码简单,很容易构建复杂的网络,缺点是太占内存/显存),再用Caffe做更细致的调节(计算快,占内存/显存少,但用起来复杂,比如不支持RMSProp,需要手动把数据分成Training/Testing,要手动用proto buffer来构建神经网络)。此外,用GPU加速是必须的,我观察到GPU速度至少是CPU速度的60倍。

从模型上讲,基本上建立不同的Convolution Network。套路似乎都是重复这个结构:Convolution Layer + ReLU Layer + MaxPooling Layer。模型可以更Deep (层数更多),比如VGG网络,但是代价是参数数量多,计算量太大。对于MNIST这个整理好(Pre-processed)的数据,这些网络的效果很相似,一般默认设置都能达到98%正确率以上(用42k数据训练,28k数据测试)。如果要更好的成绩,基本上需要更多的输入数据(比如把原始的图像上下左右移动),更深的模型(比如增加层数),更多的模型(比如重复5遍,再Ensemble结果),这样一般能稳定达到99%以上的正确率。

目前能达到99.486%正确性的模型是:把已有的42k数据按照80%, 20%的比例分成training,testing数据,然后用Lenet(120 convolution layers [kernel = 5] + ReLU + MaxPooling + 200 convolution layer [kernel = 3] + ReLU + MaxPooling + InnerProduct [param = 200] + InnerProduct [param=10]), 迭代60,000次。然后这个过程重复15次,然后选多数Ensemble。

最后放几个神经网络识别起来费劲的图案供欣赏,还有一个介绍性质的Slides。

 

DigitRecognizerRank

DigitRecognizerRank

probably 4 looks like 9 3 or 2 7 or 2

Apr 012015
 

R的诡异错误(2)
A pitfall in R

R里面的陷阱一个接一个,继上次说的诡异错误,最近我又跌进了另一个坑。

先创造一个叫df的data.frame, 如下所示:

options(stringsAsFactors = FALSE)
df <- data.frame(
  group = c("group1", "group2"),
  val = c(100, 9)
  )
df
# df content:
#
#   group val
#1 group1 100
#2 group2   9

我们统计一下val列在每组(group)的均值:

## get mean value per group
library(plyr)
res <- ddply(df, .(group), function(x) {
  c(group = unique(x$group), 
    val = mean(x$val))
  })
res

res的结果如下:

> res
   group val
1 group1  100
2 group2  9
> res[1,2] > res[2,2]
[1] FALSE

困惑的地方来了,res[1,2]显示的是“100”, res[2,2]现实的是“9”, 为什么res[1,2] > res[2,2]的结果是FALSE?
R肯定是认为100比9大。这显然是不对的。

那上面哪个地方有陷阱呢?再仔细看一下高亮显示的两行,看出问题了么?

## get mean value per group
res <- ddply(df, .(group), function(x) {
  c(group = unique(x$group), 
    val = mean(x$val))
  })

陷阱在于第3-4行,当用c()来构造一个向量时,val被自动从数值型转成字符型。也就是说,存储的是“100”和“9”,不是100和9。
在字符型数值的比较中,“9”是会比“100”大的。因此order输出的是2 1.

从这个坑里爬出来,不由得吐槽一下,R的自动类型转换能减少程序的代码量,但是这种类型转换很容易引入错误。
也就是说R编程时不需要程序员去检查类型,但在底层R会自动的(有时也是自作聪明的)转换类型。
这就要求程序员额外费心考虑类型之间的自动转换,实际上加重了编程的负担。
如果不然,程序员放任R做自动的类型匹配和转换,就可能得到错误的结果,这在科学研究中实在是很危险。

既然R的坑这么多,我想过去试试新出的Julia或者Go语言。可惜现在这些语言都没有很好的支持统计计算,比如Julia目前缺少对缺失数据(missing data)的支持,
比较有希望的NumPy写起来语法繁杂,近期吸引不到大多数的统计学家(统计的用户太少,机器学习的用户多)。
要是能有一种计算工具,可以让统计学家很容易的使用,同时又有较强的类型系统,还能把大多数R包迁移过去,那这种工具必有远大前程。

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).

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

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要好一点。

Sep 262013
 

Ridge regression 相关的资料

以前学Ridge regression,大概记住了如何求回归系数 \(\beta^{ridge}\),但有些细节不了解。我找到下面两个slide感觉又学到了新东西。

 

第一个讲义 【Link or RidgeRegression – Ryan Tibshirani

http://www.stat.cmu.edu/~ryantibs/datamining/lectures/16-modr1.pdf

Slide 6:Shrink之后的回归系数看起来什么样?

Slide 7:为什么要center和normalize predictors?

Slide 13:当真正的回归系数很大的时候,Ridge regression还有好的表现么?或者说当true coefficient很大的时候,应该向什么方向调整\(\lambda\)可以减少prediction error?

 

第二个讲义 【Link or RidgeRegression-Rudyregularization

http://www-stat.stanford.edu/~owen/courses/305/Rudyregularization.pdf

Slide 17:Ridge regression的回归系数的bias具体是多少?

Slide 23-24:Ridge regression和Principal component analysis (PCA)的关系是什么?

Slide 37:选择Ridge regression的tuning parameter时一般用cross validation。如果用GCV(Generalized cross validation)的话,可以有close form.

Sep 082013
 

矩阵最大特征值的求法

Largest Eigenvalue of Matrix

 

矩阵的特征值分解是在科学计算常见的线性代数方法。在一些情况下,我们不需要全部特征值,只需要一部分特征值。或者因为矩阵的一些性质,我们只对一部分特征值感兴趣。比如对于矩阵 B = A’ * A, 如果A是m by n的实矩阵,m很大,n很小。根据线性代数理论,A, B矩阵的Rank最多是n,并且B矩阵是半正定矩阵。对矩阵B做特征值分解时,其特征值都是非负实数,并且最多有n个正数。因此,如果求解B矩阵的特征值,只需要知道其最大的n个特征值,因为其余特征值都是零。

求解一般矩阵的特征值主要有两种:(1)特征多项式(2)Power Method(包括 Arnoldi方法和QR分解)。 对于实对称正定矩阵,可以使用Choleskey分解。这些方法可以求出所有的特征值。而对于我们特别关心的实对称(Hermitain/Symmetric)正定(Positive Definite)矩阵,比如协方差矩阵,这里我们关心的是其最大(或者最小)的特征值。这里常用的方法的方法叫做Linczos方法,可以看做是Arnoldi方法的特例。Linczos方法计算A的特征值是(1)先随机生成一个向量v,(2)再利用一系列向量: v, A*v, A*A*v, …, A^(n-1) * v 来得到一个特殊矩阵T,这个矩阵的大小是m by m,这里m往往远大于n。(3)得到T之后,可以以O(m^2)的复杂度来的到特征值。

实现Lanczos方法时会使用Multiple restart,最有名的实现是ARPACK。ARPACK是用Fortran实现的,对Fortran不熟悉就不容易直接调用ARPACK的函数。不过很多其他软件都对ARPACK的功能进行了封装。在R里面使用ARPACK,可以用igraph包,R-blogger有一个这方面的例子(链接)。其他软件中一般也都有用ARPACK来实现求一部分特征值的功能:Matlab的eigs,Scipy的eigs和eighs(分别对应实对称矩阵和Hermitian矩阵, 链接),Julia的eigs(链接)。

后话:现在开源的科学计算软件很多,当我们想寻找合适高效的数值计算工具时,可以从源代码出发。比如我们知道Matlab有eigs函数,就可以用“julia eigs”来看julia语言里是怎么实现eigs的。