MIT课程 15.071x The Analytics Edge 第九单元的学习记录。

Integer Optimization

第九单元的主题是整数优化。

1.理论

整数优化

整数优化，即所有解都是整数。
它们有可能是0或者1。这适用于回答是Yes／No的情况。
它们有可能是1，2，3……这适用于回答是具体的数值的情况。

2.实战

做法和线性优化是一样的。只是在条件里面要加一个Integer／Binary的限制。所以就不细讲了。

MIT课程 15.071x The Analytics Edge 第八单元的学习记录。

Linear Optimization

第八单元的主题是线性优化。

1.理论

线性优化

线性优化，其实是用Excel／LibreOffice求解一个简单的多元1次多项式的最大值。
使用LibreOffice是这样做的：

1. 在单元格中填写多项式和约束条件。
2. 选取Tools->Solver，指定多项式，以及各种约束条件，当然也要选择［Linear Solver］这个方法。
3. 点击［Solve］就可以得到结果啦。

在Excel中，需要在option - addin - Excel addin 选择加载solver addin，这样才会在data菜单栏中显示出solver按钮。

sensitivity analysis

sensitivity analysis用来展示，结果是如何随数据（变量／约束条件）的变化而变化的。
shadow prices：表示当需求增加时，将（总量增加量／需求增量）的值定义为影子价格。

如图，纵坐标和横坐标表示两种不同的需求。暗红色阴影表示可能的取值范围。
如果不断提高需求R，从100到125到150，影子价格都保持不变；但是如果需求提高到170，影子价格就会发生变化。
如果不短提高需求D，从150到100，影子价格都为0，总量也不变。

影子价格有可能在需求增加的一个范围内保持不变。也有可能一直为0。

2.实战

当然，用R也能解决这样的问题。

# 先安装pkg
install.packages("lpSolveAPI")
library(lpSolveAPI)

# 创建模型
# 说明：
# 第一个参数是约束条件的个数。
# one capacity constraint: 载客量有一个最大值（飞机座位数）
# two demand constraints : 每种票价的数目（regular seats／discount seats）都有一个最大值
# 所以这个参数的取值是3
# 第二个参数是变量的个数。
# decision variables : 我们有两种票价（regular seats／discount seats）
# 所以这个参数的取值是2

AirlineSimple = make.lp(3,2)
# 创建出来的AirlineSimple是这样子的：
Model name: 
            C1    C2         
Minimize     0     0         
R1           0     0  free  0
R2           0     0  free  0
R3           0     0  free  0
Kind       Std   Std         
Type      Real  Real         
Upper      Inf   Inf         
Lower        0     0         

# 那下面我们就来指定多项式和约束条件
# 最终的效果是这样的：
# max         617*R + 238*D
# subject to    1*R +   1*D <= 166
#               1*R    +   0*D <= 100
#               0*R +   1*D <= 150    

# 特别注意：执行顺序，set.objfn()不能放在前面，我被坑了。。。
# 指定约束条件（跟效果竖着对比着看）
set.column(AirlineSimple, 1, c(1,1,0))
set.column(AirlineSimple, 2, c(1,0,1))
set.constr.type(AirlineSimple, c("<=","<=","<="))
set.rhs(AirlineSimple, c(166,100,150))
# 指定两个变量的参数
set.objfn(AirlineSimple, c(617,238))
# 默认的是最小值，我们改为最大值
lp.control(AirlineSimple,sense='max')

# 这样就创建好了：
Model name: 
            C1    C2         
Maximize   617   238         
R1           1     1  <=  166
R2           1     0  <=  100
R3           0     1  <=  150
Kind       Std   Std         
Type      Real  Real         
Upper      Inf   Inf         
Lower        0     0

# 变量的取值是上面最后两行Upper和Lower，可以通过函数set.bounds()来修改

# 现在可以来运行了
# 如果正确运行，返回值是0
solve(AirlineSimple)
# 查看取得的最大值
get.objective(AirlineSimple)
# 查看取最大值时，两个变量的取值
get.variables(AirlineSimple)

# JFK 从DFW中转，到LAX的场景
# 有8个约束条件，6个变量：
AirlineConnecting = make.lp(8,6)
set.column(AirlineConnecting, 1, c(1,1,1,0,0,0,0,0))
set.column(AirlineConnecting, 2, c(1,1,0,1,0,0,0,0))
set.column(AirlineConnecting, 3, c(1,0,0,0,1,0,0,0))
set.column(AirlineConnecting, 4, c(1,0,0,0,0,1,0,0))
set.column(AirlineConnecting, 5, c(0,1,0,0,0,0,1,0))
set.column(AirlineConnecting, 6, c(0,1,0,0,0,0,0,1))
set.constr.type(AirlineConnecting, rep("<=",8))
set.rhs(AirlineConnecting, c(166,166,80,120,75,100,60,110))
set.objfn(AirlineConnecting, c(428,190,642,224,512,190))
lp.control(AirlineConnecting,sense='max')

# 模型稍微有点大
Model name: 
        C1    C2    C3    C4    C5    C6         
Maximize   428   190   642   224   512   190         
R1           1     1     1     1     0     0  <=  166
R2           1     1     0     0     1     1  <=  166
R3           1     0     0     0     0     0  <=   80
R4           0     1     0     0     0     0  <=  120
R5           0     0     1     0     0     0  <=   75
R6           0     0     0     1     0     0  <=  100
R7           0     0     0     0     1     0  <=   60
R8           0     0     0     0     0     1  <=  110
Kind       Std   Std   Std   Std   Std   Std         
Type      Real  Real  Real  Real  Real  Real         
Upper      Inf   Inf   Inf   Inf   Inf   Inf         
Lower        0     0     0     0     0     0

solve(AirlineConnecting)
get.objective(AirlineConnecting)
get.variables(AirlineConnecting)

今天早上读到两个乌托邦人这篇文章，下面这段很有感触。

1
2
3
4
5
6
7
8
9

赛斯的故事是关于一个人寻找到自己在这世界上的位置。我想起来老舍说，每人在这世间，就像八百尊罗汉，各有各的位置。你实在不能像对桌上的尘土一样，随手拿抹布就把我抹掉了。太平洋也真是座温情的海洋。它允许一个奇葩，这样合理地长大，实现梦想，拥有书、拥有沿海顶楼的好风光、拥有盛满三文鱼的冰箱。这无限孤独又无限美满的人生。

每当我怀疑，我的人生究竟有没有意义的时候，我就会想想赛斯这奇葩。我会诧异如何从一开始每个人想要的都只不过是一个幸福的人生。结果走着走着，就变成了一个“独上高楼、高处不胜寒”的人生，变成了“十年生死两茫茫”的人生，变成了日复一日、年复一年，期待着总有一天能够赢取回报、得偿所愿，然后再终于快乐起来的人生。

我们等下去的时候，他蹦跶着跳进海里捞螃蟹，不知不觉就跑远了。我有些羡慕。

我逐渐发觉心理学以平均值导出关于人性若干结论的好笑之处：人与人之间的差异如此巨大，使得一千个人与一个人的样本，距离真理都同样远。曾有一万人的数据表明勤奋带来成功，再加上三百回研究支持金钱与幸福的钟型曲线关系，可世上奇葩那样多，这些适用于群体的结论就必定不可能与任何一个个体完美匹配。这便是人类群体与小白鼠的区别，而赛斯最早看清。

他让我看到，这样活过一生也OK。

一个人非常立体的生活，在不那么熟悉的眼里都被扁平化了。好比看一场比赛的直播，过程惊心动魄，这是亲身体验；而晚间新闻中，却只有几句湖带过，这是别人眼里的我。更加可怕的是，普通人的一生，投影到历史上，简直微不足道。所以，还是把自己的生活过精彩吧。

孔庙孔府孔林

可能是宗教建筑意外，中国唯一历经各朝战火却能幸存的建筑了。建筑从宋到民国，并且在末尾一个小型的博物馆中看到了乙瑛碑，意外收获。

泰山

传统的红门线路，2个多小时到南天门。虽然不是特别辛苦，但是人挤人真的很无聊，而且并没有多少风景可以看。

泰山石刻

在岱庙见到了向往已久的李斯泰山石刻（仅存9个字），以及明朝的临摹（还有十几个字），原件据说是清朝时从山顶的水池中找到的。真是唏嘘历史的变迁。

日常生活

沉醉于旅途中这样的瞬间，窥见别人的日常生活。许久没有品尝到的人间烟火。
另外一个捷径是逛当地的菜市场，婺源清华镇、绩溪、大同、苏州……明显能感受到勃勃的生机。

1.专注的状态

1
2
3

流体验，指的是当人完全专注地投入一项事业中时，可以从中体验到的极度的快乐。沉浸在这种快乐当中的感受与一般的高兴不同，它是一种难得的高峰体验。人们在这时注意力会完全集中在当下的事物上，感觉不到自我的存在，感觉时间在不知不觉中飞逝，创造力和灵感得到激发，内心非常专注而平静，没有任何冲突，并且从自己所从事的事物中获得极大的满足。所谓的“废寝忘食”，就是沉浸在心流体验中的人们会出现的忘我状态。经常拥有这种高峰体验，会让人的内心得到升华，工作和生活进入到更高的境界，也是人充分实现自身潜能的一个重要途径。

这种高峰体验的出现又一些必要条件：首先就是需要人能毫无杂念地投入，全神贯注；另外，需要所做的事情和个人的能力很好地匹配，若是太过简单，人容易感觉无聊，太过复杂，人则容易感觉受挫；第三，人需要能在这个过程中不断地取得进展，比较明确的目标和及时的反馈都很重要，可以帮助人不断地调整自己的行动。有趣的是，人通常都是在非常积极地为事物努力的时候才会体验到这种极度的快乐。如果只是从事一些被动的活动，比如看电视等，哪怕它们会让人愉快，也不会激发心流。除此之外，精力涣散和心不在焉也会减少人们从事某项事物中获得的乐趣。

2.进入状态／保持状态

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

程序员们都知道，任务切换需要耗费许多额外的花销，通俗地来讲，首先需要保存当前上下文以便下次能够顺利切换回来，然后要加载目标任务的上下文。如果一个系统不停地在多个任务之间来回倒腾，就会耗费大量的时间在上下文切换上，无形中浪费很多的时间。

相比之下，如果只做一件任务，就不会有此损失。这就是为什么专注的人比不专注的人时间利用效率高得多的原因。任务切换的暗时间看似非常不明显，甚至很多人认为“多任务”是件很好的事情（有时候的确是），但日积月累起来就会发现，消耗在切换上的时间越来越多。

另外，大脑开始一件任务的时候必须要有一定时间来“热身”，这个时间因人而异，并且可以通过练习来改变。举个例子，你看了一会书之后，忽然感到一阵无聊，忍不住打开浏览器，十分钟后你想起来还要继续看书，但要回复到当时理想的状态，却需要一段时间来努力去集中精力，把记忆中相关的知识全都激活起来，从而才能进入“状态”，因为你上了十分钟网之后这些记忆已经被抑制了。如果这个“热身”状态需要一刻钟，那么看似十分钟的上网闲逛其实就花费了二十五分钟。

如果阅读的例子还不够生动，对于程序员来说其实有更好的例子：你写程序写得正high，忽然被叫去开了一通会，写到一半的代码搁在那儿。等你开完会回来你需要多久能够重新进入状态？又或者，你正在调试程序，你已经花了二十分钟的时间把与这个bug可能相关的代码前前后后都理解了一遍，心中构建了一个大致的地图，就在这时，呃，你又被叫去开了个会(:D)，开完会回来，可想而知，得花上一些时间来回想一下刚刚弄清的东西了。

迅速进入状态的能力是可以锻炼的，根据我个人的经验，至少可以缩短到3-5分钟。但要想完全进入状态，却是很难在这么短的时间实现的。所谓完全进入状态，举个例子：你看了3个小时的书，或者调试了半个小时的程序之后，往往满脑子都是相关的东西，所有这些知识都处在活跃状态，换言之你大脑中所有相关的记忆神经网络都被激活了，要达到这样一种忘记时间流逝的“沉浸”状态（心理学上叫做“流体验”），不是三两分钟的事情。而一旦这种状态被破坏，无形间效率就会大打折扣。这也是为什么我总是倾向于创造大块的时间来阅读重要的东西，因为这样有利于“沉浸”进去，使得新知识可以和大脑中与其相关的各种既有的知识充分融合，关联起来，后者对于深刻的记忆非常有帮助。

要充分利用暗时间，不仅要能够迅速进入状态，另一个很重要的习惯就是能够保持状态多久（思维体力）。

能够迅速进入专注状态，以及能够长期保持专注状态，是高效学习的两个最重要习惯。
值得庆幸的是，专注力和耐力与才能不同，可以通过训练于后天获得可以不断提升其资质。只要每天坐在书桌前，训练将意识倾注于一点，自然就能掌握。这同前面写过的强化肌肉的做法十分相似。每天不间断地写作集中意识去工作，这些非做不可----将这样的信息持续不断地传递给身体系统，让它牢牢地记住，再悄悄移动刻度，一点一点将极限值将上提升，注意不让身体发觉。这跟每天坚持慢跑，强化肌肉，逐步打造出跑步者的体型，乃是异曲同工。给它刺激，持续。再给它刺激，持续。这一过程当然需要耐心，不过一定会得到相应的回报。

优秀的侦探小说家雷蒙德.钱德勒在私信中说过：“哪怕没有什么东西可写，我私吞也肯定在书桌前坐上好几个小时，独自一人集中精力。”他这么做是为了什么，我完全能理解。钱德勒通过这么做，来提高职业作家必需的膂力，静静地提高士气。这样一种日常训练对他必不可缺。

我认为写作长篇小说是一种体力劳动。写文章属于脑力劳动，然而写出一本大部头来，更近于体力劳动。诚然，写书并不需要举起沉重的物体，也不需要飞速地奔来跑去，高高地蹿上跳下。世间的很多人似乎只看到表面，将作家的工作视为宁静而更改的书斋劳动，以为有了足以端起一只咖啡杯的力量，就能写小说了。试它一试，立即就会明白，写小说并非那么安逸的工作。坐在书桌前，将神经如同激光束一般集于一点，动用想象力，从“无”的地平线上催生出故事来，挑选出一个个正确的词语，让所有的流程准确无误----这样一种工作，与一般人想象的相比，更为长久地需要远为巨大的能量。这固然不必运动身体，劳筋动骨的劳动却在体内热火朝天地展开。固然，思索问题的是脑子，小说家却需披挂着叫“故事”的全副装备，动用全身进行思考，这要求作家彻底地驱使----在许多时候是奴役----肢体能力。

3.排除干扰

1

干吗还会欣羡我能坐稳在咖啡馆里呢？到咖啡馆原本就为着隔离而来，隔离自己的家，隔离善良的声音，隔离掉所有熟悉、舒适、温暖的东西；正在写长篇的小说家林俊颖一人独居，如今却也冲出到咖啡馆来，他笑着说，书架一直在那里叫你，你一碰到困难，借口翻翻资料，寻找感觉，接下来你就发现自己又埋进某本书、某部小说里两小时了。所以，所有像回事的作家最终几乎都在早上书写，趁着整个世界才刚醒来，还跟你暂时处在一种相互隔离的状态，你还有能力把它当在外头——就连海明威这种浮夸好热闹的人都告诉我们，在早晨进入写作之前，不做其他任何有企图心的事；纳博科夫一致工作到下午，知道黄昏散步时才找报纸看，才放世界溜进来；在淡水写作的舞鹤甚至不读报，他只在喂食镇上街猫时顺便瞄一眼头条，知道没发生战争，末日还没来就可以了。

4.养成习惯

1

书中主人翁流落赌城，在绝望时刻偶然从一个老头手上得到一个必然赢钱的赌方，但这个最后一定大赢的赌方非常诡异非常磨人，它必须先挨过一定阶段的输钱，只能输不能赢，而且明知是输亦一步也不能省——写小说的格林迷朱天心尤其喜欢这个例子，她在新小说顺利开笔之前，一样总要经历这同样的短则数日长可数星期的枯坐思索（在小说题材乃至内容已经完全锁定备妥的情况下），明明知道一定空手而回仍得每天带着书、草稿本和笔到写作的咖啡馆报到，她出门时的口头禅便是：“去输钱”。

善哉行

曹丕

上山采薇，薄暮苦饥。

溪谷多风，霜露沾衣。

野雉群雊，猿猴相追。

还望故乡，郁何垒垒！

高山有崖，林木有枝。

忧来无方，人莫之知。

人生如寄，多忧何为？

今我不乐，岁月如驰。

汤汤川流，中有行舟。

随波转薄，有似客游。

策我良马，被我轻裘。

载驰载驱，聊以忘忧。

MIT课程 15.071x The Analytics Edge 第七单元的学习记录。

Visualization

第七单元的主题是可视化。

1.简介

plot和ggplot2的比较
plot：只有简单的点和线，不容易添加其他元素。
ggplot2：引入图层，很容易添加其他元素

ggplot2

ggplot2三要素：

1. Data
   数据，使用data.frame。
2. Aesthetic mapping
   指定如何将 data.frame里的变量映射到图形属性上。比如，颜色，形状，比例，x／y坐标，分组等等。
3. Geometric objects
   决定数据以什么样的形式显示。比如，点，线，箱线图，条形图，多边形等等。

结合下面这条命令，参数WHO就是提供数据的data.frame，参数aes()就是Aesthetic mapping，后面用加号连结的类似geom_point()就是Geometric objects。

# 形式
# ggplot(data = NULL, mapping = aes(), ..., environment = parent.frame())
# 例子
ggplot(WHO, aes(x = GNI, y = FertilityRate, color = Region)) + geom_point()

ase()即可以作为ggplot()的参数，又可以作为geom_XXXX()的参数

Aesthetic mapping

坐标相关

aes(x, y, xmin, xmax, ymin, ymax, xend, yend)
# 当然就是x，y坐标分别指定data.frame的某一列

注：坐标相关的，一般作为ggplot()的参数，其他的都可以作为geom()的参数。

Geometric objects

颜色相关

aes(colour, fill, alpha)
# colour 颜色
# fill   填充指标，data.frame的某一列。也类似于分类，比如该列有两个因子，那么会用两种不同的颜色填充
# alpha  透明度，0到1之间的小数

分组相关

aes(group)
# group 分组指标，可以指定为1，那所有数据都在1组。也可以指定data.frame的某一列

形态相关

aes(linetype, size, shape)
# linetype 即lty，线段的类型
# size     点的大小，线的粗细。指定整数数值。
# shape    图形的类型

图形的类型，即geom_point(shape = n)中n的取值

线段的类型，即geom_point(lty = n)中n的取值

描绘形状

geom_point()  点
geom_line()   线
geom_tile()   条形图
geom_bar()    直方图
geom_ploygen()多边形

注：
binwidth = 5 :粒度？
geom_bar(stat=”identity”) :use the value of the y variable as is
geom_histogram(position = “identity”) :not to stack the histograms

2.实战

绘图

# Read in data
WHO = read.csv("WHO.csv")
str(WHO)

# Plot from Week 1
plot(WHO$GNI, WHO$FertilityRate)

# Let's redo this using ggplot 
# Install and load the ggplot2 library:
install.packages("ggplot2")
library(ggplot2)

# Create the ggplot object with the data and the aesthetic mapping:
scatterplot = ggplot(WHO, aes(x = GNI, y = FertilityRate))

# Add the geom_point geometry
scatterplot + geom_point()

# Make a line graph instead:
scatterplot + geom_line()

# Switch back to our points:
scatterplot + geom_point()

# Redo the plot with blue triangles instead of circles:
scatterplot + geom_point(color = "blue", size = 3, shape = 17) 

# Another option:
scatterplot + geom_point(color = "darkred", size = 3, shape = 8) 

# Add a title to the plot:
scatterplot + geom_point(colour = "blue", size = 3, shape = 17) + ggtitle("Fertility Rate vs. Gross National Income")

分组

# 因子，以颜色区分    
# Color the points by region: 
ggplot(WHO, aes(x = GNI, y = FertilityRate, color = Region)) + geom_point()

# 数值，以颜色深浅区分
# Color the points according to life expectancy:
ggplot(WHO, aes(x = GNI, y = FertilityRate, color = LifeExpectancy)) + geom_point()

拟合

# Is the fertility rate of a country was a good predictor of the percentage of the population under 15?
ggplot(WHO, aes(x = FertilityRate, y = Under15)) + geom_point()

# Let's try a log transformation:
ggplot(WHO, aes(x = log(FertilityRate), y = Under15)) + geom_point()

# Simple linear regression model to predict the percentage of the population under 15, using the log of the fertility rate:
mod = lm(Under15 ~ log(FertilityRate), data = WHO)
summary(mod)

# Add this regression line to our plot:
ggplot(WHO, aes(x = log(FertilityRate), y = Under15)) + geom_point() +     stat_smooth(method = "lm")

# 99% confidence interval
ggplot(WHO, aes(x = log(FertilityRate), y = Under15)) + geom_point() + stat_smooth(method = "lm", level = 0.99)

# No confidence interval in the plot
ggplot(WHO, aes(x = log(FertilityRate), y = Under15)) + geom_point() + stat_smooth(method = "lm", se = FALSE)

# Change the color of the regression line:
ggplot(WHO, aes(x = log(FertilityRate), y = Under15)) + geom_point() + stat_smooth(method = "lm", colour = "orange")

热力图

热力图（数据越多颜色越深）的效果，依靠scale_fill_gradient()来实现，可以通过low和high指定深浅区域的颜色，然后自动形成渐变效果。旁边的图例通过参数guide = “legend”来指定。
最终的命令如下，如何生成数据的，就不啰嗦了。

# Change the color scheme
ggplot(DayHourCounts, aes(x = Hour, y = Var1)) + geom_tile(aes(fill = Freq)) + scale_fill_gradient(name="Total MV Thefts", low="white", high="red") + theme(axis.title.y = element_blank())

地理热力图

顾名思义，地理热力图就是在地图上显示热力图。
包map内置了美国地图、世界地图、法国地图、意大利地图等。地图的原理跟图片类似，图片就是按照某个粒度分成很多个像素点，然后保存像素点的颜色信息；地图就是按照经纬度分成很多点，保存每个点的信息（比如这个点位于哪个州，这样就形成一个美国地图）。
对比刚才的 ggplot() + geom_tile() + scale_fill_gradient()
我们现在使用 ggmap() + geom_point() + scale_fill_gradient()

# Install and load two new packages:
install.packages("maps")
install.packages("ggmap")
library(maps)
library(ggmap)

# Load a map of Chicago into R:
chicago = get_map(location = "chicago", zoom = 11)

# Look at the map
ggmap(chicago)

# Plot the first 100 motor vehicle thefts:
ggmap(chicago) + geom_point(data = mvt[1:100,], aes(x = Longitude, y = Latitude))

# Round our latitude and longitude to 2 digits of accuracy, and create a crime counts data frame for each area:
LatLonCounts = as.data.frame(table(round(mvt$Longitude,2), round(mvt$Latitude,2)))

str(LatLonCounts)

# Convert our Longitude and Latitude variable to numbers:
LatLonCounts$Long = as.numeric(as.character(LatLonCounts$Var1))
LatLonCounts$Lat = as.numeric(as.character(LatLonCounts$Var2))

# Plot these points on our map:
ggmap(chicago) + geom_point(data = LatLonCounts, aes(x = Long, y = Lat, color = Freq, size=Freq))

# Change the color scheme:
ggmap(chicago) + geom_point(data = LatLonCounts, aes(x = Long, y = Lat, color = Freq, size=Freq)) + scale_colour_gradient(low="yellow", high="red")

# We can also use the geom_tile geometry
ggmap(chicago) + geom_tile(data = LatLonCounts, aes(x = Long, y = Lat, alpha = Freq), fill="red")

云图

# 先准备下数据，我们需要很多单词。
# 跟文本处理类似，依旧使用tweets推文，只是我们这次不抽取词干。
library(tm)
tweets = read.csv("tweets.csv", stringsAsFactors=FALSE)
corpus = Corpus(VectorSource(tweets$Tweet))
corpus = tm_map(corpus, tolower)
corpus = tm_map(corpus, PlainTextDocument)
corpus = tm_map(corpus, removePunctuation)
corpus = tm_map(corpus, removeWords, stopwords("english"))
frequencies = DocumentTermMatrix(corpus)
allTweets = as.data.frame(as.matrix(frequencies))

# 我们需要的单词就是列名
colnames(allTweets)
# 我们需要的另一个指标是单词的频率
colSums(allTweets)

# 现在加载wordcloud这个包
library(wordcloud)
wordcloud(colnames(allTweets), colSums(allTweets), scale=c(2, .25))

# 参数 scale 指定了文字的大小
# scale=c(2, .25) 表示出现频率最高的单词，显示的字号为2，出现频率最小的单词，显示的字号为0.25
wordcloud(colnames(allTweets), colSums(allTweets))
# 等效于
wordcloud(colnames(allTweets), colSums(allTweets), scale=c(4, 0.5))

# min.freq
# 只显示出现频率大于指定值的单词

# max.words
# 最多只显示指定数目的单词

# random.order == FALSE
# 最先显示出现频率最高的单词

# rot.per = 0.5
# 有一半的单词垂直显示。默认值是0.1。

# random.color == TRUE
# 使用随机颜色

颜色
包RColorBrewer支持下面这些调色板，可以输入 display.brewer.all() 看到下面这张图。

ibrary(RColorBrewer)
display.brewer.all()

# 像这样使用
colors=brewer.pal(9, "Blues")[5:9]
wordcloud(colnames(allTweets), colSums(allTweets), colors)

保存

# Save our plot:
fertilityGNIplot = scatterplot + geom_point(colour = "blue", size = 3, shape = 17) + ggtitle("Fertility Rate vs. Gross National Income")
pdf("MyPlot.pdf")
print(fertilityGNIplot)
dev.off()

附录

R中星期的显示

在中文系统上，weekdays()返回的结果是 “星期二 星期六 星期日 星期三 星期四 星期五 星期一”，如果希望输出的结果是“Friday Monday Saturday Sunday Thursday Tuesday Wednesday”，应该怎么做？

# Convert the Date variable to a format that R will recognize:
mvt$Date = strptime(mvt$Date, format="%m/%d/%y %H:%M")
mvt$Weekday = weekdays(mvt$Date)

table(mvt$Weekday)
星期二 星期六 星期日 星期三 星期四 星期五 星期一 
26791  27118  26316  27416  27319  29284  27397 

Sys.getlocale()
"zh_CN.UTF-8/zh_CN.UTF-8/zh_CN.UTF-8/C/zh_CN.UTF-8/zh_CN.UTF-8"
Sys.setlocale("LC_TIME", "en_US.UTF-8")
"en_US"
Sys.getlocale()
"zh_CN.UTF-8/zh_CN.UTF-8/zh_CN.UTF-8/C/en_US.UTF-8/zh_CN.UTF-8"

table(mvt$Weekday)
Friday    Monday  Saturday    Sunday  Thursday   Tuesday Wednesday 
29284     27397     27118     26316     27319     26791     27416

另外注意到，不管是中文还是英文，都是按照字母表顺序排列的，不是按照实际中有意义的顺序排列的。

WeekdayCounts = as.data.frame(table(mvt$Weekday))
WeekdayCounts$Var1 = factor(WeekdayCounts$Var1, ordered=TRUE, levels=c("Sunday", "Monday", "Tuesday", "Wednesday", "Thursday", "Friday","Saturday"))

factor转数字

先把factor转成character，再转成数字

# Convert the second variable, Var2, to numbers and call it Hour:
DayHourCounts$Hour = as.numeric(as.character(DayHourCounts$Var2))

参考：
形状和线段的类型
颜色)

MIT课程 15.071x The Analytics Edge 第六单元的学习记录。

Clustering

第六单元的主题是集群。它用来找到数据内的相似性。

1.理论

Recommendation Systems

Collaborative filtering:
过滤出用户间的共同特征／相似性。只使用了用户信息，跟电影内容本身无关。

Content filtering:
利用电影本身的信息，过滤出有共同导演／演员／类别的电影。跟其他用户无关。

clustering

clustering 集群是一种非监督学习，”unsupervised learning”，将有共同特征的数据分在同一组。

Hierarchical clustering

Hierarchical clustering的步骤：

1. 计算距离
2. 生成集群
3. 生成cutree

注意1:计算距离时，有可能造成内存溢出。计算每两点间的距离，得到的结果是n*(n-1)/2个，我们需要保存这个结果，如果n很大，保存结果的矩阵也很大，可能会导致内存溢出。
注意2:计算距离的三种方法：
Euclidean distance：点与点之间的欧几里得距离
Manhattan Distance：绝对值之和
Maximum Coordinate：偏离最严重的点

K-means clustering

K-means clustering的步骤：

1. 指定集群数目k
2. 随机分配所有的点
3. 计算每个集群的中心点
4. 计算每个点到这些中心点的距离，选择最近的，重新分配点到离他最近的集群
5. 重新计算每个集群的中心点
6. 重复4和5多次，直到没有提升

注意：centroid distance 集群中所有点的平均值间的距离。

normalize

如果不同列的数值不是同样的数量级，那么运算后较小的值可能会被忽略，所以需要调整到同样的数量级。

library(caret)
preproc = preProcess(airlines)
airlinesNorm = predict(preproc, airlines)

效果就是，所有列的平均值都是0。

2.建模和评估

Hierarchical clustering

# After following the steps in the video, load the data into R
movies = read.table("movieLens.txt", header=FALSE, sep="|",quote="\"")
# Add column names
colnames(movies) = c("ID", "Title", "ReleaseDate", "VideoReleaseDate", "IMDB", "Unknown", "Action", "Adventure", "Animation", "Childrens", "Comedy", "Crime", "Documentary", "Drama", "Fantasy", "FilmNoir", "Horror", "Musical", "Mystery", "Romance", "SciFi", "Thriller", "War", "Western")
# Remove unnecessary variables
movies$ID = NULL
movies$ReleaseDate = NULL
movies$VideoReleaseDate = NULL
movies$IMDB = NULL
# Remove duplicates
movies = unique(movies)

# Compute distances
distances = dist(movies[2:20], method = "euclidean")

# Hierarchical clustering
# clusterMovies = hclust(distances, method = "ward") 
clusterMovies = hclust(distances, method = "ward.D")

# Plot the dendrogram
plot(clusterMovies)

# Assign points to clusters
clusterGroups = cutree(clusterMovies, k = 10)
# Create a new data set with just the movies from cluster 2
cluster2 = subset(movies, clusterGroups==2)

K-means clustering

healthy = read.csv("healthy.csv", header=FALSE)
# 注意
# data.frame->matrix->vector 变成一个2500的vector
# data.frame->vector 还是一个50*50的data.frame
healthyMatrix = as.matrix(healthy)
healthyVector = as.vector(healthyMatrix)

# Specify number of clusters
k = 5
# Run k-means
set.seed(1)
KMC = kmeans(healthyVector, centers = k, iter.max = 1000)

# Extract clusters
healthyClusters = KMC$cluster

# Plot the image with the clusters
dim(healthyClusters) = c(nrow(healthyMatrix), ncol(healthyMatrix))

image(healthyClusters, axes = FALSE, col=rainbow(k))

# Apply to a test image
tumor = read.csv("tumor.csv", header=FALSE)
tumorMatrix = as.matrix(tumor)
tumorVector = as.vector(tumorMatrix)

# Apply clusters from before to new image, using the flexclust package
# kcca K-Centroids Cluster Analysis
install.packages("flexclust")
library(flexclust)
KMC.kcca = as.kcca(KMC, healthyVector)
tumorClusters = predict(KMC.kcca, newdata = tumorVector)

# Visualize the clusters
dim(tumorClusters) = c(nrow(tumorMatrix), ncol(tumorMatrix))
image(tumorClusters, axes = FALSE, col=rainbow(k))

MIT课程 15.071x The Analytics Edge 第五单元的学习记录。

Text Analytics

第五单元的主题是文本分析。

1.理论

Bag of Words

一段文本，可以看作是多个单词的集合。
统计这些单词的特征，可以归纳文本的倾向。

首先，我们需要对文本进行下面这几步预处理：

1. clean up irregularities(统一大小写)
2. remove punctuations(去掉标点或者特殊符号)
3. remove stop words(去掉the／who／is／do这些单词)
4. stemming(获取词干，也就是去除动词变形，比如agrued，agrues，agruing，都变成agru)

然后，我们统计文本中剩下这些单词的出现次数，生成一个矩阵，类似这样的格式：

在实际中，生成的矩阵是个稀疏矩阵（有很多0），我们只选取出现次数比较多的，忽略那些不常见的单词。
比如选取至少出现过20次的单词，其他的忽略。
这样的矩阵，每列的列名就是自变量，矩阵的值就用做自变量的取值。

最后，手动添加一列，作为因变量，这样就可以根据这些单词的出现次数，预测因变量的取值了。
所以，这一列因变量的数值如何定义，它的实际意义是什么，其实是比较复杂的。

在课程的例子中，它定义了”好感度”，并且只有下面五种取值，{-2,-1,0,1,2}，最终要建立模型预测哪些文本暗示发推的人对苹果公司很没有好感（好感度是－2）。最终发现，文本中含有”hate”,”wtf”的情况，推主对苹果公司很没有好感。😄

IBM Watson

Watson的工作步骤是这样的：

1. Find LAT
   首先得搞明白问题是什么，也就是要找到问题的LAT(Lexial Answer Type)。
   问题”Mozart’s last and perhaps most powerful symphony shares its name with this planet.”的LAT是”this planet”，因为把答案”Jupiter”替换进原来的句子，
   “Mozart’s last and perhaps most powerful symphony shares its name with Jupiter”
   仍然是说得通的。
2. Generate Hypothesis
   在数据库中搜索上百个候选答案，替换掉LAT，生成很多假说。
3. Score Hypothesis
   对每个假说，进行文本搜索，可以将搜索的到的结果数目作为评分。
4. Rank Hypothesis
   对评分进行排序，选取评分最高的那个作为答案。

2.建模和评估

预处理

# Read in the data
# 不要把文本转化为因子
tweets = read.csv("tweets.csv", stringsAsFactors=FALSE)

# Create dependent variable
tweets$Negative = as.factor(tweets$Avg <= -1)

# Install new packages
install.packages("tm")
library(tm)
install.packages("SnowballC")
library(SnowballC)

# Create corpus
corpus = Corpus(VectorSource(tweets$Tweet))

# Convert to lower-case
corpus = tm_map(corpus, tolower)
corpus = tm_map(corpus, PlainTextDocument)

# Remove punctuation
corpus = tm_map(corpus, removePunctuation)

# Remove stopwords and apple
corpus = tm_map(corpus, removeWords, c("apple", stopwords("english")))

# Stem document 
corpus = tm_map(corpus, stemDocument)

统计，生成单词出现次数的矩阵

# Create matrix
frequencies = DocumentTermMatrix(corpus)

# Look at matrix 
inspect(frequencies[1000:1005,505:515])

# Check for sparsity
# 找出出现次数至少有20次的单词
findFreqTerms(frequencies, lowfreq=20)

# 忽略99.5%的稀疏数据，只选取0.5%作为有效数据
# Remove sparse terms 
sparse = removeSparseTerms(frequencies, 0.995)

# Convert to a data frame
tweetsSparse = as.data.frame(as.matrix(sparse))
# Make all variable names R-friendly
colnames(tweetsSparse) = make.names(colnames(tweetsSparse))

# Add dependent variable
tweetsSparse$Negative = tweets$Negative

建模和评估

# Split the data
library(caTools)
set.seed(123)
split = sample.split(tweetsSparse$Negative, SplitRatio = 0.7)
trainSparse = subset(tweetsSparse, split==TRUE)
testSparse = subset(tweetsSparse, split==FALSE)

# Build a CART model
library(rpart)
library(rpart.plot)
tweetCART = rpart(Negative ~ ., data=trainSparse, method="class")
# Evaluate the performance of the model
predictCART = predict(tweetCART, newdata=testSparse, type="class")
table(testSparse$Negative, predictCART)


# Random forest model
library(randomForest)
set.seed(123)
tweetRF = randomForest(Negative ~ ., data=trainSparse)
# Make predictions:
predictRF = predict(tweetRF, newdata=testSparse)
table(testSparse$Negative, predictRF)

后人看历史的时候，不可避免地站在上帝视角思考问题。
读史记，一个人将来能够成就什么，从一开始就预设好了。
又或者大部分历史书，一个王朝如何兴起、成长、强盛、衰落，仿佛就是按照剧本来演的。

而实际上，当时的人，哪知道后来会发生什么，哪知道做了某件事情的结果，甚至不知道大国的另外一个角落，此刻发生了什么。

就好像看比赛直播，心情各种复杂，岂止是晚间新闻报道里三言两语能够概括的。

亚马逊链接：眼前：漫游在《左传》的世界