用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