Beta

有点标题党的嫌疑，实际是介绍如何使用 R 绘制 heatmap 的文章。

今天无意间在Flowingdata看到一篇关于如何使用 R 来做 heatmap 的文章（请移步到这里）。虽然 heatmap 只是 R 中一个很普通的图形函数，但这个例子使用了2008-2009赛季 NBA 50个顶级球员数据做了一个极佳的演示，效果非常不错。对 R 大致了解的童鞋可以直接在 R console 上敲

？heatmap

直接查看帮助即可。

没有接触过 R 的童鞋继续围观，下面会仔细介绍如何使用 R 实现 NBA 50位顶级球员指标表现热图：

A heat map is a graphical representation of data where the values taken by a variable in a two-dimensional map are represented as colors.Heat maps originated in 2D displays of the values in a data matrix. Larger values were represented by small dark gray or black squares (pixels) and smaller values by lighter squares.

先解释一下数据：

这里共列举了50位球员，估计爱好篮球的童鞋对上图右边的每个名字都会耳熟能详。这些球员每个人会有19个指标，包括打了几场球（G)、上场几分钟（MIN)、得分（PTS)……这样就行成了一个50行×19列的矩阵。但问题是，数据有些多，需要使用一种比较好的办法来展示，So it comes, heatmap!

简单的说明：

比如从上面的热图上观察得分前3名（Wade、James、Bryant）PTS、FGM、FGA比较高，但Bryant的FTM、FTA和前两者就差一些；Wade在这三人中STL是佼佼者；而James的DRB和TRB又比其他两人好一些……

姚明的3PP（3 Points Percentage）这条数据很有意思，非常出色！仔细查了一下这个数值，居然是100%。仔细回想一下，似乎那个赛季姚明好像投过一个3分，并且中了，然后再也没有3p。这样本可真够小的！

最后是如何做这个热图（做了些许修改）：

Step 0. Download R

R 官网：http://www.r-project.org，它是免费的。官网上面提供了Windows,Mac,Linux版本（或源代码）的R程序。

Step 1. Load the data

R 可以支持网络路径，使用读取csv文件的函数read.csv。

read.csv("http://datasets.flowingdata.com/ppg2008.csv", sep=",")

Step 2. Sort data

按照球员得分，将球员从小到大排序：

nba <- nba[order(nba$PTS),]

当然也可以选择MIN,BLK,STL之类指标

Step 3. Prepare data

把行号换成行名（球员名称）：

row.names(nba) <- nba$Name

去掉第一列行号：

nba <- nba[,2:20] # or nba <- nba[,-1]

Step 4. Prepare data, again

把 data frame 转化为我们需要的矩阵格式：

nba_matrix <- data.matrix(nba)

Step 5. Make a heatmap

R 的默认还会在图的左边和上边绘制 dendrogram，使用Rowv=NA, Colv=NA去掉

heatmap(nba_matrix, Rowv=NA, Colv=NA, col=cm.colors(256), revC=FALSE, scale='column')

这样就得到了上面的那张热图。

Step 6. Color selection

或者想把热图中的颜色换一下：

heatmap(nba_matrix, Rowv=NA, Colv=NA, col=heat.colors(256), revC=FALSE, scale="column", margins=c(5,10))

延伸阅读：

来自于kerimcan和krees这些人的讨论：

http://sekhon.polisci.berkeley.edu/stats/html/heatmap.html http://enotacoes.wordpress.com/2007/11/16/easy-guide-to-drawing-heat-maps-to-pdf-with-r-with-color-key/

补充：

早上起来发现 David Smith 同样更新了博客。唉，这厮嗅觉也忒灵敏！哈哈

本场比赛前7分钟火箭发挥还不错，最高取得了10分的领先（11-21），但受上一场力拼森林狼三个加时影响，火箭诸将体能逐渐不支，慢慢失去优势。虽然巴丁格整场替补发挥出色，无奈，随着阿里扎上篮不进，比分定格在了115-106。

相比热火发烫的53.6%投篮命中率，火箭发挥比较正常，48.7%。火箭唯一问题出现在了失误方面，8-15，如果火箭失误控制的好，也许结果可能是另外一个结果。

比赛过程中，杨毅提到（大致意思）：由于没有超级球星，火箭必须比其他球队付出更多的努力才能获得胜利，也就是说其他球队会从容地为季后赛调整状态。火箭没有这个资本，为了常规赛的成绩，火箭进入疲劳期的时间要更早。要保证后面的球队成绩，要么莫雷尽最大能力交易，补充火箭；要么阿德尔曼通过更加合理细致的轮转，让每个核心球员得到充分休息。

问题来了：在不改变现有火箭球员结构的前提下，火箭的战术轮转体系中，球员的位置如何？

引子：

本场比赛火箭一共20个助攻，Brooks 和 Battier 分别助攻了最高的5个和4个。每次助攻都会涉及两位球员，那么本场比赛所有助攻结果综合在一起，即我们将助攻者和被助攻者之间的关系使用社会网络关系表现出来，会有一些有趣的现象：

注释：

箭头方向是助攻方向，比如最下面的是 Shane Battier 给 Chuck Hayes 的助攻。

1. Brooks 无疑是比赛的发起者，我们发现他的助攻几乎包括了中锋、前锋位置的所有人（但不包括阿里扎，好像我记得有个镜头阿里扎要球，Brooks 没有理会）。
2. Battier 既是助攻的受益者，又是助攻的发起者。个人一直比较喜欢的球员，篮球智商非常高。
3. Budinger 这场比赛发挥出色，同队友给予其的帮助分不开。我们看到很多个球员对其都有直接帮助。
4. Andersen 从助攻网络关系角度看，属于一个策应型中锋，而且是由里向外策应的那类。从比赛中观察，似乎球风有些偏软（本赛季我第一次看直播比赛 ^_^）
5. Ariza 接受的助攻并不多，只有 Battier 的一次，其他都是给别人的助攻，和 Brooks 一样，属于个人能力比较强，擅于自己创造得分机会的球员。

假如：

1. Brooks、Battier、Ariza、Andersen 在火箭进攻体系中位置比较重要，轮转的时候尽量保证其中的两人或三人同时在场。
2. Budinger 属于绝好的替补球员，但似乎不适合同 Battier、Hayes 同时在场。
3. Lowry 在组织进攻方面能力欠佳，使用上须谨慎。

更新分割:

40场比赛助攻数据同时考虑，结果有些凌乱：

注：这个赛季火箭队有一些球员实际上并没有真正的进入轮转，比如"Tracy McGrady","Mike Harris","Jermaine Taylor","Pops Mensah-Bonsu","Brian Cook"。虽然有些球员（比如麦蒂）的确对球队的（被）助攻仍有帮助，但贡献非常有限。出于结果整洁性的考虑，上图已将这些球员因素剔除。

由于绘图算法使用的是 Force-based_algorithms，也就是说这种算法做出的图，边（edges）会尽可能的少。解释为，对球员关系的影响就是：

同其他球员关系比较多的球员将绘制的比较靠近中心，而关系较少的球员会绘制在相对靠外的位置。

重新观察火箭队助攻网络图，发现：

1月16日对热火比赛中，Kyle Lowry 和 Trevor Ariza 发挥的确出了问题，尤其是 Lowry 这点上。

Brooks、Ariza、Scola、Lowry、Battier 在助攻重要性角度上，占据球队的前五位。如何使用“田忌赛马”的策略制胜，则是教练组的问题。

让我欣慰的是 Hayes 的重要程度要比 Andersen 要好，单单从上一场比赛上看，Andersen 发挥的有些超常。

P.S. 维基百科上没有区别 “有向网络”和“无向网络”的 page rank ，上个表中的 page rank 值属于“无向网络”值，同上面的图略有区别（有向网络中，Lowry 的关键性仅比 Hayes 高，有些无奈）。

以前就想写一篇博客，讲述 Google 给我们生活带来的便利，这不，再不写，也许再也没有机会了。为什么说没机会了，是因为 Google 官方博客的一篇文章——A new approach to China（抱歉，因为是订阅的内容，已经被墙了，我实在不能找到链接，不过可以访问这里，中文翻译的和原文），也许真的有一天 Google 会撇下我们。

说到 Google，感慨颇深。以前在人大读书的时候，受舍友影响，从来不知道有其他的搜索引擎，不论干啥，第一反应就是 Google 之。后来，Google 进入中国，正式提出”谷歌“的中文名称，当时我还笑话 Google 的中文名太傻。而现在呢，患有严重的 Google 依赖综合征，算是交代了。

Google 最适宜对英文资料的搜索，很多童鞋对 google.com 和 google.cn 感觉一样，实际上是有区别的。用一句很低俗很低俗的话来说就是，google.cn 是被阉割过的 Google（虽然仍然比 baidu 好很多）。举个最简单的例子，google.cn 是没有账户信息的，也就是说，你不能通过 google.cn 来登录 Google 的服务。而且似乎在中国，在浏览器中使用 Google 会默认指向 Google 中国。再换句话说，我们一般会把 Google 认为是一个很普通的搜索引擎，但实际上如果登录 Google，会发现别有洞天。

• Gmail：Google 的 一款优秀的 mail。还记得最早 Google 开放邮件系统，没有独立注册的地方，必须通过其他人邀请。
• Google Reader：每天开电脑后，要做的第二件事（第一件事是 foxmail 收取 gmail）。信息需要捕捉，使用书签记录互联网信息那是 web1.0 时代，现在我们有 RSS，Google Reader 可以带领我们翻墙去学统计，sigh！不过有个问题就是，如果你恰好在某个站点看到了一个很好的文章，而上面又有一个pdf链接……点击，浏览器报告错误链接……噢，那是在墙外。哪位童鞋有好办法解决，请告之。
• Google Group：顾名思义，是小组讨论的论坛，很多志同道合的童鞋发言交流思想的地方。比如，申请了 TopLanguage，COS R Team 等，不过我更倾向把它归为 maillist 的一种，就和 R 的Mailing Lists 一样。
• Google Alert：Google 中国翻译成“快讯”，恰到好处。它能告诉你每天最新发生的事情，当然发生的事情是你用关键词来定义的。

办公类：

• Google Notebook：真正意义的互联网笔记本。
• Google Docs：美国华盛顿特区政府官方使用的办公软件。挺好的，国内不知为啥又被封掉。
• Google Talk：非常适合办公环境的即时聊天软件，聊天记录保存在 Gmail 帐号里，支持语音功能，同时有 gmail 邮件通知。
• Google code：直接参考 R 的 sqldf 包 http://code.google.com/p/sqldf/

网站类：

• Google Analytics：做网站的朋友肯定对它非常熟悉。
• Google calendar：日程，合理规划时间是成功的必要条件。
• Google Site 和 Google pages：功能上感觉比较类似，都是用来做站点的，但都不能用了，sigh again!
• Goolge Picasa：图片分享，也不能用了
• Google blog：其实挺好的 blog，由于众所周知的原因，被封了（偶尔也会能上）。唯一可惜的是，上面有很多不错的统计资源。

软件类：

• Google 输入法：拼音输入法，表现中规中矩，可以同步用户词典。
• Goolge 浏览器：感觉和 firefox 差不多，当然远远比 Internet explore 好很多很多辈（强调一下——不是倍）。
• Google 词霸：自从有个这个，我就不再买正版的金山词霸。配合 neospeech 的 TTS （text-to-speech），挺舒服 ^_^
• Google Earth：这个就更不用说了吧，我能通过它找到我家屋顶。

我这里只是简单了列了一下常用的 Google 服务，像一些比如 Google 学术搜索、Google 生活搜索、地图、桌面之类的我都没有提到，但它们都在或多或少的影响着我们的生活。

还是那句话——我患有严重的 Google 依赖综合征。

最近《阿凡达》应该是大陆电影市场最为火爆的电影了，媒体（包括博客）铺天盖地的报道，终于让我没按耐住性子。不过想到大早起排队去买 imax+3D 的痛苦，只好退而求其次去看 3D 版。

在影院看的时候还是非常震撼（虽然架着两个眼睛很累）：

• 阿凡达讲述的故事依旧很俗套，正义战胜邪恶，男女主人公收获爱情。
• 特效相当赞，虽然没有感受 imax 版本，但普通 3d 效果让人印象非常深刻
• 召集外星部族抵抗地球人的那段演讲，让我想起了《勇敢的心》里华莱士那段

为什么《阿凡达》这么火，网上诸多评论。有人说是今年影迷对国内电影的绝望（三枪），也有上升到综合国力的问题（冯导的）。我也凑个热闹，跟风总结一下不靠谱的：

• 《阿凡达》实际上是一部日本（美国）动漫的电影版：好像依稀记得小时候看过的《太空堡垒》这部漫画，男主人公最后就是爱上了一个女性外星人。
• 地球人是为了一种“磁悬浮矿石”才到潘多拉行星的，噢，这是传说中的“沙丘”的引子……
• 《阿凡达》讲述了一部外星人原住民，即钉子户抵抗暴力拆迁的故事……

最后一个更为不靠谱，看完之后我一直在怀疑卡梅隆卡导是不是魔兽争霸的粉丝，要么就是收受了暴雪公司的贿赂，替人家做商业宣传。这里总结一下国内某著名w3论坛中对《阿凡达》的故事情节叙述：

刚开始，人族（地球人）经济科技都有绝对的优势，依靠 100 人口高攻防的飞机（船）、坦克、火枪手打爆了暗夜精灵的生命古树。但是在骑乘着奇美拉的恶魔猎手的带领下，暗夜精灵集结了 100 人口的角鹰骑士+女猎+弓箭手的立体混合兵种，但因没有升级攻防而处于下风。关键时刻，由于人族不小心勾到重甲野怪，导致大量中立生物攻击人族。在暗夜大军和中立生物的联手攻击下，人族无奈溃败。

最终英雄的单挑结局是：虽然拥有 6 级技能——天神下凡的山丘之王打败了同样早早拥有 6 级变身能力的大恶魔，但是由于厮杀过度，没血没魔，最后被一个红血的弓箭手连射两箭给搞定。

最后还附了一下卡梅隆导演打魔兽的一些技巧感想：

• 最后决战的时候千万要小心，别打到中立生物。
• 角鹰骑士的攻击值太低，打不过直升机；但是角鹰能打过。
• 女猎打不过火枪兵。
• 在重甲野怪面前，火枪不堪一击。
• 骑着奇美拉的恶魔猎手非常拉风，估计卡导肯定用的这个游戏bug......

前一段时间，John Baez 在自己的主页上更新了一篇文章名为he Beauty of Roots ，这篇文章之后在“科学松鼠会”上被《多项式的根之美》 转载。上面提到了曼德布洛特集，根据其发明者法国数学家Benoît Mandelbrot 而命名。

曼德布洛特集是一种分形，从一般分形性质来说：

客观自然界中许多事物，具有自相似的“层次”结构，在理想情况下，甚至具有无穷层次。适当的放大或缩小几何尺寸，整个结构并不改变。不少复杂的物理现象，背后就是反映着这类层次结构的分形几何学。

常见的曼德布洛特集是这个样子（分辨率原因，部分细节显示不够）：

假如我们把这个集合的下半部分（最下边的小块）分割出来，就是这个样子（8倍放大）：

由于分辨率的提高，所以显示了第一幅图中并没有显示的细节。

继续放大，上图的左上部分的那个小枝（6倍放大）：

再把上图最靠近左边的那个小枝——放大（50/3倍放大）：

继续放大最左边的小枝，似乎在末端又出现了一个类似的小枝（5倍放大）：

如果继续放大下去可能还是这个样子 ：）

注释：

1. 最后一张图相比第一张图来说相当于局部放大了 4000 倍。
2. 高质量的矢量绘图数据量比较大，R 处理起来有些问题，只好使用局部放大的方式。
3. 更多的使用其他软件绘制图形可以见：http://commons.wikimedia.org/wiki/Mandelbrot_set

附代码：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108

setwd("D:\\doc\\Mandelbrot\\pic")
# Another neat animation of Mandelbrot Set
jet.colors <-
  colorRampPalette(
    c(
      "#00007F",
      "blue",
      "#007FFF",
      "cyan",
      "#7FFF7F",
      "yellow",
      "#FF7F00",
      "red",
      "#7F0000"
    )
  ) # define "jet" palette
m = 400
C = complex(real = rep(seq(-1.8, 0.6, length.out = m), each = m),
            imag = rep(seq(-1.2, 1.2, length.out = m), each = m))
C = matrix(C, m, m)
Z = 0
for (j in 1:20) {
  Z = Z ^ 2 + C
  X = exp(-abs(Z))
}
png("scale1.png")
par(mar = c(0, 0, 0, 0))
image(X, col = jet.colors(1000))
dev.off()



m = 400
C = complex(real = rep(seq(-1.6, -1.3, length.out = m), each = m),
            imag = rep(seq(-0.15, 0.15, length.out = m),      m))
C = matrix(C, m, m)
Z = 0
for (j in 1:20) {
  Z = Z ^ 2 + C
  X = exp(-abs(Z))
}
png("scale2.png")
par(mar = c(0, 0, 0, 0))
image(X, col = jet.colors(1000))
dev.off()


####

m = 400
C = complex(real = rep(seq(-1.4, -1.3, length.out = m), each = m),
            imag = rep(seq(-0.15, -0.05, length.out = m),      m))
C = matrix(C, m, m)
Z = 0
for (j in 1:20) {
  Z = Z ^ 2 + C
  X = exp(-abs(Z))
}
png("scale3.png")
par(mar = c(0, 0, 0, 0))
image(X, col = jet.colors(1000))
dev.off()


m = 400
C = complex(real = rep(seq(-1.35, -1.30, length.out = m), each = m),
            imag = rep(seq(-0.13, -0.08, length.out = m),      m))
C = matrix(C, m, m)
Z = 0
for (j in 1:20) {
  Z = Z ^ 2 + C
  X = exp(-abs(Z))
}
png("scale4.png")
par(mar = c(0, 0, 0, 0))
image(X, col = jet.colors(1000))
dev.off()


m = 400
C = complex(real = rep(seq(-1.3280, -1.3250, length.out = m), each = m),
            imag = rep(seq(-0.1225, -0.1195, length.out = m),      m))
C = matrix(C, m, m)
Z = 0
for (j in 1:20) {
  Z = Z ^ 2 + C
  X = exp(-abs(Z))
}
png("scale5.png")
par(mar = c(0, 0, 0, 0))
image(X, col = jet.colors(1000))
dev.off()


m = 400
C = complex(real = rep(seq(-1.3276, -1.3270, length.out = m), each = m),
            imag = rep(seq(-0.1221, -0.1215, length.out = m),      m))
C = matrix(C, m, m)
Z = 0
for (j in 1:20) {
  Z = Z ^ 2 + C
  X = exp(-abs(Z))
}
png("scale6.png")
par(mar = c(0, 0, 0, 0))
image(X, col = jet.colors(1000))
dev.off()