大自然的力量永远让人敬畏，地震、海啸、陨石、雷击……因为我们在大自然的面前如此渺小，所以我们必须协作，必须发展科技，必须懂统计 ：）
说到地震，我可能比较敏感，因为我是唐山人。虽然 76 年唐山大地震时，我还没有形成碳水化合物形态，但后来，每每听到老一辈讲起当时的惨烈，心常戚戚。
讲几则作为唐山人的小故事：
- 小时候对地震的初相识：有一次刚刚从床上爬起来，正在懒洋洋的坐着穿衣服，就发觉床开始做规则的前后晃动。当时年龄小，很无知，只知道傻乎乎地沉浸在如秋千般的跳动中，那叫个带劲……
- 有一次上课，感觉课桌在晃动，没法看书，于是停下来查看周围哪个同学在晃。检查一圈发现周围几个同学没有一个再晃！结果，紧张地直接拍案而起，大呼——地震啦（事后新闻证明是真的）！在我们那儿，这点比较好：即使是课上误判地震的这种事儿，一般老师都是笑笑而过 ；）
- 每年我们那都会有地震的谣言，而且说的神乎其神，俺老爸一般会守夜（感谢俺老爸！）。或者天气好的话，干脆去广场之类空旷的地方，找地方打地铺。当然一般都是打牌、聊天到 24 点，然后回家睡觉。
恩，不多扯了，言归正传。自从汶川大地震以后，国人对地震明显敏感很多。且不说海地，单单前两天（24日）山西河津、运城地震就让然琢磨不懂：有人说，21日山西省地震局辟谣称不会有地震，但运城震感明显。为什么地震局会出来辟谣，仔细一读，原来才知——地震局指的是“破坏性”地震。
但有个问题：
国务院1995年颁布的《破坏性地震应急条例》，破坏性地震指“造成一定数量的人员伤亡和经济损失的地震事件”， 并没有规定特定的级数。
这破坏级地震可不是闹着玩的，得仔细瞧瞧最近这地震都发生在哪里了，震级多大？是不是会对我们构成威胁！？于是有了下面这张图——最近一周中国及周边版图地震情况（1月20日至1月25日共计六天）：
数据童鞋们可以在这里查看，里面的震级需要注意一下，有Ms和ML两种，换算关系如下。但具体什么意思大家直接 wiki 好了。
ml=(1.17mb+0．67)/1.13 ml=(ms+1.08)/1.13
一些说明（不是写商业报告，偷工减料啦）：
蓝色的背景是地震点的密度——也许是喜马拉雅造山运动，也许是三峡工程，不管怎样，四川地区不太平啊！弟兄们小心！
红色的点代表地震的位置，其大小表示震级的大小。
从1月20日至1月25日，版图周边共计有901条地震记录（有点吓人）！其中大于ML5级的一共两次：
2010-01-24 10:36:13.8 35.45 110.70 15 Ms4.8 天然地震 山西河津
2010-01-21 10:02:02.8 13.70 125.85 33 Ms5.1 天然地震 菲律宾群岛地区
最后我们再回头看一下，最近一周地震的震级（ML）分布：
至少可以长舒一口气，原来大部分都是小震，不具“破坏性”的居多。
有点标题党的嫌疑，实际是介绍如何使用 R 绘制 heatmap 的文章。
今天无意间在Flowingdata看到一篇关于如何使用 R 来做 heatmap 的文章（请移步到这里）。虽然 heatmap 只是 R 中一个很普通的图形函数，但这个例子使用了2008-2009赛季 NBA 50个顶级球员数据做了一个极佳的演示，效果非常不错。对 R 大致了解的童鞋可以直接在 R console 上敲
？heatmap
直接查看帮助即可。
没有接触过 R 的童鞋继续围观，下面会仔细介绍如何使用 R 实现 NBA 50位顶级球员指标表现热图：
关于 heatmap，中文一般翻译为“热图”，其统计意义wiki上解释的很清楚：
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.
下面这个图即是Flowingdata用一些 R 函数对2008-2009 赛季NBA 50名顶级球员指标做的一个热图（点击参看大图）：
先解释一下数据：
这里共列举了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 同样更新了博客。唉，这厮嗅觉也忒灵敏！哈哈
R 各个镜像中的 Contributed Packages 越来越多，截至今日，已经达到1950个，单单拉动鼠标把所有的 包名 从 A 到 Z 过一遍也得 10 几秒。随便考你一道：最后一个 R 包是啥？
zoo？
呵呵，我的印象里一直是它，仔细瞧了瞧发现是个叫 zyp 的包。
又一次领略了 R 强大的扩展能力撒？这个特点给我们带来了一些烦恼，因为人类的大脑能够理解的概念是有限的，对于没有任何关联的概念，我们的识别能力一般不超过 7，而且 R 的涵盖范围实在太广。从我们的有限性（7个概念）和 R 的无限性这一角度讲，逐一认识这些包几乎是不可能的！不过还好，至少我们可以可以参考 CRAN 上的 Task Views，大致了解 R 包的使用方向。
我们换个思路，不是从 R 的使用方向上，而是从 R 包的依赖关系上？
这些 R 包并不是相互独立的。比如说，MASS 包依赖于 R (>= 2.5.0), grDevices, graphics, stats, utils 这些基础包；而又会有包依赖于 MASS 包，比如 yihui 的 animation ，当然还有可能有包依赖于 animation ……
遍历所有的包，我们就看到了一个网络，一个 R 包的网络。
为了简化起见，这里忽略了同其他包没有关系的包（当然并不是完全没有关系，所有的包都和 R 或 R 的基础包有关，如果这样计量的话，会导致所有的包都会指向 R）。
首先抽取了这个庞大网络的一部分：
从上图我们可以看到，标记点为215、271的两个包是我们研究的包网络中的两个关键点，这两个包分别是lattice、mvtnorm。
关于这两个包：
-
lattice：网格绘图的基础包。很多包基于它扩展并不惊讶吧；
-
mvtnorm：多元正态分布和t分布的概率密度函数、累计分布函数、分位数函数、分布随机数。多元分布的基础。
从 271（mvtnorm）向左上，又会有一个小的聚集。那个小的聚集中心（110），是 fBasics 包，如果各位对金融领域关注的话，应该知道它在其中的地位吧。
当然，由于抽取的是一个子网络，很多的连接都被生硬地隔断，因此出现了大量的孤立点。
如果我们把 CRAN 上的1950个包都放到我们的网络中会是这样：
1、第一张图的 包 id 换成 包名称 会导致 演示的视觉效果很差，网页又不支持 pdf 直接显示，只好把带包名的图放这（pdf）。
2、带包名的 ，1950 个包的全图就算了吧，单绘图就得 2 分钟，更别提调整参数了 ……
还记得第一次看到水立方时的惊讶么？
是什么这么吸引我们？是有如天空般的颜色？还是那气泡似的形状？
从水立方的外墙结构上看，不但外观美观，而且十分紧凑。水立方外墙为什么会有这样的性质，是因为它上应用了一项最优化的技术，即Voronoi 原理。
Voronoi 图也常常被称为 Dirichlet 格局（Dirichlet tessellation）。通俗讲，其原理是一项从点到面的技术。它的每个多边形只有一个”生成点”，而这个多边形上的每个点到”生成点”的距离总是比到其他”生成点”的距离要小（是不是想到了 K-means 算法？）。
在建筑设计上，有设计人员争论这类方法定义为“参数化设计”。认为这种方法不能同传统的、依靠灵感的设计方式相比，因为这种方法高度依赖计算机，只需要简单改变若干参数就能得到设计方案。但这个论断，恰恰忽略了“参数化设计”背后的数学意义。
既然 Voronoi 是一种最优化的算法，那么除在建筑学上给我们带来的美轮美奂的视觉效果外，它在空间统计上同样也有应用。
下面，我根据各个省会城市（包括香港、澳门）的地理位置，利用 Voronoi 原理，计算每个省最佳控制范围（使用红色的线条标记）：
虽然理论值（最优）和现实值（行政区划、地理）总有差距，但是，比较一下会发现一些值得探讨的现象：
-
内蒙古应该好好的规划一下，从东边到西边实在太远了，把西边的划给宁夏可能好点；东边划给北京、东三省；
-
河北北部，不论是属于北京还是天津都会好些，记得我小的时候，宁可去北京、天津，也不乐意去遥远的省会–石家庄；
-
青海应该把甘肃的北部包括进去；
-
上海、香港、澳门有一部分管辖区也也不错么：）
整体上看，大部分省的行政区划还是符合 Voronoi 原理。也就是说，单纯从空间距离的角度来看，我国的行政区划还是比较不错的。，
## the R code ## library(diagram) pos <- coordinates(pos=c(3,3,3,3,3,3)) cc <- c("Start",LETTERS[2:16],"End") openplotmat(main="Flow Chart") for(i in seq(1,15,by = 3)) straightarrow(from = pos[i,] ,to = pos[i+3,]) for(i in c(2,5,8)) straightarrow(from = pos[i,] ,to = pos[i+6,]) segmentarrow(from = pos[16,],to=pos[2,],path="RVL",dd=0.15) bentarrow(from = pos[8,], to=pos[6,],path='H') bentarrow(from = pos[6,], to=pos[2,],path='V') straightarrow(from=pos[14,],to = pos[17,]) for(i in c(2,7,13,14,16)) textrect(pos[i,],radx=0.08,rady=0.04,lab = cc[i]) for(i in c(1,17)) textround(pos[i,],radx=0.08,rady=0.04,lab = cc[i]) textdiamond(mid = pos[8,],radx = 0.15,rady = 0.08,lab = c("Has","Detect?")) textmulti(mid = pos[4,],radx=0.1,rady=0.05,nr=6) textmulti(mid = pos[6,],radx=0.1,rady=0.05,nr=5) textellipse(mid = pos[10,],radx = 0.1,rady = 0.05) text(pos[8,1]+0.2,pos[8,2]+0.03,"YES",cex = 0.8) text(pos[11,1]+0.04,pos[11,2],"NO",cex = 0.8)
以前听歌的时候，一直认为一首 mp3 大概就是3、4M的样子。今天突然想：
到底我常听的歌，大概都有多大呢？说做就做——用 R 跑个结果：
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.7019 3.3540 4.1090 4.4320 5.2390 14.3500
应该是4M左右.
