A short-cut formula for Chi-square test a 2 by 2 table

$$\begin{array}{|c|c|c|} \hline a & b & (a+b)\\ \hline c & d & (c+d)\\ \hline (a+c) & (b+d) & n\\ \hline \end{array}$$

$$\chi^2_1 = \frac{n(ad-bc)^2}{(a+b)(c+d)(a+c)(b+d)}$$

－－ February 17, 2017更新 －－

$$|MCC| = \sqrt{\frac{\chi^2}{n}}$$

1. Target-specific variable的作用域

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

target … : variable-assignment 

target … : prerequisites … recipe … …

2. 除错的方法

 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

3. 管理大型项目

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 （来自链接 更详细的介绍可以参考［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) 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)  • August 12, 2015 • Tagged with: 深度学习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 不要用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++代码吧。 • June 16, 2015 • Tagged with: , 安装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个进程同时编译 • May 24, 2015 • Tagged with: 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包迁移过去，那这种工具必有远大前程。 • April 1, 2015 • Tagged with: , 设置代理服务器 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 .


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

•  March 5, 2015
•  Tagged with: ,

2015 春节
2015 Chinese New Year

——王安石的《元日》

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

•  February 19, 2015