Author Archive for 写长城的诗

一份数据挖掘小书评

书看得多了,渐渐有了品味,所以看那些英文版的时候也不再有很膜拜的感觉。不过再差的书,总会有些许长处或教训,值得借鉴。下面列出一些在数据挖掘实现方面的书,主要是基于R和python的。

  • Data mining with R Learning with Case studies: 我的启蒙书,很好的案例教学,这种书应该出更多些。
  • Data Mining Applications with R: 风格类似上面,多人合作的案例讲解,非常不错,就是每个例子略为简短了一点。
  • Machine Learning for Hackers: 也是以案例为主,但是内容质量不能保持在同一水准上,代码讲的比较细。
  • An Introduction to Statistical Learning: 比较系统的讲解了统计学习的内容,以及相应的R包函数,优点在于有习题。
  • Machine Learning with R: 同样是系统的列出各机器学习算法对应的R包函数,直观理解很好,而且讲了一些caret包的使用。比较新的书。
  • Applied predictive modeling: 对预测算法讲的最细,涉及R包最多,由于作者是caret的开发者,所以要学习caret则必看此书。
  • Big Data Analytics with R and Hadoop: 刚出的大数据挖掘书籍,讲解了如何用RHadoop来整合R和hadoop,什么时候会有R和spark的书出来啊。

  • programming_collective_intelligence: 有些年头的书了,但确实是经典,原理讲的简洁,就讲python代码实现,喜欢这种风格。
  • Machine Learning in Action: 和上面风格类似,不过实现是主要基于python的numpy库,这样代码量要精简很多,也很好。
  • Python for data analysis: 主要讲pandas库,精彩例子不多,可以直接去看帮助文档。
  • Learning scikit-learn Machine Learning: 主要以scikit-learn扩展库为工具做机器学习,比较简单,可以直接去看帮助文档。
  • Building Machine Learning Systems with Python: 和上面那本类似,讲函数功能为主,原理直觉没涉及到。

用模拟来理解混合效应模型之一:random intercept model

混合效应模型(Mixed-effect Model)在之前文章提到过,但感觉仍是雾里看花。此番又研究了一些资料,准备来做一个系列讲讲清楚,也算是自己学习的一个总结。

普通的线性回归只包含两项影响因素,即固定效应(fixed-effect)和噪声(noise)。噪声是我们模型中没有考虑的随机因素。而固定效应是那些可预测因素,而且能完整的划分总体。例如模型中的性别变量,我们清楚只有两种性别,而且理解这种变量的变化对结果的影响。

那么为什么需要 Mixed-effect Model?因为有些现实的复杂数据是普通线性回归是处理不了的。例如我们对一些人群进行重复测量,此时存在两种随机因素会影响模型,一种是对某个人重复测试而形成的随机噪声,另一种是因为人和人不同而形成的随机效应(random effect)。如果将一个人的测量数据看作一个组,随机因素就包括了组内随机因素(noise)和组间随机因素(random effect)。这种嵌套的随机因素结构违反了普通线性回归的假设条件。

你可能会把人员ID(组间的随机效应)看作是一种分类变量放到普通线性回归模型中,但这样作是得不偿失的。有可能这个factor的level很多,可能会用去很多自由度。更重要的是,这样作没什么意义。因为人员ID和性别不一样,我们不清楚它的意义,而且它也不能完整的划分总体。也就是说样本数据中的路人甲,路人乙不能完全代表总体的人员ID。因为它是随机的,我们并不关心它的作用,只是因为它会影响到模型,所以不得不考虑它。因此对于随机效应我们只估计其方差,不估计其回归系数。

混合模型中包括了固定效应和随机效应,而随机效应有两种方式来影响模型,一种是对截距影响,一种是对某个固定效应的斜率影响。前者称为 Random intercept model,后者称为 Random Intercept and Slope Model。Random intercept model的函数结构如下:

Yij = a0 + a1*Xij +  bi + eij

a0: 固定截距
a1: 固定斜率
b: 随机效应(只影响截距)
X: 固定效应
e: 噪声

下在的代码就是人工生成这种数据,再用nlme包建模展示。

参考资料:
Regress by simulation
A very basic tutorial for performing linear mixed effects analyses


如何善用R的语法:三段示例

在改进R代码的运算效率时,我们常说的第一句话就是:要善用R的语法。那么问题是:HOW?下面用一个具体的问题来表达这句比较抽象的话。问题就是用R代码来模拟一个二维布朗运动,一个人从原点出发,随机的向任何一个方向前进或后退一步,不会停留在原地。这个问题是一位英国留学的朋友问我的,我觉得可能比较典型,就顺便写在博客上供大家参考比较吧。

下面就这个问题,给出了三个解答,slowwalk是运行最慢的,在我的电脑上大约要运算8秒,这段代码慢的原因在于用c来补进新的向量元素。而在R中,对一个向量进行新增的操作实际上是进行了复制重建,所以速度是很慢的。为了改进这一点,speedwalk1函数中,先用一个空向量占好了位置,再向这个向量中放东西,这样不需要反复新建,速度会增加10倍,只花了0.7秒左右。但R最重要的运算思路是向量化,也就是同时对向量各元素进行计算,而不是把向量元素一个个放到循环中。speedwalk2函数就是向量化计算的实现,这样下来代码速度又加快了十几倍,大约只需0.05秒,而且还没有用到多核心计算。

奇异值分解和潜在语义分析

奇异值分解是一种矩阵分解技术,如果将列看作特征,SVD也能看作是一种降维技术。如果将数据先进行标准化之后,再进行SVD操作,其左奇异向量U就是主成分得分,右奇异向量V就是主成分负载。但SVD的优势恰恰在于不进行标准化,当我们处理高维稀疏数据(例如文本数据)时,标准化会使矩阵不再稀疏,影响计算效率。所以一般是用SVD直接对原始数据进行操作。得到U的第一列是平均值的意义,U的后几列才有区别各样本的功能。

在文本挖掘中,潜在语义分析(LSA)假设在词项和文档之间有一个潜在语义层,词汇和语义产生关联,文档也和语义产生关联。这样通过语义概念,可以将词项和文档放在一个语义空间中观察分析。进一步我们可以将语义所为词项的“主成分”放到数据挖掘中进行分析。LSA的一般过程是先构造词项-文档矩阵(其转置为文档-词项矩阵,视需要而定),再计算TFIDF值,然后用SVD来分解这个矩阵,得到的奇异向量就是潜在语义。U的第一列可以认为是词项在文档中的平均出现频率,V的第一列可以认为是文档中包含的词项频率,我们更有兴趣的是后几列。关于SVD和LSA有一份精彩的入门文档,下面的R代码就是对该文档的模仿,但是SVD得到的结果却和文档不同,我很奇怪为什么会这样。(使用的数据

说英雄谁是英雄之光荣三国

本文的缘起是知乎的这一篇文章,老马这位三国历史专家让大家都知道了有这么一份三国志游戏数据的存在。这份数据很容易从网上找到,而本人只关心三国志四的各人物数据,为什么呢,因为只玩过三国志四嘛。游戏是多年前玩的,如今来玩玩数据吧。下面的问题主要围绕一个焦点,谁是三国志四中文武兼资第一人。

网上下载的数据是一个Excel表格,整理之后只取了三国志四的部分,去除了分值为0的人物,整理后数据你可以从这里下载。整个数据有454名人物,包括了六个变量,分别为人物’名字’,’统率’,’武力’,’智力’,’政治’,’魅力’。我们可以使用plyr包中的arrange函数找到各指标分别对应的前十名。不出意外,武神关二爷获得统率第一,人中吕布获得武力第一,诸葛亮获得智力第一,张昭获得政治第一,而会哭的刘备获得魅力第一,晕!

arrange(data,desc(tongshuai))[1:10,]
arrange(data,desc(wuli))[1:10,]
arrange(data,desc(zhili))[1:10,]
arrange(data,desc(zhengzhi))[1:10,]
arrange(data,desc(meili))[1:10,]
 
由于每个人物对应多个指标,我们需要一个综合指标来比较,最简单的,上PCA。
prinModel <- prcomp(x=data[,2:6],
retx=TRUE,
center=TRUE,
scale=TRUE)
 
prinModel$rotation[,1:2,drop=F]
summary(prinModel)
 
对数据进行主成份分析,发现前两个主成份的累积方差贡献率达到88%,还不错,第一PC以智力和政治为主,达到50%贡献率,第二PC以武力和统率为主,达到30%贡献率。那么第一主成份可以认为是文韬因素,第二主成份可以认为是武略因素。好的,下面可以根据两个主成份得分列出文韬前十和武略。

# 武略榜
data$name[order(prinModel$x[,2],decreasing=T)][1:10]
[1] "關羽" "趙雲" "曹操" "孫策" "孫堅" "黃忠"
[7] "周瑜" "陸遜" "夏侯惇" "姜維"
# 文韬榜
data$name[order(prinModel$x[,1],decreasing=T)][1:10]
[1] "郭嘉" "張紘" "張昭" "諸葛亮" "荀彧" "諸葛瑾"
[7] "荀攸" "曹植" "魯肅" "馬良"
 
下面要将这两个主成分合为一个综合指标,采用加权平均方法,权重是方差贡献率,大致上0.6*文韬+0.4*武略,下面就取出兼资文武之名将前十名,当当当当,他们是:
final <- prinModel$x[,1]*0.6 + prinModel$x[,2]*0.4
data$name[order(final,decreasing=T)][1:10]
 
[1] "諸葛亮" "曹操" "周瑜" "陸遜" "司馬懿" "趙雲"
[7] "魯肅" "姜維" "龐統" "徐庶"
 
看来在三国志的游戏数据中,诸葛亮是最强人物啊。最后来看看蜀国五虎排序,通常民间排序是关张赵马黄,这种序列显然没有考虑到智力因素问题啊,下面我们用前面的评价体系来为五虎重新排序。
shu5 <- c('趙雲','馬超','張飛','關羽','黃忠')
data.frame(round(prinModel$x[data$name %in% shu5,1:2],2),
final=round(final[data$name %in% shu5],2),
name=data$name[data$name%in% shu5])
 
PC1 PC2 final name
437 -1.51 2.22 -0.02 馬超
488 -2.65 1.31 -1.07 張飛
658 0.59 2.92 1.52 黃忠
741 1.68 3.31 2.34 趙雲
922 1.20 3.48 2.12 關羽
 
从上面这个表格可以看到,根据游戏数据正确的排序应该是赵关黄马张,赵云君,我看好你啊。我们还可以看看五虎将在群雄中的位置,将PC1和PC2做散点图如下,这张图中处于右上角的就是文韬武略兼备的人才。
library(ggplot2)
ggdata <- data.frame(prinModel$x[,1:2],name=data$name)
p <- ggplot(ggdata,aes(PC1,PC2))+
geom_point(size=3,alpha=0.5)+
geom_text(data=subset(ggdata,name %in% shu5),
aes(label=name))+
labs(x='文韬',y='武略')+
theme_bw()
print(p)


 
最后还可以用雷达图看这五虎的技能值


使用PCA来进行综合指标评价不一定是最好的方法,各位有什么好主意,请不吝赐教。

生存分析函数小结

生存分析(survival analysis)适合于处理时间-事件数据。例如中风病人从首次发病到两次复发,其中就涉及到时间和事件。此例中时间就是复发的时间间隔,事件就是是否复发。如果用普通的线性回归对复发时间进行分析,就需要去除那些没有复发的病人样本。如果用Logistic回归对是否复发进行分析,就没有用到时间这个因素。而生存分析同时考虑时间和事情这两个因素,效果会更好些。

在R语言中我们可以使用survival包进行生存分析,其中主要的函数功能罗列如下:

Surv:用于创建生存数据对象
survfit:创建KM生存曲线或是Cox调整生存曲线
survdiff:用于不同组的统计检验
coxph:构建COX回归模型
cox.zph:检验PH假设是否成立
survreg:构建参数模型


下面是使用一个实例来使用R中的生存分析函数,其中用到的数据集可以在这里下载