之前,corrplot包(部分效果见此)只能通过Rforge下载:

install.packages("corrplot", repos="http://R-Forge.R-project.org")

目前小bug都找的差不多了,加上近来比较忙碌,故打算提交到CRAN(大约五天之内会到CRAN上露脸吧),需要此包的朋友们就不用发email给我原始数据让我代劳了

此包以后的更新方向主要是变量的重排序方法:

1. Robinsonian
2. Dimension reduction
3. Heuristics
4. Block modeling
5. TSP

现在已经实现了主成分排序和各种系统聚类排序,其他的还得边学边卖,慢慢更新。相关矩阵可视化竟然能扯出这么多数学、统计甚至图论的东西,之前从没想到过,真是好玩。

注1:最初是在R会议上看见bjt大哥用椭圆图来表示相关矩阵,那时觉得很新鲜、很好玩,记忆很深刻。后来随便想了一阵子,写了个小函数来娱乐,却没想到滚雪球滚成了一个小package。

注2:曾经觉得自己折腾得太久了,很无聊,不过现在又觉得很好玩了,因为还有很多有趣的工作要做。

注3:corrplot包在Rforge上最近不太好用,等我忙完手头的事立即更新。

 
前两周在北大上可视化的暑期班,有幸和五湖四海的朋友们一起聆听Kwan-Liu Ma、Han-Wei Shen、Alex Pang、Michelle Zhou、Hua-min Qu、Jean-Daniel Fekete、Jian Huang、田捷等老师的教诲,这些老师、研究人员在各自的领域内都非常优秀,部分还是界内大牛,更可敬的是他们都对学生很有感情、很有耐心——标准的德艺双馨。

整个暑期班的学习中,课程包含流体可视化、张量可视化、医学影像、信息可视化、时变可视化、智能可视化、并行可视化等很多方面,其中我最感兴趣的是:Jean-Daniel的Visualizing Social Networks using Hybrid Matrix/Node-Link Representations,因为和我之前的工作颇有渊源。

Jean-Daniel在做social network的时候,用到了类似相关矩阵可视化的东西,就是将两两之间的关系数字化,得到一个相似度矩阵,然后可视化这个矩阵。social network的传统做法是画个网络图,用节点和连线来表示,但这样很容易使整个图变得乱七八糟,什么也看不清。可视化相似度矩阵的方法则不存在这个问题,当然也会带来新的麻烦。

对于相似度矩阵的可视化,主要存在两个问题:
1. 如何用颜色、图形、线条表示这个矩阵,;
2. 如何对矩阵对应的变量进行重排序,使得相似的变量聚在一起,不相似的分开,这样我们可以通过可视化的图形直观地发掘变量之间内在的关系。

其中,第一点已经很成熟了,就是用方块、圆,再辅之以渐变色等,corrplot包的初期也就是做这些工作;而第二点,即如何重现排序变量,这会涉及到统计、数学知识,也是本问题的精髓之所在。之前,我仅知道用PCA、聚类等方法重排变量,也在corrplot包中实现了。而现在,我发现重排变量是个不小的问题,因为它本身非常重要,而PCA、聚类方法有时在效果或者速度上并不占优势,这就需要我们探索其他方法。Jean-Daniel在课堂上介绍了两个, Robinsonian和TSP。Robinsonian是个很数学的东西,我暂时还没有翻看论文,但是TSP(Travelling salesman problem)大家都再熟悉不过了,把这个东西灵活地用在变量排序中,的确是别出心裁,匪夷所思!!相似度、相关系数等本身都是距离,而TSP问题恰恰求最短路的。

TSP是个NP问题,但很幸运的是,我们目前已经有很多算法可以快速得到不错的解,R中有相关的包(TSP),包含了常见算法并提供了concorde软件(解TSP问题的优秀开源软件)的接口。这样一来,写个用TSP排列变量函数就方便很多了。

当然,除了TSP问题的求解难度问题之外,它在变量排序中还存在一个问题,就是TSP问题求出的最短路是个环线,所以在重排序的变量中,第一个和最后一个可能很相似,但在图中,它们一个最上、一个最下,离得最远。这个问题可以这么解决:

1. 不在一张图上吊死。用至少两张图,第二张图的变量顺序是第一个图的水平移动,比如第一个图中(A, B, C, D, E),而第二个图则是(D, E, A, B, C),这样第二个图中E和A就在一起了。当然,我们也可以通过观察图形,得到一个最容易接受的排序。这虽然是个解决方法,但是人总是贪婪而又懒惰的,一张图能看清楚的,绝不会看两张图,因此还需要探索一张图的方法。

2. 不在一个算法上吊死。既然TSP可以,那么图论中经典的Dijkstra 、Floyd算法也很可能适用,虽然这两个算法不是穷尽各个节点的,而是求各个节点之间的最小路程。比如,我们可以通过这两个算法辅助TSP算法确定起点和终点:我们可以求出网络中任意两点间的距离,然后找出最大距离所对应的两个节点。然后,将距离矩阵中这两个节点所对应的距离修改为0,这样得到的结果中这两个节点肯定挨在一起,这相当于将TSP环路算法转换为非环路算法。然后,将这两个节点分别设置为头尾,就可以得到一个粗糙的结果了。我在R中试了试,基于经典的mtcars数据,将得到的图展示出来:
vis-tsp-pca
从上图可以看到,两种排序方法还是比较相似的,并且效果都不错。这样确定起点终点的好处是,起点和终点对应线路是所有两两路程中最长的,这样再用非环线的TSP算法就不容易使排序失去意义。当然,这种方法还是很粗糙,比如计算量过大。实际中,我们可以通过别的方法更快地确定起始点和终点。

3. 不在一种介质上吊死。常见的纸、屏幕是平面的,如果我们有圆柱形式的立体介质,那么TSP得到的变量重排序就很有舞台了。将图绘制在圆柱上,首尾相接,看的时候转动圆柱即可。这个方法听起来的确有些扯,但是我觉得这种介质的出现不是没有可能(实际上,一些路边的广告就是这样的),当然这种方法的局限性也很大。

等手头的杂事忙完的话,corrplot也会逐步更新,添加一些变量排序的新方法。可视化不是简单的画图,背后的算法、模型非常重要

注:
1. 已经有很多文献讨论了矩阵的重排序,不过我都没看,先自己折腾一番。
2. 本文之前写得不太明了,因此重新修改了,2009-08-27,19:16。

 

据楚天金报7月31日报道,前日,老河口市第二期经济适用房公开摇号,从1,138名具有资格的申请人中摇出了514名住户。当晚,老河口市民发现这514户中出现了“14连号”的现象。

很多网友都算了这个概率,我也凑个热闹。首先计算仅仅出现14连号的概率,不妨先如此定义:1.不限制出现14连号的串数;2. 不得出现15连号。

下面的式子是把 1125种可能的14连号分为两类(两头的和中间的)再粗略计算。其中分子上的1123是指1138-14-1(不选两头14连号挨着的1个数字,比如选了1-14,就不选15),1122是指1138-14-2(不选中间14连号左右挨着的2个数字,比如选了3-16,就不选2和17),这样可以减少重复和15++的连号,但远不是精确解,这是我见到本题的第一想法。

0.004624725

这仍然是个比较粗略的值,有不少重复事件以及15++连号事件算了进去,但是对结果影响似乎不是很大,同计算机模拟的结果差不多,可能是和真值比较接近吧。模拟代码如下(R 环境,更多代码参见我在cos统计之都论坛的帖子),模拟结果是 0.0045 左右:

## 该函数模拟出现"cont连号"的概率,但不限串数
function(n = 50000, cont = 14){
  for(i in 1:n){
    nums <- sort(sample(1138,514))
    nums1 <- nums[1:(514 - cont + 1)]
    nums2 <- nums[cont:514]
    flag <- nums2-nums1
    if(any(flag==(cont-1))) x <- x + 1 #这里不区分14连号的串数

    nums3 <- nums[1:(514-cont)]
    nums4 <- nums[(cont + 1):514]
    flag <- nums4-nums3
    if(any(flag==cont)) x <- x - 1 #去掉15++连号事件
  }
  x/n
}
cont()

模拟14++连号(定义为至少出现了14连号,可能还有15、16或者更多数目的连号,且不限串数)的概率,结果约为0.0082左右,这个结果很多网友都得出来了:

## 该函数模拟出现"con++连号"的概率
con <- function(n = 50000, con = 14){
  x <- 0
  for(i in 1:n){
    nums <- sort(sample(1138,514))
    nums1 <- nums[1:(514-con + 1)]
    nums2 <- nums[con:514]
    flag <- nums2-nums1
    if(any(flag==(con-1))) x <- x + 1
  }
  x/n
}
con()

模拟单串14连号出现的概率,结果和不限串数的14连号概率很接近,都在0.0045左右,代码如下:

## 该函数模拟仅仅出现1串"con连号"的概率
con2 <- function(n = 50000, con = 14){
  x <- 0
  for(i in 1:n){
    nums <- sort(sample(1138,514))
    nums1 <- nums[1:(514-con + 1)]
    nums2 <- nums[con:514]
    flag <- nums2-nums1
    if(sum(flag==(con-1))==1) x <- x + 1
  }
  x/n
}
con2()

此外,有网友曾经得到了下面更粗糙的解(我不确定他是如何定义问题的,但显然不论单串还是多串,不论P(14)还是P(14++),其结果都粗糙了点,因为重复事件太多了),尽管式子和上面的式子很像。算是差之毫厘,失之千里吧。

0.01499478

如果出现14连号事件(不包括15及以上连号事件)的正确答案子0.45%附近,算是个稍微小了一点的概率事件,我们顶多可以怀疑(这样说主要是心理因素而已),但是不能因此断定他们的随机抽样是做了手脚的。因为小概率事件是完全可能发生的,并且实验只有一次,任何一组号都是独一无二的,都是小概率事件。记得有个买体彩P3的朋友坚信出号机有规律,告诉我888三个号很不容易出现,但是我查了历史,这个号照样出现了,并且出现的频率很符合统计规律。我不知道为什么他说不会出现,但是似乎在888出现的时候人们的反响比较大,其中不乏一些质疑的声音。

在1138个数字中随机选择514个,任何一组数字出现的概率都是:

0.008336363

概率都是非常小的,只不过本次的数字有14连号,引起了大家的注意,我在一些论坛上看见了很多荒唐的评价,很多质疑都是没有根据的谩骂。实际上,即使抽出的是一堆杂乱的数字,我们也不能说抽样没有做手脚。

怎么来证明他们的随机抽样有没有做手脚?我认为在这里单纯算算概率仅仅可以让我们心神不安地怀疑一下,但根本不足以下结论。关键是要从常理出发,并且考虑到经济适用房的销售模式、历史经验,然后客观地进行分析,最起码要考虑这几点:

1. 申请人的ID是怎么确立的,是随机的吗?还是按照什么规律。
2. 这14个申请人之间有什么联系?是不是相互认识,是不是有相同的后台或者单位?他们个人或者所在的单位和负责分配房子的机构有没有值得怀疑的来往?这一点尤其关键。
3. 随机抽样的过程是不是公开的?采用什么样的技术和设备?是不是有信得过的监管机构?
4. 负责分配房子的机构里有哪些人?他们有没有违法违纪或者不得民心的前科?
5. 不要把注意力仅仅放在这14连号之上,如果抽样有猫腻,其他500个申请人也必须仔细盘查。这一点很容易被人忽略,但是却是非常关键的一点。
6. ……

要调查的可能还有很多很多,总之,看问题要以系统的眼光来看,俯瞰全局,任何和该系统有关的东西都不应该放过。而在千头万绪之中,还要善于抓住事物的主要矛盾,去粗取精、去伪存真。

统计永远不能告诉我们所有的答案,它可以给我们一个思路、一个方法或者一些灵感,但不能替代对事物本身的机理分析。工具永远只是工具,尽管它的确很重要,但它不是问题本身。我们的最终目的不是炫耀工具,而是解决问题。当然,在本案中,不存在炫耀工具的问题,也没有什么技术性的东西。

我们应该相信科学,而不是凭空觉得这个概率比月亮撞上地球的概率还小;也希望某些试图仅仅通过统计手段就想得到结论的朋友们不要太“唯统计论”。写到这里,突然想起今天见到的阔别三年的朋友的一句话:任何东西,只要犯上个“唯”字,就很容易走极端而出错。诚然,任何手段任何理论任何观点,都有它本身的局限性,没有放之四海而皆准的真理。这也就是哲学中的矛盾论所阐述的精华了。

© 2010 优秀是一种习惯 taiyun.wei@cos.name Suffusion theme by Sayontan Sinha