新书推荐：脏数据手册

当你学完一本数据分析软件教程，在电脑上做完了所有的练习题，志得意满地准备去处理实际问题时候，你会被真实世界的“脏数据”所震惊。例如那些随处可见的缺失和格式不一的数据会让分析工作举步维艰，但脏数据的陷阱远不止这些。初入数据江湖的白板青年很需要一本江湖经验手册来帮助成长，而《Bad Data Handbook》正好满足了这种需要。

该书由十九位江湖高手合力联手写成，作者包括了数据分析领域的数据科学家和统计学、物理学、经济学的学术专家，他们分享了自己在处理非常规数据问题的经验和技巧。例如对数据进行分析前的观察测试，如何处理棘手的表格数据，应付各种编码问题，进行网页数据抓取，处理文本数据，以及一些数据质量控制的方法。这些秘诀在一般的传统教科书上是难以见到的。正所谓江湖险恶，要行走数据江湖则必看此书。

来玩玩QQ群的数据

上周COS论坛上有位老兄发布了一个关于QQ群的数据，正好拿来玩玩。这批数据并不复杂，只有两列，一列是用户名，一列是用户发言时间，不过从这批数据中仍然可以得出一些好玩的东西，且让本人一一道来。

推荐两本python书

python虽然不是专门的数据分析工具，但是它的库超多。随着数据分析相关各种库的日益完善，也可以用它来处理一些数据方面的工作，特别在数据预处理方面。这门书号称是数据分析，实际上大部分就是讲的数据处理。介绍了用numpy, pandas等库来实施数据读入、清理、转换、合并等工作。不得不说，pandas的语法真的和R好象啊。不过分析方面没有什么很出彩的地方了。

上面这本是个小册子，专门讲了numpy和scipy两个库，特别还介绍了scikit的一些机器学习用法，值得看一下吧。
连接1，连接2

决策树之三国争霸

决策树是一种简洁实用的数据挖掘方法。在R中通常可以用rpart包和party包来实现两种算法的决策树。最近著名的C4.5决策树算法的升级版本C5.0已经可以在官网下载到。对于这三种决策树算法，本文来做一个预测效果的简单对比。

对比用的数据集是C50包中自带的churn数据，它是用来预测顾客流失的数据集，其中样本量为3333个，变量数为20个。为不平衡数据，没有缺失值存在。对比基本步骤是用10重交叉检验，将数据随机分为10份，用9份训练决策树，用1份来检验结果。循环后求出10个预测准确度的均值。然后在外面再套一个100次大循环，产生三个决策树算法各100个准确率。最后绘制为提琴图，从图中可以观察到C5.0的表现最好，而party次之，rpart的效果最差。在本例实验中最大的差距虽然不过0.02，但如果放在kaggle的数据挖掘比赛中，就相当于是一百位名次的差距了。

如何批量处理文本文件

最近数据堂为了弄数据挖掘比赛提供了一批用户行为日志数据。对于以前没玩过的数据，我是特别的好奇。处理这批文本文件确实花了不少时间。数据以不同的日期作文件夹分别存放，每个文件夹中又有近一千个文本文件，每个文件都是一个用户的行为日志。为了分析这些数据，首先需要将这两万个文本文件读入R中，再用字符串函数进行处理成结构化的数据。处理方法如下：

新书推荐：数据之魅

在amazon书店里头，如果将统计类和数据挖掘类书籍除外的话，还真没有一本正经八百讲数据分析的书。不过《Data Analysis with Open Source Tools》倒是填补了这个空白。一般说到数据分析，可能要么是概念和公式推导居多，要么是软件教程和计算代码居多，这本书的写作风格和内容则大不一样，作者是以一种诚恳的态度在和你聊数据分析的事情，而不是以一种教科书的形式讲技术细节。它的重点在于作者的行业经验和心得。该书提出了很多很好的行业见解，而且知识覆盖范围极广，甚至包括了一些数据建模和BI的东西，不然也不会号称是程序员和数据科学家的指南手册了。

三门问题的模拟

有一个著名的蒙提霍尔问题，亦称为三门问题（英文：Monty Hall problem），大致出自美国的电视游戏节目Let's Make a Deal。问题的名字来自该节目的主持人蒙提·霍尔（Monty Hall）。这个游戏的玩法是：参赛者会看见三扇关闭了的门，其中一扇的后面有一辆汽车，选中后面有车的那扇门就可以赢得该汽车，而另外两扇门后面则各藏有一只山羊。当参赛者选定了一扇门，但未去开启它的时候，知道门后情形的节目主持人会开启剩下两扇门的其中一扇，露出其中一只山羊。主持人其后会问参赛者要不要换另一扇仍然关上的门。问题是：换另一扇门会否增加参赛者赢得汽车的机会率？

这条问题亦被叫做蒙提霍尔悖论：虽然该问题的答案在逻辑上并不自相矛盾，但十分违反直觉。这问题曾引起网络上一阵热烈的讨论。有人认为概率还是1/3，不用换。有人认为概率增大为1/2，可以换。实际上，如果严格按照上述的条件的话，赢得汽车的正确机率是2/3。我们可以用一种极端的方法来考虑，如果有100扇门，主持人打开了998扇空门，你换不换呢？显然是要换的。

还有许多概率推理的方法可以获得正确答案。但我们有计算机和R，让我们来写个小程序模拟一下嘛。

用Shiny包快速搭建基于R的交互网页应用

RStudio是我最喜欢用的R语言IDE，其开发团队最近又推出了一个新的产品，即Shiny包。它的作用是快速搭建基于R的交互网页应用。使得那些对代码不熟悉的人士在工作中也可以应用统计模型。对于R和web的交互，之前已经有一些相关的包，例如：rApache, Rhttpd, Rack, Rook。不过这些工具都需要开发者不仅要熟悉R，还要熟悉网页编程语言（html,CSS,JS）。而Shiny包的特点在于不需要了解网页语言，可以用纯R来搭建。生成的网页应用是动态交互，而且是即时更新的。Shiny还提供了现成组件方便快速在网页上展示数据、图表和模型，的确是非常的炫。本例将用ggplot2包来绘制iris数据集的散点图，并将图形放到网页中。

首先安装Shiny包：
options(repos=c(RStudio='http://rstudio.org/_packages', getOption('repos')))
install.packages('shiny')

参加上海第五届R会议的PPT和代码

上海的这次R会议来了很多知名公司和嘉宾，参与听众反应也非常热烈。感觉比北京的要好一些哦。而我只不过讲了一些业余玩的东东，附上本次演讲的PPT和代码。

果壳中的R第二版新鲜出炉

《R in Nutshell》是O'REILLY公司出版的果壳系列图书之一。该系列图书的特点是知识覆盖面广，讲解全面细致，索引、参考资料以及进一步阅读都包括在内，是非常难得的桌头参考书籍。 《R in Nutshell》也继承了该系列的特点，从简单的R入门知识、R语法、数据整理和可视化，到统计回归、机器学习均有涉及，还包括了代码优化、生物计算和Hadoop的相关内容。

之前统计之都团队正组织翻译出版此书，不过英文第二版正酝酿出来，所以翻译工作暂时停下来。近日在网上看到英文第二版终于出来了，抽时间翻看一下，发现改进不少。特地搬到网盘上供各位下载。

相对于第一版来说，新版本增加了两个完整的章节，即ggplot2绘图和Hadoop。这样使可视化部分终于完整的包括了基本绘图、Lattice包和ggplot2包。而且在大数据背景下，hadoop相关知识的引入也称得上是与时俱进。

新版本还修改完善部分代码，并重新安排了篇章结构。例如原来的第11章high-performance改写后成为了第24章optimizing R progrmas。

此外还增加了一些热门R包的介绍，例如plyr和reshape。在回归模型中则增加了弹性网glmnet。这样总页数从原来第一版的636增加到第二版的722。

总而言之，有此书在手则基本不需要其它的R语言资料了。不过对于新人仍建立先选用较薄的入门小册子。例如《R for beginners》和《R导论》，它们都有中文版。

在R语言中使用SQL

数据分析经常需要从外部获得数据。很多情况下数据存放在关系型数据库中。一般我们可以用SQL来提取需要的数据，存为文本再由R来读入。这种方式结合了数据库的储存能力和R的分析能力，速度也非常快。但是如果要形成一套可重复性的自动工作流程，则可以将R与外部数据库连接，直接在R中操作数据库，并生成最终结果，这也是一种可行的方法。

在R中连接数据库需要安装其它的扩展包，根据连接方式不同我们有两种选择：一种是ODBC方式，需要安装RODBC包并安装ODBC驱动。另一种是DBI方式，可以根据已经安装的数据库类型来安装相应的驱动。因为后者保留了各数据库原本的特性，所以个人比较偏好用DBI连接方式。有下面这几种主要的包提供了DBI连接：RMySQL，RSQLite，ROracle，RPostgreSQL。由名字看得出它们分别对应了几种主流的数据库。

Economist风格的统计绘图

《Economist》（经济学人）是一份由伦敦经济学人报纸有限公司出版的杂志，于1843年9月由詹姆士·威尔逊创办。杂文章写得机智，幽默，有力度，严肃又不失诙谐，并且注重于如何在最小的篇幅内告诉读者最多的信息。杂志主要关注政治和商业方面的新闻，但是每期也有一两篇针对科技和艺术的报导，以及一些书评。从2012年1月28日的那一期杂志开始，《经济学人》杂志开辟了中国专栏，为有关中国的文章提供更多的版面。

杂志不仅文章内容精彩，排版和图片也非常具有特色。其中的统计图片均是以蓝色为基调，简洁而大气。如果你想绘制出这种风格的统计图，现在也很简单了。只需要使用R语言中大名鼎鼎的ggplot2包，以及Jeffrey Arnold新开发的ggthemes包，就能轻易的实现了。

ggthemes包还不能从CRAN库中安装，所以需要从github上下载原始文档，编译成zip本地安装包，最后在R中进行安装加载就可以使用了。如果你已经安装了Rtools，也可以使用devtools包中的install_github()函数来直接安装github上的R包。

下面就是用ggplot2+ggthemes包作的一个小例子，其中重要的就是要用到theme_economist()和scale_color_economist()这两个函数来实现economist风格的统计图。很简单吧！ggthemes包中还有其它主题可供使用，慢慢研究。

p <- ggplot(data=mpg, mapping=aes(x=cty,y=hwy))
p + geom_point(aes(color=factor(year)),
               alpha=0.5,position = "jitter")+
    stat_smooth()+
    theme_economist() + scale_color_economist()+
    labs(title = '汽车型号与油耗',
         y = '每加仑高速公路行驶距离',
         x = '每加仑城市公路行驶距离',
         colour = '年份')

用Parallel和foreach包玩转并行计算

众所周知，在大数据时代R语言有两个弱项，其中一个就是只能使用单线程计算。但是在2.14版本之后，R就内置了parallel包，强化了R的并行计算能力。parallel包实际上整合了之前已经比较成熟的snow包和multicore包。前者已经在之前的文章中介绍过了，而后者无法在windows下运行，所以也就先不管了。parallel包可以很容易的在计算集群上实施并行计算，在多个CPU核心的单机上，也能发挥并行计算的功能。我们今天就来探索一下parallel包在多核心单机上的使用。

parallel包的思路和lapply函数很相似，都是将输入数据分割、计算、整合结果。只不过并行计算是用到了不同的cpu来运算。下面的例子是解决欧拉问题的第14个问题。

使用GitHub进行版本控制的傻瓜方法

不论是团队合作还是单打独斗，代码和文档的版本控制是数据极客不可缺少的工具。高阶极客能随心所欲的用 linux终端+Git+编辑器完成这类任务。但是对于像本人一样的Git入门者来讲，图形工具还是略微让人心安一点。我们下面就来示例，用RStudio结合GitHub for Windows来完成这项任务。

首先你需要在GitHub中建立一个帐号，然后安装上述两种软件。试着登录GitHub for Windows看是否正常。此时你的本地库和远程库应该都是空的。之后在RStuio中点击新建项目，以后这个项目中的代码都会同步到远程代码库中。

抓取网页数据的几种套路

没有数据就没有乐趣。有的数据提供者心肠很好，会直接给出txt或是csv文档。这个时候我们可以直接在R里头用read.table()函数把数据读进来。

有的时候我们需要的数据在网页上以一个表格呈现，例如前面文章遇到过的地震数据。此时可以用XML包中的readHTMLTable()函数读取数据，后续再配合一些字符串处理一般就OK了。

如果你对R不大熟悉，抓取这些表格也有更方便的法子，就是利用Chrome的扩展。有两个扩展值得推荐使用：一个扩展叫作table capture，它会自动找出网页中的若干表格，你只需选择所需的那个将其拷贝到剪贴板即可，然后再用下面的命令就可以读入到R中。
data <- read.table('clipboard',T)

另一个扩展叫作scraper。先选择你所需要的部分内容，然后右键选择scraper similar也能抓取表格，不过它会存到一个google doc中去。在天朝这玩意儿不大方便。

有些数据不是以表格方式出现的，例如用XML或是JSON方式储存的数据。在R中都有对应的包来处理。下面的示例即是用XML包来处理XML数据。在此之先你需要有一点关于XML和XPath的知识，首先处理的对象是这样一个页面：http://www.w3schools.com/xml/plant_catalog.xml

library(XML)
xml.url <- "http://www.w3schools.com/xml/plant_catalog.xml"
# 解析xml页面
xmlfile <- xmlTreeParse(xml.url)
# 观察对象属性
class(xmlfile)
# 获取根结点
xmltop <- xmlRoot(xmlfile)
# 用xmlValue函数获取叶结点处的值
xmlValue(xmltop[[1]][[1]])
xmlValue(xmltop[['PLANT']][['COMMON']])
# xmlSApply类似于sapply函数，取出第一个子结点中的所有叶结点值
xmlSApply(xmltop[[1]],xmlValue)
# 进一步可以取出所有子结点中的叶结点值
plantcat <- xmlSApply(xmltop, function(x) xmlSApply(x, xmlValue))
# 将数据转为数据框
plantcat_df <- data.frame(t(plantcat),row.names=NULL)
plantcat_df[1:5,1:4]

有时候会遇到更为复杂的XML页面，此时的节点内含有参数值。如果要获取这些数据则需要使用getNodeSet()函数配合xmlValue()函数。当遇到更为复杂的数据，那我们只能用readLines读进来，再用字符串函数配合正则表达式来加以处理了。

参考资料：

http://www.omegahat.org/RSXML/Tour.pdf
http://www.stat.berkeley.edu/~statcur/Workshop2/Presentations/XML.pdf

如何在WIN下写一个简单的R包

虽然玩了一段时间的R，但很惭愧的是一直没有学着自己编写一个R包。一个是觉得R本身的包已经是浩如烟海了，另一个感觉好象写包非常麻烦。在参加完北京的这次统计夏令营后，看太云讲写包也不是太难。于是回来找了些资料尝试了一下，于是就形成了这篇笔记。

在WIN下面除了要安装好R之外，还要安装Rtools，这个在官网上可以找到。如果不写复杂的帮助文档和Vignette的话也不需要装MikTeX。此外在WIN中还要确保PATH的设置正确，这些可以在后面的参考文献中得到详细说明。在编写R包时还要知道自己的工作目录在哪，可以用getwd()函数来了解。

笨办法学R编程（6）

有时候用R来解一些Project Euler的题目会非常简单，今天就来三题连解（6、7、8）。题目就不再这里复述了，可以查看官方网站。用到函数和表达式大部分在前面都已经熟悉过了，不过还是会接触到一些新的函数。废话不多说直接上代码吧。

# Euler 6  
x <- 1:100
sum(x)^2 - sum(x^2)
 
# Euler 7 
 
n <- 0
i <- 1
m <- rep(0,10001)
while (n <10001) {
    if (findprime(i)) {
        n <- n +1 
        m[n] <- i}
    i <- i + 1
}
m[10001]
 
# 预备练习，熟悉一些字符串操作函数
text <- c('hello','world','I','love','code')
gsub('o',' ',text)
gsub('o','*',text)
gsub('o','',text)
 
(temp1 <-paste(text,collapse=' '))
paste(text,collapse='*')
paste(text,collapse='')
 
(temp2 <- strsplit(temp1," "))
class(temp2)
(temp3 <- unlist(strsplit(temp1," ")))
class(temp3)
 
# Euler 8 
 
web <- 'http://projecteuler.net/problem=8'
# 用readLines函数来抓取网页
raw <- readLines(web)
raw <- raw[54:73]
# 删除多余字符串
data <- gsub('<br />','',raw)
# 粘合成一个字符串
num <- paste(data,collapse='')
# 分割后转为数值向量
temp <- as.numeric(unlist(strsplit(num,'')))
 
p <- numeric()
for ( i in 1:(1000-4)) {
    p[i] <- prod(temp[i:(i+4)])
}
max(p)

第六题用R的向量化计算非常简单，第七题需要用到第二题中建立的findprime函数。第八题的解决是用字符串方法，先将那一长串数字作为字符串切开，再转为数值型向量，最后用循环求乘积。为了偷懒，是直接抓取的网页，没有输入那个长长的数字。

EM算法的R实现和高斯混合模型

EM（Expectatioin-Maximalization）算法即期望最大算法，被誉为是数据挖掘的十大算法之一。它是在概率模型中寻找参数最大似然估计的算法，其中概率模型依赖于无法观测到的隐变量。最大期望算法经过两个步骤交替进行计算，第一步是计算期望（E），也就是将隐藏变量象能够观测到的一样包含在内，从而计算最大似然的期望值；另外一步是最大化（M），也就是最大化在 E 步上找到的最大似然的期望值从而计算参数的最大似然估计。M 步上找到的参数然后用于另外一个 E 步计算，这个过程不断交替进行。对于信息缺失的数据来说，EM算法是一种极有效的工具，我们先用一个简单例子来理解EM算法，并将其应用到GMM（高斯混合模型）中去。

幼儿园里老师给a，b，c，d四个小朋友发糖吃，但老师有点偏心，不同小朋友得到糖的概率不同，p(a)=0.5, p(b)=miu, p(c)=2*miu, p(d)=0.5-3*miu 如果确定了参数miu，概率分布就知道了。我们可以通过观察样本数据来推测参数。设四人实际得到糖果数目为(a,b,c,d)，可以用ML（极大似然），估计出miu=(b+c)/6*(b+c+d)，假如某一天四个小朋友分别得到了(40,18,0,23)个糖。根据公式可求出miu为0.1。在信息完整条件下，ML方法是很容易估计参数的。

假设情况再复杂一点：知道c和d二人得到的糖果数，也知道a与b二人的糖果数之和为h，如何来估计出参数miu呢？前面我们知道了，如果观察到a，b，c，d就可以用ML估计出miu。反之，如果miu已知，根据概率期望 a/b=0.5/miu，又有a+b=h。由两个式子可得到 a=0.5*h/(0.5+miu)和b=miu*h/(0.5+miu)。此时我们面临一个两难的境地，a，b需要miu才能求出来，miu需要a，b才能求出来。有点类似岸上的八戒和河里的沙僧的对话：你先上来，你先下来，你先上来......

那么一种思路就是一方先让一步，暂且先抛出一个随机的初值，然后用对方算出的数值反复迭代计算。直到计算结果收敛为止。这就是EM算法的思路，我们来看看用R来实现这种思路：

笨办法学R编程（5）

随着教程推进，基本的语法都接触得差不多了。当要解决某个具体问题时，只需要考虑用什么样的算法来整合运用这些函数和表达式。今天来解决Project Euler的第五个问题，该问题可以用很笨的暴力搜索法子来作，但是更聪明的作法是采用质因子分解的思路。即任何一个合数都可以分解为质数的乘积。为了完成这个题目，还需要学习一点点矩阵，以及和sapply函数相似的另一个函数apply。

# 预备练习
mat <- matrix(1:12,ncol=4)
print(mat)
t(mat)
colnames(mat) <- c('one','two','three','four')
rownames(mat) <- c('a','b','c')
print(mat)
apply(mat,1,sum)
apply(mat,2,sum)
sum(apply(mat,2,sum))
prod(apply(mat,2,sum))

笨办法学R编程（4）

看到各位对“笨办法系列”的东西还比较感兴趣，我也很乐意继续写下去。今天的示例将会用到数据框（data.frame）这种数据类型，并学习如何组合计算两个向量，以及如何排序。我们将用所学的东西来解决Project Euler的第四个问题，就是找出一个集合中最大的回文数。回文数是指一个像1534351这样“对称”的数，如果将这个数的数字按相反的顺序重新排列后，所得到的数和原来的数一样。开始啦！

# 预备练习
 
x <- y <- 1:9
data <- expand.grid(x=x,y=y)
print(data)
z <- data$x * data$y
# 一个九九乘法表
z <- matrix(z,ncol=9)
 
set.seed(1)
x <- round(runif(10),2)
print(x)
order(x)
x[order(x)[1]]
which.min(x)
x[which.min(x)]
x[order(x)]

笨办法学R编程（3）

经历了前面两个小挑战，你应该对R有点理解了。我们继续推进，今天的问题有点点复杂，复杂的不是R，而是一个数学概念：质数和质因子。任何一个合数都可以被几个质数所分解，这个性质很重要，我们将用它来解决Project Euler的第三个问题。还是和之前一样的，你需要自己在R控制台中敲打下面这些命令，根据结果自行揣摩其用处。

# 预备练习，学习for循环、建立自定义函数和其它一些函数
 
for (n in 1:10) {
    print(sqrt(n))
}
 
x <- c('hello','world','I','love','R')
for (n in x) {
    print(n)
}

新书推荐：数据新闻手册

如果你和我一样是数据爱好者，那么一定会经常造访卫报的数据博客栏目，去看看他们是如何用不同来源的数据制作有趣的新闻。你一定想知道，他们是如何提出问题的？如何采集数据？又如何整理和分析数据，并最终用可视化方法来展现到读者面前的。答案就在《数据新闻手册》这本书中。

笨办法学R编程（2）

本例将介绍R语言中的while循环和if条件。最终用它来解决Project Euler的第二个问题。除了练习之外你还需要了解一些斐波纳契数列的知识。废话不多说了，打开R控制台，跟着输入下面的代码，自行琢磨吧。

# 预备练习，while循环和if判断
x <- 1:10
print(x)
print(x[10])
print(x[-10])
 
i <- 1
while (i <= 10) {
  print(x[i])
  i <- i + 1
}

笨办法学R编程（1）

在倚天屠龙记中，有一人唤作火工头陀。此人练功不靠心法，只靠模仿他人招式，由外而内，自成一家。练习编程也有如此的法门，不看文字描述，只观察和模仿别人的代码。这样也可以由外而内学会编程。《笨办法学python》的作者Zed Shaw 就说过这种笨办法入门其实更简单。阳志平在他的文章《如何学习一门新的编程语言》中也讲到，初学编程要在学习区刻意的大量练习，少看理论书。

TED上一位教育家同样谈到这么一个故事，他把一个计算机扔在一个偏远的印度小村子里不去管它，在那里没有上过学的小孩就能自己学会英语和计算机的用法。实际上人脑是非常善于自我探索和学习的。因此本系列教程的特点就是只有演示代码加少量注释。通过反复模仿和练习，揣摩代码的变化和结果，你就能自行领悟其含义，并打下坚实的编程基础。

本系列每篇文章的目的都是用R语言编程来解决一个Project Euler的问题。Project Euler是一系列由易到难的计算机编程挑战，它提供了一个平台来激发我们解决问题的灵感和思路。本人写这个教程的目的有三：一是为了好玩，二是提高编程水平，三是示范说明以提供给需要的R初学者。另外从R-Blogger上了解，已经有两位高人用R在计算Project Euler，各位也可以参照他们的文章（博客1、博客2）。

Let's Go

2012伦敦奥运的金牌项目分布

伦敦奥运会尘埃落定，又到了数金牌的时候了。上图还是用的R中的treemap包绘制了这个大game中各项目的金牌分布情况，矩形的大小反映了项目的重要性。比如说在所有302块金牌中，Athletics（田径）所占的金牌比例就很高。与之相当的则是Aquatics（水上项目）。而中国所获金牌数则用不同的颜色深度来表示，例如说在Badminton（羽毛球）上拿到了5枚金牌，颜色就非常深。可以观察到在很多项目中，中国仍是一片空白。

参考资料：
http://www.london2012.com/country/china/medals/index.html
http://en.wikipedia.org/wiki/2012_Summer_Olympics#Sports
整理后的数据文档

代码：

在R语言中绘制Treemap

Treemap是一种流行的可视化技术，它用不同尺寸的嵌套矩形来表现层次数据。最常见的层次数据包括了文件目录、文章结构、组织结构等等。Treemap技术是通过矩形的大小来显示节点的重要性，并通过嵌套层次来显示结构的层级。R语言中的treemap包就可以实现这种可视化技术，下图就是对苹果公司财务报表的可视化。

中国的环境状况在全球的位置

由美国耶鲁大学和哥伦比亚大学联合推出的“年度全球环境绩效指数”（EPI）排名在2012年年初放出。此排名旨在评估一个国家的环境政策，环境卫生与生态系统之平衡的状态，涵盖10项领域共22项环境指标，涉及领域包括气候变化、农业、渔业、森林、水源、空气污染及环境负担等。详细报告和数据都可以在网站主页找到。本文就使用这个数据集来观察中国的环境状况在全球的位置。

上面的小提琴图即是22项指标的分布，叠加的红色圆点表示了中国得分位置。这些指标已经被研究人员进行了转换，所有的指标均是得分越高越好。从这些得分点可以看出，中国在森林保护(FORGROINV)和生物多样性(PACOV)方面做得还不错，CO2排放(CO2GDP)和饮用水(WATSUPINV)等方面做得一般，在PM2.5和室内空气(INDOOR)等方面处于垫底情况。更多的指标说明请参看报告主页。其中的EPI是环境绩效评价的综合得分。

用igraph包探索世界航空网络

本文使用的数据仍然是上篇博文中用到的世界航班数据，不过本例不再仅限于中国国内航班。如果用社交网络的角度来观察数据，一个机场可以看作是一个人，而机场之间的来往航班可以看作是人与人之间的某种联系。整体世界的航线可以看作是一个社交网络。那么用R语言的igraph包来简单探索一下这个社交网络，看能不能得到什么发现。现在星图真得很热门，所以本文最后也会搞一个山寨出来。

中国国内航线信息的可视化

上图是对国内机场和航线信息进行了一个简单的可视化。圆点表示了中国163个机场的位置，线条显示了5381条航线。之前曾在这个网站上见到了作者用R语言来对全世界的航线进行可视化。正所谓见贤思齐，本图就是模仿山寨的结果。但是这个图的生成没有原文那么复杂，所用到的地理图形包和步骤也与原例略有不同，比较失败的是没有展现出原图的夜景效果。具体实施的步骤如下：

• 从这个网站下载到机场数据和航线数据；
• 从中挑选出中国的机场和国内航线，并加以整理；
• 用ggmap包读取谷歌地图；
• 将机场和航线信息绘制在地图上。

用stringr包处理字符串

《Machine Learning for Hackers》一书的合著者John Myles White近日接受了一个访谈。在访谈中他提到了自己在R中常用的几个扩展包，其中包括用ggplot2包来绘图，用glmnet包做回归，用tm包进行文本挖掘，用plyr、reshape、lubridate和stringr包进行数据预处理。这些包本博客大部分都有所介绍，今天就来看看这个遗漏的stringr包。

从名字就看得出，stringr包是用来处理字符串的。R语言本身的字符处理能力已经不错了，但使用起来并不是很方便。stringr包将原本的字符处理函数进行了打包，统一了函数名和参数。在增强功能基础上，还能处理向量化数据并兼容非字符数据。stringr包号称能让处理字符的时间减少95%。下面将其中的一些主要函数罗列一下。

library(stringr)
 
# 合并字符串
fruit <- c("apple", "banana", "pear", "pinapple")
res <- str_c(1:4,fruit,sep=' ',collapse=' ')
str_c('I want to buy ',res,collapse=' ')
 
# 计算字符串长度
str_length(c("i", "like", "programming R", 123,res))
 

来玩一玩全球500强排行榜数据

金融时报在7月20日公布了全球500强排行榜。根据这个数据尝试回答下面的一些问题。

1. 哪个行业的上榜公司最多？

看得出来，银行、石油、制药是前三强。

用gbm包实现随机梯度提升算法

中国有句老话：三个臭皮匠，顶个诸葛亮。这个说法至少在变形金刚中得到了体现，没有组合之前的大力神只是五个可以被柱子哥随手秒掉工地苦力。但组合之后却是威力大增。在机器学习领域也是如此，一堆能力一般的“弱学习器”也能组合成一个“强学习器”。前篇文章提到的随机森林就是一种组合学习的方法，本文要说的是另一类组合金刚：提升方法(Boosting)。提升方法是一大类集成分类学习的统称。它用不同的权重将基学习器进行线性组合，使表现优秀的学习器得到重用。在R语言中gbm包就是用来实现一般提升方法的扩展包。根据基学习器、损失函数和优化方法的不同，提升方法也有各种不同的形式。

自适应提升方法AdaBoost
它是一种传统而重要的Boost算法，在学习时为每一个样本赋上一个权重，初始时各样本权重一样。在每一步训练后，增加错误学习样本的权重，这使得某些样本的重要性凸显出来，在进行了N次迭代后，将会得到N个简单的学习器。最后将它们组合起来得到一个最终的模型。

梯度提升方法Gradient Boosting
梯度提升算法初看起来不是很好理解，但我们和线性回归加以类比就容易了。回忆一下线性回归是希望找到一组参数使得残差最小化。如果只用一次项来解释二次曲线一定会有大量残差留下来，此时就可以用二次项来继续解释残差，所以可在模型中加入这个二次项。

同样的，梯度提升是先根据初始模型计算伪残差，之后建立一个基学习器来解释伪残差，该基学习器是在梯度方向上减少残差。再将基学习器乘上权重系数(学习速率)和原来的模型进行线性组合形成新的模型。这样反复迭代就可以找到一个使损失函数的期望达到最小的模型。在训练基学习器时可以使用再抽样方法，此时就称之为随机梯度提升算法stochastic gradient boosting。

随机森林及其副产品

随机森林(Random Forest)方法是Leo Breiman于2001年提出的一种集成学习（Ensemble Learning）方法，它是传统决策树方法的扩展，将多个决策树进行组合，来提高预测精度。随机森林利用分类回归树(CART)作为其基本组成单元，也可称之为基学习器或是子模型。CART在之前的文章中我们已经介绍过，就不再详细说明了。而集成学习的思路是试图通过连续调用单个学习算法，获得不同的学习器，然后根据规则组合这些学习器来解决同一个问题，可以显著的提高学习系统的泛化能力。组合多个学习器主要采用加权平均或投票的方法。常见的集成学习算法还包括了装袋算法（Bagging)和提升算法（Boosting）。

1. 随机森林计算步骤

• 从原始训练样本中随机有放回抽出N个样本；
• 从解释变量中随机抽出M个变量；
• 依据上述得到的子集实施CART方法（无需剪枝），从而形成一个单独的决策树；
• 重复上面步骤X次，就构建了有X棵树的随机森林模型。
• 在对新数据进行预测分类时，由X棵树分别预测，以投票方式综合最终结果。

O'Reilly新书推荐：用R和Ruby来探索万物

这本书绝对是个另类，它并不以严肃的学术研究或商业项目作为主题，而是以好玩为宗旨。用R和Ruby这两种免费工具，来探索我们身边的各种数据资源。

首先作者用两章篇幅对Ruby和R作了一个介绍。之后第三章来解决办公室内的洗手间数量问题，用Ruby来模拟人们上洗手间的次数，然后用R来绘制各种可能的情况。第四章建立了一个简单的经济动态系统，其中包括了生产者、消费者、价格和市场，并对这些因素进行了模拟。第五章比较有趣，作者用Ruby中的mail库获取了安然丑闻中的电邮数据，然后用R进行了电邮的时间分布描述和文本挖掘。最有Geek味道的是第六章，作者自制声频拾取器测量自己的心跳，并用Ruby来处理音频文件，最后用R绘制出了心跳波形，还分析了自己的心率数据。最后两章同样是模拟，简单模拟分析了生物迁徙和人类社会的进化。

这本书立意非常新奇有趣，但称不上非常有Discovery的感觉，因为它并非用大量的篇幅来介绍如何挖掘真实世界的数据，很多章节内容是用Ruby来进行动态模拟，然后用R中的ggplot2包来可视化展现。不过其中的第五章和第六章还是很精彩的。下载本书电子版

此外，如果想学习用R来抓取真实数据进行分析，我建议看这个小册子Data_Mashups_in_R

在R中进行基于稳健马氏距离的异常检验

我们研究的数据中经常包含着一些不同寻常的样本，这称之为异常值(Outlier)。这些异常值会极大的影响回归或分类的效果。异常值产生的原因有很多，其中可能是人为错误、数据测量误差，或者是实际确实存在这样的异常。为了使模型能够反映大部分数据的规律，所以在数据预处理阶段要进行异常值检测，为下一步分析奠定基础。还有一类情况是，当研究人员希望发现不平凡的事物时，异常值检测本身就是分析的首要目的。例如在信用卡欺诈、计算机入侵检测等问题中。此时由于样本的不平衡性，导致一般的分类方法无法使用，必须转而考虑异常检测方法。

一种常用的异常检验思路是观察各样本点到样本中心的距离。如果某些样本点的距离太大，就可以判断是异常值。这里距离的度量一般使用马氏距离(Mahalanobis Distance)。因为马氏距离不受量纲的影响，而且在多元条件下，马氏距离还考虑了变量之间的相关性，这使得它优于欧氏距离。

但是传统的马氏距离检测方法是不稳定的，因为个别异常值会把均值向量和协方差矩阵向自己方向吸引，这样算出来的样本马氏距离起不了检测异常值的所用。所以首先要利用迭代的思想构造一个稳健的均值和协方差矩阵估计量，然后计算稳健马氏距离(Robust Mahalanobis Distance)。这样使得异常值能够正确地被识别出来。

在mvoutlier包中提供了基于稳健马氏距离的异常值检验方法。我们首先构造一个二维变量的人工数据，其中80个样本是标准正态分布，另一小撮别有用心的样本是均值为5，标准差为1的观测值。我们首先使用uni.plot函数在一维空间中观察这个数据。

谈一谈支持向量机分类器

支持向量机（Support Vector Machine）名字听起来很炫，功能也很炫，但公式理解起来常有眩晕感。所以本文尝试不用一个公式来说明SVM的原理，以保证不吓跑一个读者。理解SVM有四个关键名词：分离超平面、最大边缘超平面、软边缘、核函数。

• 分离超平面（separating hyperplane）：处理分类问题的时候需要一个决策边界，好象楚河汉界一样，在界这边我们判别A，在界那边我们判别B。这种决策边界将两类事物相分离，而线性的决策边界就是分离超平面。
• 最大边缘超平面（Maximal Margin Hyperplane）：分离超平面可以有很多个，怎么找最好的那个呢，SVM的作法是找一个“最中间”的。换句话说，就是这个平面要尽量和两边保持距离，以留足余量，减小泛化误差，保证稳健性。或者用中国人的话讲叫做“执中”。以江河为国界的时候，就是以航道中心线为界，这个就是最大边缘超平面的体现。在数学上找到这个最大边缘超平面的方法是一个二次规划问题。
• 软边缘（Soft Margin)：但世界上没这么美的事，很多情况下都是“你中有我，我中有你”的混杂状态。不大可能用一个平面完美的分离两个类别。在线性不可分情况下就要考虑软边缘了。软边缘可以破例允许个别样本跑到其它类别的地盘上去。但要使用参数来权衡两端，一个是要保持最大边缘的分离，另一个要使这种破例不能太离谱。这种参数就是对错误分类的惩罚程度C。
• 核函数(Kernel Function)，为了解决完美分离的问题，SVM还提出一种思路，就是将原始数据映射到高维空间中去，直觉上可以感觉高维空间中的数据变的稀疏，有利于“分清敌我”。那么映射的方法就是使用“核函数”。如果这种“核技术”选择得当，高维空间中的数据就变得容易线性分离了。而且可以证明，总是存在一种核函数能将数据集映射成可分离的高维数据。看到这里各位不要过于兴奋，映射到高维空间中并非是有百利而无一害的。维数过高的害处就是会出现过度拟合。

所以选择合适的核函数以及软边缘参数C就是训练SVM的重要因素。一般来讲，核函数越复杂，模型越偏向于拟合过度。在参数C方面，它可以看作是LASSO算法中的lambda的倒数，C越大模型越偏向于拟合过度，反之则拟合不足。实际问题中怎么选呢？用人类最古老的办法，试错。

朴素贝叶斯分类与贝叶斯网络

朴素贝叶斯分类（Naive Bayes Classifier）是一种简单而容易理解的分类方法，看起来很Naive，但用起来却很有效。其原理就是贝叶斯定理，从数据中得到新的信息，然后对先验概率进行更新，从而得到后验概率。好比说我们判断一个人的品质好坏，对于陌生人我们对他的判断是五五开，如果说他做了一件好事，那么这个新的信息使我们判断他是好人的概率增加了。朴素贝叶斯分类的优势在于不怕噪声和无关变量，其Naive之处在于它假设各特征属性是无关的。而贝叶斯网络（Bayesian Network）则放宽了变量无关的假设，将贝叶斯原理和图论相结合，建立起一种基于概率推理的数学模型,对于解决复杂的不确定性和关联性问题有很强的优势。下面我们用mlbench包中的一个数据集来看看如何用这两种方法进行学习训练。

PimaIndiansDiabetes2是美国一个疾病研究机构所拥有的一个数据集，其中包括了9个变量，共有768个样本。响应变量即是对糖尿病的判断，它是一个二元变量。其它各解释变量是个体的若干特征，如年龄和其它医学指标，均为数值变量。其中有一些缺失值的存在，虽然朴素贝叶斯分类对于缺失值并不敏感，我们还是先将其进行插补，再进行建模以观察准确率。R语言中的e1071包中就有可以实施朴素贝叶斯分类的函数，但在本例我们使用klaR包中的NaiveBayes函数，因为该函数较之前者增加了两个功能，一个是可以输入先验概率，另一个是在正态分布基础上增加了核平滑密度函数。为了避免过度拟合，在训练时还要将数据分割进行多重检验，所以我们还使用了caret包的一些函数进行配合。

# 加载扩展包和数据
library(caret)
data(PimaIndiansDiabetes2,package='mlbench')
# 对缺失值使用装袋方法进行插补
preproc <- preProcess(PimaIndiansDiabetes2[-9],method="bagImpute")
data <- predict(preproc,PimaIndiansDiabetes2[-9])
data$Class <- PimaIndiansDiabetes2[,9]
# 使用朴素贝叶斯建模，这里使用了三次10折交叉检验得到30个结果
fitControl <- trainControl(method = "repeatedcv", number = 10, repeats = 3,returnResamp = "all")
model1 <- train(Class~., data=data,method='nb',trControl = fitControl,tuneGrid = data.frame(.fL=1,.usekernel=F))
# 观察30次检验结果，发现准确率在0.75左右
resampleHist(model1)
# 返回训练数据的混淆矩阵
pre <- predict(model1)
confusionMatrix(pre,data$Class)

用lubridate包来处理时间数据

人生有一道难题，那就是如何使一寸光阴等于一寸生命。在数据分析中也有一道难题，那就是如何自如的操作时间数据。R语言的基础包中提供了两种类型的时间数据，一类是Date日期数据，它不包括时间和时区信息，另一类是POSIXct/POSIXlt类型数据，其中包括了日期、时间和时区信息。一般来讲，R语言中建立时序数据是通过字符型转化而来，但由于时序数据形式多样，而且R中存贮格式也是五花八门，例如Date/ts/xts/zoo/tis/fts等等。用户很容易被一系列的数据格式所迷惑，所以时序数据的转化和操作并不是非常方便。所幸的是，我们有了lubridate包。lubridate包主要有两类函数，一类是处理时点数据（time instants），另一类是处理时段数据（time spans）。

时点类函数，它包括了解析、抽取、修改。

# 从字符型数据解析时间，会自动识别各种分隔符
> x <- ymd('2010-04-08')
# 观察x日期是一年中的第几天
> yday(x)
# 修改x日期中的月份为5月
> month(x) <- 5

用ggmap包进行地震数据的可视化

最近又发现了一个比较好玩的包ggmap。从名字上可以猜测出来，它的作用就是将ggplot2和map相结合。这样R语言用户能方便的获取各种静态地图数据，并在其基础上使用强大的ggplot绘图工具。ggmap包整合了四种地图资源，分别是Google、OpenStreetMaps、Stamen和Cloudmade。为了演示ggmap的作用，本例是从地震信息网获取最近一周的地震数据，得到其经纬度，然后以散点形式绘制在google地图上，另外也显示地震发生的密度估计。这个思路本质上和之前的一篇博文是一致的，但用ggmap包实现起来更为简单。

如何用R来处理图片

做为“会电脑”的人，除了“友情”帮别人装系统杀杀毒之外，时常会承担一些图片处理的活。也就是对一些照片施加缩放、旋转、裁剪之类的事情。这类小事自然无需动用photoshop这种庞然大物了。在R语言中就有一个很好玩的biOps包，它可以很方便的进行各类图片处理操作。而且用程序来处理图片有一个天然的优势，就是很容易进行批量处理和自动化处理。下面以本人的头像机器人瓦利来做一些简单的示范。

# 加载扩展包
library(biOps) 
# 读取本地的图像文件，观察到此图像是399*399像素的rgb图片
x <-  readJpeg("d:\\xccds.jpg")
print(x)
# 先尝试缩放操作，后面的参数是缩放的比例，可以采用四种方式进行插值，这里用的是最近邻法
plot(imgScale(x,1.3,1.3,interpolation='nearestneighbor'))
#之后进行柔化降噪，此处采用的是中位数滤镜
plot(imgBlockMedianFilter(x,5))

ggplot的图形组合与添加jpeg文档的方法

最近有位朋友问了一个关于ggplot2作图的问题。涉及到多个图形组合的问题，所以还是费了一些时间来解决，自己也从中学习了一些新东西。顺手将这个画图的过程扔上来当作一篇博文吧。下图就是最后的结果。画这个图有几个障碍，一个是二维散点的置信椭圆，另一个是一维直方图的边缘显示。解决的方法是用ellipse包来生成置信椭圆数据，并且用gridExtra包来组合多个图形，还要用opts来去除图例标注的显示。

初次尝试igraph包

igraph是为了进行社会网络分析而创建的一个包。与R语言中同类包相比，它的速度更快，而且函数命令与图形展现更为丰富。它可以处理有向网络和无向网络，但无法处理混合网络。igraph中的函数非常多，本文只是初步的介绍如何创建图形，并提供一些简单的例子。

# 可以使用最基本的graph函数，用向量作为参数来创建图形，之后用plot绘制出结果
library(igraph)
g1 <- graph( c(0,1, 1,2, 2,3, 3,4))
plot(g1,layout=layout.circle(g1))

# 也可以画出一些特殊结构的图形，例如下面的星形图
g2 <- graph.star(10, mode = "in")
plot(g2,layout=layout.fruchterman.reingold(g2))

基于模型的决策树

线性模型是很常用的一类建模方法，它可以用于回归和分类问题。但有时候用全体数据来建模并不理想，这是因为数据内部存在一些局部特性。例如我们要考虑教育对收入的影响程度，那么对于不同的性别，这个影响程度可能也是不同的。所以我们会将数据先进行分割，然后对每个组分别进行回归建模。这种模型和分割相结合的方法又称为基于模型的递归分割算法（model-based recursive partitioning algorithm）。这个工作可以通过手工分组或加入虚拟变量来做，更方便的是利用party包中的mob函数来完成。
mob函数算法的基本步骤如下：

• 对所有数据构建线性模型
• 根据不同的分割变量，评价模型参数的不稳定性。选择某个分割变量使模型参数不稳定性最大。
• 计算分割变量的分割点
• 将数据根据上述结果分成两个子集
• 重复上面的步骤

下面的算例来自于faraway包中的psid数据集，其中我们关注的是不同性别的组中，教育对收入的影响。

library(party)
data(psid, package = "faraway")
mod1 <- mob(sqrt(income) ~ educ | sex, data = psid)
plot(mod1)

可以观察到男性组别中，教育对收入的影响显然要更大一些。mob函数也可以处理两分类变量，下面的算例来自于mlbench包中的PimaIndiansDiabetes2数据集，其中我们关注的是对糖尿病风险的影响因素。

data(PimaIndiansDiabetes2, package = "mlbench")
data <- na.omit(PimaIndiansDiabetes2[, -c(4, 5)])
mod2 <- mob(diabetes ~ glucose | pregnant + pressure + mass + pedigree + 
    age, data = data, model = glinearModel, family = binomial())
plot(mod2)

从上图可以观察到，首先将数据根据体重变量mass进行了分组。在体重较轻的一组即结点2中，糖尿病发病风险较低，但随着血糖浓度glucose的增加，风险也在急剧增加。在结点4中表示了体重较大但年龄小于30的人群，糖尿病发病风险整体上升。而结点5中表示了体重较大且年龄大于30的人群，糖尿病发病风险达到最大，随着血糖浓度的增加，风险也在逐步增加。
参考资料：
http://cran.r-project.org/web/packages/party/vignettes/MOB.pdf
顺便说下，本文是采用knitr包+markdown+Rstudio生成html后，贴在博客上的，语法高亮木有了，整体效果不算很好。