卡方检验的公式(2×2)

卡方检验的公式(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\)

Makefile 笔记

Some Notes about Makefile

记录一下最近用Makefile的经验和教训
内容比较凌乱。放在这里供以后参考。

1. Target-specific variable的作用域

Target-specific variable指的是变量在某些特定目标下的值。用法是:

target … : variable-assignment

注意这里的变量作用域只是在Recipe部分,而不能在目标(Target)或者依赖(Prerequisites)部分。

target … : prerequisites …
recipe

就是说这个Target-specific变量的值仅在recipe中有效。

2. 除错的方法

用warning 函数,比如$(warning msg)。它可以出现在各个部分: 目标, 依赖,recipe,顶级变量, 都行 (参见[2]):

$(warning A top-level warning)

FOO := $(warning Right-hand side of a simple variable)bar
BAZ = $(warning Right-hand side of a recursive variable)boo

$(warning A target)target: $(warning In a prerequisite list)makefile $(BAZ)
$(warning In a command script)
ls
$(BAZ):

输出:
$ make
makefile:1: A top-level warning
makefile:2: Right-hand side of a simple variable
makefile:5: A target
makefile:5: In a prerequisite list
makefile:5: Right-hand side of a recursive variable
makefile:8: Right-hand side of a recursive variable
makefile:6: In a command script
ls
makefile

还可以配合这些参数–just-print (-n), –print-data-base (-p), and –warn-undefined-variables

3. 管理大型项目

大型项目的代码有层级结构。xieMakefile要注意代码精简,善用函数,尽量自动化,避免重复代码,检查变量的有效性。

例如:
create-output-directories := \
$(shell for f in $(modules); \
do \
$(TEST) -d $$f || $(MKDIR) $$f; \
done)

例如用call, eval, foreach等函数:

# use foreach and test
.PHONY: validate_build
validate_build:
test $(foreach f,$(RELEASE_FILES),-s $f -a) -e .

local_src := $(addprefix $(subdirectory)/,playlist.y scanner.l)

# create library
define make-library
libraries += $(BINARY_DIR)/$1
sources += $2

$(BINARY_DIR)/$1: $(call source-dir-to-binary-dir, \
$(subst .c,.o,$(filter %.c,$2)) \
$(subst .y,.o,$(filter %.y,$2)) \
$(subst .l,.o,$(filter %.l,$2)))
$(AR) $(ARFLAGS) $$@ $$^
endef
$(eval $(call make-library, $(subdirectory)/libdb.a, $(local_src)))


参考
[1]Target-specific variables. GNU make. (link)
[2]Debugging Makefiles, Managing Projects with GNU Make, 3rd Edition, Chapter 12. (link)
[3]C and C++, Managing Projects with GNU Make, 3rd Edition, Chapter 8. (Link)

解决一个奇怪的网络连接错误

解决一个奇怪的网络连接错误
Solve a strange network problem

我的服务器bunny最近出现了一个奇怪的网络连接错误:只有台式机和服务器能连接到bunny,笔记本就连不上。
假设网络IP如下:

A. 服务器bunny:198.215.54.48 10G 网络
B. 服务器: 198.215.54.5 10G 网络
C. 服务器bronco: 129.112.7.169 1G 网络
D. 台式机:129.112.185.246 局域网
E. 笔记本:172.17.157.121 无线网

现在的问题是可以从B, C, D 连接到A,但不能从E连接到A。

解决思路如下:

1. 怀疑E->A的路由有问题

通过traceroute,发现E到A的最后一跳是*,状态是Host unavailable.
但进一步发现E->B,A和B在同一个子网,因此E->A的路由不大可能有问题

2. 检查A的防火墙设置

先关闭所有防火墙,用:

ufw disable

但E仍然不能连A

3. 检查A的路由

> route -n
Kernel IP routing table
Destination Gateway Genmask Flags Metric Ref Use Iface
0.0.0.0 198.215.54.254 0.0.0.0 UG 0 0 0 em1
172.17.0.0 0.0.0.0 255.255.0.0 U 0 0 0 docker0
192.168.1.0 0.0.0.0 255.255.255.0 U 0 0 0 em4
198.215.54.0 0.0.0.0 255.255.255.0 U 0 0 0 em1

> ip route --list
ip rule list
0: from all lookup local
32766: from all lookup main
32767: from all lookup default

仔细检查发现172.17.0.0这个路由有问题:因为E在172.17.0.0网络,从E发出的ICMP包到达A,而从A返回的ICMP包没有通过em1接口而是docker0发出,这样E->A就显示了Host unavailable。

值得注意的是在路由IP包时,它的终点并不是第一个可以匹配的路由,而是所有路由里的第一个匹配最好的路由。
原文:Any entry whose first field matches the destination IP address completely(a host) or partially (a network) would signal the IP address of the next router. (link)

找到原因后,我们可以让docker0的网络IP地址避开172.17.0.0子网,这就可以解决E->A的连接问题。
这个解决方案里的重要命令是 (参考链接


ip link set dev docker0 down
ip addr del 172.17.42.1/16 dev docker0 # 172.17.42.1是docker0的ip
ip addr add 10.0.0.10/24 dev docker0
ip link set dev docker0 up
ip addr show docker0 #检查docker0的ip

此外要在/etc/default/docker里加上DOCKER_OPTS (–bip参数见Docker文档Customize the docker0 bridge)。这样docker0的IP不会在服务器重启后改变。


DOCKER_OPTS="--bip=10.0.0.10/24"

题外话

1. 为什么docker和无线网都用172.17.0.0子网?

根据RFC 1918 3. Private Address Space, 保留IP除了常见的192.168.0.0外, 还有:

10.0.0.0 - 10.255.255.255 (10/8 prefix)
172.16.0.0 - 172.31.255.255 (172.16/12 prefix)
192.168.0.0 - 192.168.255.255 (192.168/16 prefix)

因为172.17.0.0是保留地址,正巧docker和无线网同时都使用了这个子网,因此造成了我遇到的网络问题。

2. 怎么监测服务器是否接收到IP包?

可以使用iptables

例如参考 http://askubuntu.com/questions/348439/iptables-log-file-and-how-change-it, 记录192.168.11.0/24发来的包:


iptables -A INPUT -s 192.168.11.0/24 -j LOG --log-prefix='[netfilter] '

可以在/var/log/kern.log 找含有“[netfilter]”的日志。

如果想监测ICMP包,可以用:


iptables -A INPUT -p icmp -j LOG --log-prefix='[netfilter] '

其他常用命令(参考这个链接)包括:


sudo iptables -L # 打印iptables内容
sudo iptables -L INPUT -v # 检查INPUT表收到了多少数据
sudo iptables -Z # 清零
sudo iptables -Z INPUT # 清零INPUT表数据清零
sudo iptables -A INPUT -p icmp -j LOG --log-prefix='[netfilter] ' # 记录ICMP包
sudo iptables -D INPUT -p icmp -j LOG --log-prefix='[netfilter] ' # 删除防火墙规则(记录ICMP包)
sudo iptables -F #防火墙规则生效

以下命令重置防火墙

sudo iptables -P INPUT ACCEPT
sudo iptables -P FORWARD ACCEPT
sudo iptables -P OUTPUT ACCEPT
sudo iptables -t nat -F
sudo iptables -t mangle -F
sudo iptables -F
sudo iptables -X

iptable是个复杂的防火墙,它使用了table,chain等概念,比如下图:

iptables Overview
iptables Overview

(来自链接
更详细的介绍可以参考[2]。

除iptables外,Ubuntu还提供了UFW,这个工具建立在iptables上,提供比iptables更简易的配置。比如:


ufw enable # 启用防火墙
ufw disable # 关闭防火墙
sudo ufw status # 显示防火墙设置
ufw deny 80/tcp # 封锁tcp的80端口
sudo ufw delete deny 80/tcp # 删除上一条(封锁tcp的80端口)

更多的信息可以参考Ubuntu UFW文档 (链接)。

 

参考

[1]7 Linux Route Command Examples (How to Add Route in Linux)

[2]Linux Firewall Tutorial: IPTables Tables, Chains, Rules Fundamentals

R的诡异错误(3)

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)

深度学习MNIST数据

深度学习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

不要用C++11特性开发R包

不要用C++11特性开发R包
Don’t use C++11 features to develop R packages

这是个教训。在我开发seqminer的时候浪费太多时间了。希望其他人不要重复我的错误。

C++11有不少有意思的特性,但为什么不要在开发R包的时候使用C++11的特性呢(这里指R 3.2)?下面有一些原因。

1. 没有清晰的文档
R对于C++11的支持写的很模糊。你需要看两个文档,Wring R ExtensionsR Installation and Administration
第一个文档提到你需要CXX_STD,第二个文档提到CXX1XFLAGS。这时候你知道应该设置一些变量,但怎么设置?没有文档告诉你。

这里我找到一个方案,可以凑合用。你需要改configure.ac文件。不过你仍会碰到下面提到的第3个问题。

dnl https://github.com/cran/dfcomb/blob/master/configure.ac
dnl copy the following lines to configure.ac

AC_DEFUN([AX_GET_R_CONF_VAR],[$1=`"${R_HOME}/bin/R" CMD config $1`])
AX_GET_R_CONF_VAR([CC])
AX_GET_R_CONF_VAR([CFLAGS])
AX_GET_R_CONF_VAR([CPPFLAGS])
AX_GET_R_CONF_VAR([CXXCPP])
AX_GET_R_CONF_VAR([CXX1X])
AX_GET_R_CONF_VAR([CXX1XSTD])
AX_GET_R_CONF_VAR([CXX1XFLAGS])
AX_GET_R_CONF_VAR([CXX1XXPICFLAGS])

CXX="${CXX1X} ${CXX1XSTD}"
CXXFLAGS="${CXX1XFLAGS} ${CXX1XPICFLAGS}"
CXXPICFLAGS="${CXX1XXPICFLAGS}"

2. 跨平台支持很艰难

当你把包提交到CRAN时,很有可能在Solaris系统上出现编译或运行的错误。
如果你使用了C++11的功能,你可能想在自己安装Solaris试试。
从我的经验来看,Solaris11 + Solaris Studio 12.4环境比较好搭建。
但CRAN上用的是很老的Solaris10 + Solaris Studio 12.3。
如果你想搭建这个系统,你需要从源代码编译R,然后你会发现,编译不了!
我发现blas库无法被链接。因此你在提交R包到CRAN时,最好明确说,我不支持Solaris。

3. 支持Linux系统也会出问题

在Linux系统上,也许gcc的版本很新,可以支持C++11,但如果安装R的时候没有设置C++11需要的一些变量,你还是用不了。
比如Ubuntu系统里, R是用apt安装的,这个版本的R是默认不知道你有没有支持C++11的编译器的。
这就是说,你在开发时使用的R,CRAN上的R都可以安装你的包,但不保证你的用户能安装。

怎么知道你的R是否支持C++11呢?你可以运行“R CMD config”:

  CXX1X         C++ compiler command for C++11 code
  CXX1XSTD      flag used to enable C++11 support
  CXX1XFLAGS    C++11 compiler flags
  CXX1XXPICFLAGS
                special flags for compiling C++11 code to be turned into
                a shared library

然后看一下CXX1X的值, “R CMD config CXX1X”,如果这个值是空的。默认情况下没法编译安装你写的R包(install.packages())。
顺便说一下: CXX1XXPICFLGS是一个笔误,正确的拼写是CXX1XPICFLAGS (不是两个XX)

解决方法:
先在shell下:

export CXX1X=”g++ -std=c++11″
export SHLIB_CXX1XLD=”g++ -std=c++11 -shared”

然后R里面就能用install.packages()安装了。
顺便说一下,SHLIB_CXX1XLD 没有写在手册里……

现在总结一下,如果你用了C++11的特性写R包,你会花很多时间读手册让R去选用合适的编译器,很难去支持CRAN上的老旧的Solaris,不能保证最终用户能够编译安装你的R包。
所以,我的建议的是,如果想让更多的人用你的R包,别折腾C++11的特性去适应R,还是写最portable的C++代码吧。

安装GCC 4.9

安装GCC
Install GCC

这是我找到的最简单的安装gcc,并且不需要root权限的方法 (https://gcc.gnu.org/wiki/InstallingGCC

这个方法适用与GCC的多个版本(4.6,4.7,4.8,4.9)

tar xzf gcc-4.6.2.tar.gz
cd gcc-4.6.2
./contrib/download_prerequisites
cd ..
mkdir objdir
cd objdir
$PWD/../gcc-4.6.2/configure --prefix=$HOME/gcc-4.6.2
make
make install

第七行:如果只需要64位编译器,不需要32位的,可以加上–disable-multilib
第八行:可以用-j N参数指定N个进程同时编译

R的诡异错误(2)

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包迁移过去,那这种工具必有远大前程。

设置代理服务器

设置代理服务器
Setting up proxy

出于对患者信息的保护,学校内网接入Internet时必须使用代理(Proxy)。
这令很多软件的安装使用过程变得复杂。

这个帖子会列出给常用软件设置代理的方法。

R

最直接的设置代理的方法是设置环境变量: http_proxy, https_proxy, 例如:

Sys.setenv(http_proxy="http://proxy.swmed.edu:3128")
Sys.setenv(https_proxy="http://proxy.swmed.edu:3128")

如果上述方法不起作用,也有可能需要:

# In R, a partial solution is to use:
options(download.file.method="wget")

这是在curl不能正常使用代理,但wget可以用的时候的解决方案。

还有一种情况是用devtools安装R的扩展包,可以用:

# For devtools
library(httr)
set_config(
  use_proxy(url="proxy.swmed.edu", port=3128)
)

如果Mac电脑上没有安装X11,装package时需要指定repository(来源),比如:

install.packages('RMySQL', repos='http://cran.us.r-project.org')

Docker

很多软件提供了Dockerfile。用户需要使用docker build来制作自己的容器(docker)。
如果在防火墙后面, 简单的方法是用(link):

docker build 
  --build-arg https_proxy=$HTTP_PROXY --build-arg http_proxy=$HTTP_PROXY 
  --build-arg HTTP_PROXY=$HTTP_PROXY --build-arg HTTPS_PROXY=$HTTP_PROXY 
  --build-arg NO_PROXY=$NO_PROXY  --build-arg no_proxy=$NO_PROXY -t java .

另外一个方法是手动修改Dockerfile,加入ENV语句,例如:

ENV HTTP_PROXY http://proxy.swmed.edu
ENV HTTPS_PROXY http://proxy.swmed.edu
ENV http_proxy http://proxy.swmed.edu
ENV https_proxy http://proxy.swmed.edu
....

2015春节

2015 春节
2015 Chinese New Year

爆竹声中一岁除,春风送暖入屠苏。
千门万户曈曈日,总把新桃换旧符。

——王安石的《元日》

除夕夜里,稍稍反思一下2014年经历的投稿。

背景:我们写的软件速度快,功能全。比如在文章中,我们强调该软件可以处理很大规模的数据,这是已知其他软件达不到的。

审稿:除去离题万里的审稿意见以及错误的审稿意见,我们得到的反馈是对自己的软件过分称赞,对有竞争关系的其他软件过分贬低。同时这些意见往往对我们有强烈的负面情绪。

反思:比较软件从某种意义上就像是比较孩子。说自己的孩子好没问题,但要是说自己的孩子比别人的好,就易引起强烈的情绪反扑。即使是数字上的比较(比如身高说自己的孩子高于其他小孩XX厘米),有意见的人也可能会说:身高不是重要的标准,你比较的方式不对等等。这个经历告诉我指出其他人/事在科研上的不足应该非常谨慎。或许寻找“和气”、容易接受的比较方式是有必要的。

2015年我们会继续努力,写有价值的文章,做好科研,也希望能通过温和的方式,得到其他科研同行的肯定。此外,希望看到这里的各位 新春快乐!