Author Archive for 写长城的诗

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

奇异值分解是一种矩阵分解技术,如果将列看作特征,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中的生存分析函数,其中用到的数据集可以在这里下载

电影爱好者的R函数

作为一个伪影迷,经常纠结一些电影该不该下,要不要看。毕竟吾生也有涯而片源无涯。还好可以去豆瓣一类的地方看看大家的评分择优录用。去豆瓣查分需要登录网站搜索再鼠标点点点,如果要查好几部电影就有点费事儿。其实可以用R写个函数,先抓取相应的网页,再筛选返回需要的分值。这样在R里头就可以批量查分了,恩,走起来。

library(RCurl)
library(XML)
movieScore <- function(x) {
stopifnot(is.character(x))
# 提交搜索豆瓣表单
search <- getForm("http://movie.douban.com/subject_search", search_text = x)
searchweb <- htmlParse(search)
# 解析搜索结果页面
resnodes <- getNodeSet(searchweb, "//div[@id='wrapper']//table[1]//a")
if (is.null(resnodes))
return(NULL) else resurl <- xmlGetAttr(resnodes[[1]], name = "href")
# 得到影片页面后第二次解析
resweb <- getURL(resurl, .encoding = "UTF-8")
content <- htmlParse(resweb, encoding = "UTF-8")
resnodes <- getNodeSet(content, "//div[@id='interest_sectl']//p[@class='rating_self clearfix']//strong")
namenodes <- getNodeSet(content, "//div[@id='content']//h1//span")
# 得到影片评分
score <- xmlValue(resnodes[[1]])
name <- xmlValue(namenodes[[1]])
return(list(name = name, score = score))
}

看看天机这部大烂片多少分。

movieScore("天机")
## $name
## [1] "天机·富春山居图"
##
## $score
## [1] "2.9"

抓网页比较慢,豆瓣为人民群众着想提供了API,我们也可以使用API来调取分数,函数也比较简单。

library(RCurl)
library(XML)
library(RJSONIO)
movieScoreapi <- function(x) {
api <- "https://api.douban.com/v2/movie/search?q={"
url <- paste(api, x, "}", sep = "")
res <- getURL(url)
reslist <- fromJSON(res)
name <- reslist$subjects[[1]]$title
score <- reslist$subjects[[1]]$rating$average
return(list(name = name, score = score))
}
movieScoreapi("僵尸世界大战")
## $name
## [1] "僵尸世界大战"
##
## $score
## [1] 7.5

有了这个查分函数,我们可以在R中批量查阅电影评分了。但是豆瓣对于频繁的访问会有限制,对于没有认证的API使用是每分钟10次,超过就会暂时封IP。对于网页抓取,肖楠在第六次R会议上有个很棒的演讲,有兴趣的同学可以去统计之都看看。

实现可重复的统计slides

制作幻灯片是数据分析师的必备技能之一,优秀的slides对外可以忽悠住客户,对内可以震慑住领导。精良的slides要秀外慧中,有逻辑有内容,还要有外形有风格。达到这种标准并不容易。在过去,你可能使用office中的ppt来做传统意义的幻灯片,将图片和代码费力的copy到一张张的slides上去,然后到处找模板。但这种方式已经凹特了。

在HTML5的发展背景下,已经出现了大批以网页形式的slides框架,例如:Google IO 2012\HTML5 slides\HTML5 Rocks\Shower\Deck.js这个名单可以很长。而这里只是其中六种演示框架的介绍。在这些slides框架下,你只需要懂一点点web知识,将图片和数据嵌入到一个html模板中,就可以生成一个动态可交互的slides。

实际上对于统计分析的slides,有一种工具可以使我们的制作效率更高。这就是谢益辉的knitr包。knitr像一把钥匙,打开了可重复统计报告的大门,同时也打开了通向html的大门。可以使R代码、运算结果和分析文字自然的融为一体,能自由的输出为LaTex、Markdown或者是Html。在借助knitr打通了这些奇经八脉之后,可以惊奇的发现,knitr可以将以上所有东西粘在一起,有助于实现可重复的统计分析幻灯片。

以后的工作流应该是这样的:数据分析的过程和结果均存放在rmd文件中,可根据不同需要转换为pdf或是html,当需要制作slides的时候,只要将生成的html进行一点加工,加入一些js/css效果即可生成美观炫目的slides。目前已经有很多这类现成的例子,例如coursera的公开课data analysis就是使用了这种slides来讲解R和数据分析。作者利用的是名为slidify的R包,直接从rmd文档生成html5的slides页面,非常推荐各位去下载来研究一下。

slidify使用的默认框架是io2012,比较淡雅朴素,另一种比较炫目的框架是reveal.js,这里有一个例子值得参考下。这份slides的原码可以在此找到。这份slides是使用reports包来生成的,基本流程如下:

install.packages(“reports”)
library(reports)
presentation(‘example’)
setwd(‘~/example’)
reveal.js()

安装reports包后使用presentation建立一个目录,这个目录会在当前R工具目录下自动生成,之后用户在目录下的rmd文件中编辑内容,保存后生成html文件,最后使用reveal.js函数生成slides,你会在reveal.js目录下找到index.html,就是最终的文件。由于reports包的瑕疵问题,在reveal.js前一步你可能需要手动切换一下工作目录。而且最好是再手动编辑一下index.html文件,以完全自定义实现reveal.js的各种特效。笔者也做了一个简单的slides可以看一看。关于reveal.js的完整资料看这里。 

相比reports包,slidify要显得稳定些,使用方法也类似,但它还放在github上,所以要用devtools包来下载它:

install_github(“slidify”, “ramnathv”)
install_github(“slidifyLibraries”, “ramnathv”)
library(slidify)
author(‘example’)
slidify(‘index.Rmd’)

加载slidify后使用author,它会自动在当前工作目录下建立example目录,并打开index.Rmd文档,供用户编辑,当然在写slides的时候需要遵从一定的格式才会让它识别。写完之后使用slidify即可生成网页。目前slidify支持的框架不多,但在开发版本中已经支持了包括reveal.js在内的很多框架了,所以尽请期待吧。