之前,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。

 
发现第一届 R 中国会议上了The R Journal第一期(Conference Review: The 1st Chinese R Conference,第69页),很有历史意义,谢兄辛苦了啊。

此外,还发现文中还介绍各个演讲者和对应的演讲课题,我的也放了进去(Statistical Computation :optimization in R by Taiyun Wei),这是我的名字第一次上正式学术期刊,呵呵。
不知第二届 R 中国会议何时何地举办,消息灵通的看客了不妨透漏一下

 
下午拿 R 去忽悠师兄师姐们,主题是回归和分类在R中的实现,讲得乱七八糟,幸好许青松老师宽宏大量,师兄师姐们都很给面子。不过可能由于我耿耿于怀没有讲好,晚上吃饭时闹了个笑话:
晚上去吃蒸菜,愣是要份回归肉(这两天回归树、回归机说得多了些,可怎么又冒出了个回归肉啊),搞得服务员摸不着头脑,整个半天才弄明白我要的是梅菜扣肉(不是回锅肉)。
 
截一个游戏结束时的图如下:

slide-puzzle

目前该函数已被fun包收录:

install.packages(“fun”, repos=”http://R-Forge.R-project.org”)


代码如下:

function (size = NULL, bg = "lightblue", z = NULL, ...)
{
    if (!is.null(size)) {
        n <- size[1]
        m <- size[2]
    }
    if (!is.null(z)) {
        n <- dim(z)[1]
        m <- dim(z)[2]
        if (!is.null(size))
            warning("Because \"z\" is specialized, parameter \"size\" will be omitted.")
    }
    z.right <- matrix(1:(n * m), n, byrow = TRUE)
    z.right[n, m] <- 0
    neg_seq.length <- function(x) {
        len <- 0
        for (i in 1:(length(x) - 1)) {
            tmp <- x[(i + 1):length(x)] - x[i]
            len <- len + sum(tmp < 0)
        }
    }
    len.right <- neg_seq.length(as.vector(z.right)) + n + m
    if (is.null(z))
        z <- matrix(sample(0:(n * m - 1)), n)
    else {
        len.z <- neg_seq.length(as.vector(z)) + sum(which(z ==
            0, arr.ind = TRUE))
        if ((len.right%%2) != (len.z%%2))
            stop("The sliding puzzle is insoluble!")
    }
    len.z <- neg_seq.length(as.vector(z)) + sum(which(z == 0,
        arr.ind = TRUE))
    while ((len.right%%2) != (len.z%%2) | all(z == z.right)) {
        z <- matrix(sample(0:(n * m - 1)), n)
        len.z <- neg_seq.length(as.vector(z)) + sum(which(z ==
            0, arr.ind = TRUE))
    }
    z[!z] <- NA
    replot <- function(z) {
        bg <- ifelse(z, bg, "white")
        fg <- ifelse(z, bg, "white")
        par(mar = c(0, 0, 0, 0), bg = "white")
        plot(c(0, m), c(0, n), type = "n", axes = FALSE, asp = 1,
            xlab = "", ylab = "")
        segments(0:m, rep(0, m + 1), 0:m, rep(n, m + 1), col = "grey",
            lwd = 2)
        segments(rep(0, n + 1), 0:n, rep(m, n + 1), 0:n, col = "grey",
            lwd = 2)
        symbols(0.5 + rep(0:(m - 1), each = n), 0.5 + rep((n -
            1):0, m), squares = rep(0.9, n * m), add = TRUE,
            inches = FALSE, fg = as.vector(fg), bg = as.vector(bg))
        text(0.5 + rep(0:(m - 1), each = n), 0.5 + rep((n - 1):0,
            m), as.vector(z), cex = 3)
    }
    push <- function(x, begin, space) {
        tmp <- x[space]
        if (begin < space) {
            x[(begin + 1):space] <- x[begin:(space - 1)]
            x[begin] <- tmp
        }
        if (begin > space) {
            x[(begin - 1):space] <- x[begin:(space + 1)]
            x[begin] <- tmp
        }
        x
    }
    count <- 0
    mousedown <- function(buttons, x, y) {
        plx <- grconvertX(x, "ndc", "user")
        ply <- grconvertY(y, "ndc", "user")
        m.col <- ceiling(plx)
        m.row <- n - floor(ply)
        ind.NA <- which(is.na(z), arr.ind = TRUE)
        if (!xor(m.row == ind.NA[1], m.col == ind.NA[2]))
            cat("Warning: Cannot push any number!\n")
        ind.NA <- which(is.na(z), arr.ind = TRUE)
        if (ind.NA[1] == m.row & ind.NA[2] != m.col) {
            z[m.row, ] <<- push(z[m.row, ], m.col, ind.NA[2])
            cat("step = ", count <<- count + 1, "\n")
        }
        if (ind.NA[1] != m.row & ind.NA[2] == m.col) {
            z[, m.col] <<- push(z[, m.col], m.row, ind.NA[1])
            cat("step = ", count <<- count + 1, "\n")
        }
        replot(z)
        flag <- z == z.right
        if (all(flag[!is.na(flag)])) {
            paste("You win! Time:", round((proc.time() - ptm)[3],
                2), "seconds.")
        }
    }
    ptm <- proc.time()
    windows(5, 5)
    replot(z)
    if (.Platform$OS.type == "windows")
        getGraphicsEvent("Game begin!", onMouseDown = mousedown)
}
 
R有个sudoku包,不过有些粗糙,不能保证产生的数独有唯一解(看看源代码就知道它的生成机理,对一个已知的数独做变换而已),对于多解问题也只是仅仅得到一个结果(其实一般数独不允许多解)。下面是玩数独时一个问题的正确结果,但是软件包判错:
sudo-bug

其实本题至少有两个正确结果,sudoku包的结果:

sudo-ok
 

花了很长时间,终于搞定了《R软件与最优化》的初稿。

包含内容:线性规划、整数规划、目标规划、非线性规划、图与网络规划五大块,对于指派问题、运输问题、TSP问题、最大流、最短路等进行了专题讨论。

先说说LaTeX。这是我第一次拿LaTeX软件写正式东西,以前一直听人家说LaTeX如何如何强大,这次可在这强大的软件上栽跟头了,插个三线表、排个数学公式之类的问题就让我叫苦不迭。更可恶的是为了看看效果,稍稍写点东西就得编译PDF文件,讨厌的黑DOS窗口哗啦啦哗啦啦,好久才能得到PDF文件,有时还得运行两三遍。然后,把原来写错的改一改,再写,再哗啦啦……唉,还是所见即所得用着舒服,不由想起了好用的office,金山个人免费版就不错,才几十M,什么公式、表格、图形的自动编号、自动更新都可以轻松搞定,数学公式也很容易插入,并随时可以看见文章的效果。当然LaTeX的最终效果还是好看,但愿以后用它能够轻车驾熟。

再说说R。当初打算写R软件和数学建模,是因为暑假里做数模用过一点R。不过大都是画画图、分析点数据,至于优化问题,都是队友用LINGO或MATLAB搞定的。而我用R做的那点统计实在是太简单了,怎么敢好意思写给别人看。数模中,最优化是一大块内容,是很重要的一类问题,加之R有一系列关于最优化的包,因此,我就选了最优化这一块。于是乎,我将最优化的包几乎全下载下来,并翻译了一个最优化包的简单介绍,然后就刀耕火种了:一个一个尝试,将自己认为比较好用(可能别人不觉得好用)的写在文章里。遇到最严重的问题就是语言问题,我的英文极差,高考后就江河日下,而R的文档全是英文的,我只能硬着头皮看了,当然有些不用看结合以下问题就知道是什么东西。朋友介绍我用linges全文翻译,试了一下效果很差,再此绝望。看来,英文又得捡起来好好学学……

现在终于基本搞定了,就27页,当初做的时候还以为能出本书呢,不过要是把理论加上去,多写几个案例,页数很容易就上去了,只是我时间有限,还得忙统计正事呢。希望数模元老人物薛毅先生以后能写写R和最优化了。

最后要说的是,R在处理最优化问题绝对比LINGO好,很多问题都有专门的包,对症下药,解决问题自然是事半功倍。用LINGO的数模队友看到R两三句代码就轻松搞定LINGO几十行代码才能勉强解决的问题时,感到非常惊讶。其实岂止如此,好些R轻松解决的问题LINGO是无法解决的,比如画个地图、画个网络图等,因为LINGO画图很差。当然,画漂亮的图不是R的专利,LINGO可以去开发,但是LINGO公司就那么多科研人员,还要保住自己的优势,哪有余力涉及太广的领域?但是R就不同,你擅长这个,捐个这方面的包,他擅长那个,捐个那方面的包,如此一来,R就汇集了百家之所长,成为一个强大的平台。

牢骚好久了,打住,好好学习……

p s.需要本文的朋友请移步: https://github.com/taiyun/Optimization-using-R

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