一个简单的马尔可夫链

在一个随机过程中，如果事件发生概率在t时刻所处的状态为已知时，它在t + 1时刻只与t时刻的状态有关，而与之前所处的状态无关，则称该过程具有马尔可夫性。

时间和状态都是离散的马尔可夫过程称为马尔可夫链。马尔可夫链在经济学，社会学，生命科学领域有着广泛的应用。这里举例说明。

例子： 姜华平、陈海泳对某城市2002年居民出行方式所占比例进行了调查。结果如下:
公交车bus,自行车Bicycle,步行walk,其他other
19%, 14%, 56%, 11%

本时期各出行方式转移概率如下表(%)
| |bus | bicycle | walk |other|
|—|—|—|—|—|
|bus | 90 | 4 | 2 |4 |
|bicycle | 7 | 86 | 1 |6 |
|walk | 8 | 7 |80 |5 |
|other | 10 | 2 | 3 |85 |

假设该城市居民每天出行总人数为468万人次，出行人数不变，各出行方式的转移概率也不变

问题：

  1. 预测2006年该城市乘公交出行的人数
  2. 经历足够长的时间，求出行方式的比例是多少？

解：

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
## 分析：这是一个时间齐次马尔科夫过程，可根据转移矩阵的初始定义进行推断
## 第一问
## 根据题目写出转移矩阵
T <- matrix(c(90, 4 , 2 , 4 ,
              7 , 86, 1 , 6 ,
              8 , 7 , 80, 5 ,
              10, 2 , 3 , 85 )/100,
              nrow = 4, ncol = 4, byrow = TRUE)
# 初始矩阵
p <- matrix(c(19, 14, 56, 11)/100, nrow = 1, ncol = 4, byrow = TRUE)
## 下一年的概率应该为当年分配概率和转移矩阵的乘积
## 2003
p1 <- p%*%T
## 2004
p2 <- p1%*%T
## 2005
p3 <- p2%*%T
## 2006
p4 <- p3%*%T
## 2006年乘坐公交车出行的总人数应为
res <- 468 * p4[1]

## 第二问，用计算机模拟的方法，通过对转移矩阵的平衡状态近似求解
## 初始化一个空向量
s <- c()
## 假设一个人在初始时刻选择1公交车出行
s[1] <- 1
## 则其在t1时刻选择任何一种出行方式的概率如下
T[s[j-1],]
## 但是他选择的出行方式可以是随机的，故用sample按照前一个状态的概率，随机抽取一次
 res <- sample(1:4, size = 1, prob = T[s[j-1],])
## 抽取的结果，就是t1时的状态
## 而 t2时的状态只受到t1时状态的影响，因此又回到T[s[j-1],]，至此完成一次模拟
## 每一次抽样都是只受到前一次抽样的影响
for (j in 2:50000)
   s[j] <- sample(1:4, size = 1, prob = T[s[j-1],])
## 在进行多次模拟后，马尔可夫链逐渐收敛。
### 模拟50000代的概率分配如下
res <- table(s[1:50000])/50000
names(res) <- c("bus", "bicycle", "walk", "other")

题外

无论假设s[1] = 1,2,3 还是4，进行多次模拟后，所得的结果是非常接近的。这也表明，马尔可夫过程的平衡状态与初始值无关。

参考文献

待添加

二项分布的对数似然面:算法与R语言实现

对数似然面是极大似然估计的结果之一，这里介绍在R语言中如何绘制对数似然面（Log Likelihood Surface）,这里用到的对数似然函数是二项分布的，对于泊松分布，正态分布，负二项分布等，只需对对数似然函数进行适当修改即可。

绘制二项分布的对数似然面举例：
主要过程

  1. 写出对数似然函数
  2. 用nlm求对数似然函数的估计值
  3. 若绘制似然曲线，则固定一个参数值，并变化另外一个参数，计算Loglikelihood，用lines添加曲线。
  4. 若绘制趋势面，则先要计算参数组合矩阵，并求各栅格处的 LogLikelihood，用persp, image, contour等函数即可。也可以用lattice函数包的wireframe()

以下的例子给出全部的R代码和相应说明。

该例子引自:
http://www.unc.edu/courses/2010fall/ecol/563/001/docs/lectures/lecture10.htm#Rcode

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
109
110
111
112
113
114
## 实例： 假设有以下数据, 记录每个枝条上蚜虫的数量，该数据出自
##Krebs, Charles. 1999. Ecological Methodology.
## Menlo Park, CA: Addison-Wesley. 130页

###########################################################
#######################原始数据############################
###########################################################

### 数据导入

num.stems <- c(6,8,9,6,6,2,5,3,1,4)
aphids_count <- c(0,1,2,3,4,5,6,7,8,9)
aphids <- data.frame(aphids_count, num.stems)

### Raw Data,整理出二项分布拟合所需要的格式
aphids_raw <- rep(aphids_count, num.stems)

###########################################################
####################写出对数似然函数#######################
###########################################################

## 二项分布的对数似然函数，用来绘制对数似然曲线和趋势面

nnominb.loglik <- function(p) {
    res <- sum(log(dnbinom(aphids_raw, mu = p[1], size = p[2])))
    return(res)
}
## 二项分布的负对数似然函数，nlm()函数只能求最小值，因此极大

## 似然值时必须用到此函数
neg.nnominb.loglik <- function(p) {
    res <- sum(log(dnbinom(aphids_raw, mu = p[1], size = p[2])))
    return(-res)

}

###########################################################
###################估计对数似然函数极值####################
###########################################################

#利用nlm估计极值
out.NB <- nlm(neg.nnominb.loglik, c(4,4))
# estimate 是对两个参数的估计，即对neg.nnominb.loglik函数参数p的估计
out.NB

###########################################################
###################图 1 ##绘制似然曲线#####################
###########################################################

 

### x 序列
x <- seq(0.1, 10, by = 0.1)

### 计算y
y_max <- sapply(x, function(x) nnominb.loglik(p = c(out.NB$estimate[1], x)))
plot(y_max ~ x, type = "l", ylab = "Log Likelihood", xlab = "Size")

### 当mu = 2时
y_2 <- sapply(x, function(x) nnominb.loglik(p = c(2, x)))
lines(y_2 ~ x, col = 2)

### 当mu = 5时
y_5 <- sapply(x, function(x) nnominb.loglik(p = c(5, x)))
lines(y_5 ~ x, col = 3)

### 添加图例
legend('bottomright', c(expression(mu=='MLE (3.46)'),
        expression(mu==2), expression(mu==5)),
        col=c(1,2,3), lty=c(1,1,1), cex=.9, bty='n')

###########################################################
##################图 2#绘制栅格图##########################


### 绘制栅格图

## 按照值设色
image(seq(2,5,.1), seq(1,4,.1), z.matrix, col = cm.colors(30),
      xlab = expression(mu), ylab = "Number of Trials")
## 添加等高线
contour(seq(2,5,.1), seq(1,4,.1), z.matrix, add = TRUE)
## 极大似然值
points(x = out.NB$estimate[1], y = out.NB$estimate[2], pch = 3)

###########################################################
###############图 3绘制3D 对数似然面#######################
###########################################################


### 计算趋势面原始数据
### 构建网格
grids <- expand.grid(mu=seq(2,5,.1), theta=seq(1,4,.1))
### 计算z值
grids$z <- apply(grids, 1, function(x) nnominb.loglik(p = c(x[1], x[2])))
### 趋势面原始数据
z.matrix <- matrix(grids$z, nrow=length(seq(2,5,.1)))


### 绘制三维趋势面

persp(seq(2,5,.1), seq(1,4,.1), z.matrix,
      xlab=expression(mu), ylab=expression(theta), zlab='Log Likelihood',
      ticktype='detailed', theta=30, phi=20, d = 4)
###########################################################

##########图 4##绘制对数似然面的三维设色图#################
library(lattice)       
wireframe(z~mu*theta, data=grids, xlab=expression(mu),
          ylab=expression(theta), zlab="Log-nlikelihood",
          scales = list(arrows = FALSE), drape=TRUE,
          par.settings=list(par.zlab.text=list(cex=.75),
          axis.text=list(cex=.75)))

用R语言绘制似然面Likelihood Surface

本文的目的，是直观表达极大似然估计。给出一组数，若用正态分布拟合，则其平均值和标准差最可能的取值是什么？

极大似然面的顶点，就是最优估计。实际情况下，是对对数似然函数Log Likelihood取极值。这里只按照似然（Likelihood）的原始定义，直接求Likelihood 。

likelihood即为各种概率密度的乘积

注意R code 中的 prod(dnorm(dat, mean_val[i], sd_val[j]))

以下是绘制拟合正态分布的似然面。

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
### 绘制似然面： Likelihood Surface

## Author：jinlongzhang01@gmail.com

LikSurf <- function(dat, precision = 0.05){

    mean_val <- seq(min(dat), sqrt(max(dat)), by = precision)

    sd_val <- seq(0, exp(sd(dat)), by = precision)

    mat <- rep(NA, length(mean_val)*length(sd_val))

    dim(mat) <- c(length(mean_val), length(sd_val))

    for(i in 1:length(mean_val)){

       for(j in 1:length(sd_val)){

           mat[i, j] <- prod(dnorm(dat, mean_val[i], sd_val[j]))

       }

    }

    #contour(mat)

    res <- list(xaxis = axisTicks(range(mean_val), log = FALSE),

                yaxis = axisTicks(range(sd_val), log = FALSE),

                 mean_val = mean_val, sd_val = sd_val,

                 mat = mat, pre = precision)

    class(res) <- "LikSurf"

    return(invisible(res))

}

 

print.LikSurf <- function(x){

    cat("This is a likelihood surface with", "nx_range: ",

    range(x$xaxis), "ny_range:", range(x$yaxis),

    "nprecision", x$pre, "n")

}

 

contour.LikSurf <- function(x, ...){

    contour(x = x$mean_val, y = x$sd_val, z = x$mat, axes = FALSE, xlab = "Mean", ylab = "sd", ...)

    axis(1, x$xaxis); axis(2, x$yaxis); box()

}

 

image.LikSurf <- function(x, ...){

    image(x = x$mean_val, y = x$sd_val, z = x$mat, axes = FALSE, xlab = "Mean", ylab = "sd", ...)

    axis(1, x$xaxis); axis(2, x$yaxis); box()

}

 

### 以上函数拷贝到R中即可

### 应用举例

dat = c(-0.77, 6.50, 2.17, -0.51, 1.1, 2.52, 0.17, -0.01, 1.02, 0.11, -0.31, -0.12, 0.25, 0.84)

#dat = rnorm(10)

lks <- LikSurf(dat, precision = 0.05)

image(lks, col = cm.colors(20))

contour(lks, drawlabels = TRUE, add = TRUE)

似然面： 求极大似然

6亿年来的海陆变迁

img

晚元古代（600MA）6亿年前

元古代末期（560Ma） 5.6亿年前

早寒武纪（540Ma） 5.4亿年前

晚寒武纪（500MA） 5.0亿年前

中奥陶纪（470Ma） 4.7亿年前

晚奥陶纪（450MA） 4.5亿年前

志留纪（430Ma） 4.3亿年前

早泥盆纪（400MA） 4.0亿年前

晚泥盆纪（370Ma） 3.7亿年前

密西西比世（340Ma） 3.4亿年前

宾夕法尼亚世（300MA） 3.0亿年前

早二叠纪（280MA） 2.8亿年前

晚二叠纪（260Ma） 2.6亿年前

早，中三叠纪（240MA） 2.4亿年前

晚三叠纪（220MA） 2.2亿年前

早侏罗纪（200MA） 2.0亿年前

中侏罗纪（170Ma） 1.7亿年前

晚侏罗纪（150MA） 1.5亿年前

早白垩纪（120MA） 1.2亿年前

晚白垩纪早期（105Ma） 1.05亿年前

晚白垩世（90Ma）9000 万年前

K - T界限（65Ma） 6500万年前

始新世（50MA） 5000万年前

渐新世（35MA） 3500万年前

中新世（20MA） 2000万年前

img

更新世（50KA） 5万年前

现在

资料来源： http://jan.ucc.nau.edu/~rcb7/mollglobe.html

如何编写R函数

R语言实际上是函数的集合，用户可以使用base，stats等包中的基本函数，也可以自己编写函数完成一定的功能。但是初学者往往认为编写R函数十分困难，或者难以理解。这里对如何编写R函数进行简要的介绍。

函数是对一些程序语句的封装。换句话说，编写函数，可以减少人们对重复代码书写，从而让R脚本程序更为简洁，高效。同时也增加了可读性。一个函数往往完成一项特定的功能。例如，求标准差sd,求平均值，求生物多样性指数等。R数据分析，就是依靠调用各种函数来完成的。但是编写函数也不是轻而易举就能完成的，需要首先经过大量的编程训练。特别是对R中数据的类型，逻辑判别、下标、循环等内容有一定了解之后，才好开始编写函数。 对于初学者来说，最好的方法就是研究现有的R函数。因为R程序包都是开源的，所有代码可见。研究现有的R函数能够使编程水平迅速提高。

R函数无需首先声明变量的类型，大部分情况下不需要进行初始化。一个完整的R函数，需要包括函数名称，函数声明，函数参数以及函数体几部分。

  1. 函数名称，即要编写的函数名称，这一名称就作为将来调用R函数的依据。

  2. 函数声明，包括 <- function, 即声明该对象的类型为函数。

  3. 函数参数，这里是输入的数据，函数参数是一个虚拟出来的一个对象。函数参数所等于的数据，就是在函数体内部将要处理的值，或者对应的数据类型。 函数体内部的程序语句进行数据处理，就是对参数的值进行处理 ，这种处理只在调用函数的时候才会发生。函数的参数可以有多种类型。R help的界面对每个函数，及其参数的意义及所需的数据类型都进行了说明。

  4. 函数体

常常包括三部分.

  • （1）. 异常处理

输入的数据不能满足函数计算的要求，或者类型不符， 这时候一定要设计相应的机制告诉用户，输入的数据在什么地方有错误。 错误又分为两种。

第一种， 如果输入的数据错误不是很严重，可以经过转换，变为符合处理要求的数据时， 此时只需要给用户一个提醒，告知数据类型不符，但是函数本身已经 进行了相应的转换。

第二种，数据完全不符合要求，这种情况下，就 要终止函数的运行，而告知因为什么，函数不能运行。这样，用户在 使用函数的情况先才不至于茫然。

  • （2）. 运算过程

包括具体的运算步骤。 运算过程和该函数要完成的功能有关。

R运算过程中，应该尽量减少循环的使用，特别是嵌套循环。R提供了 apply，replicate等一系列函数，来代替循环，应该尽量应用这些函数， 提高效率。 如果在R中实在太慢，那么核心部分只能依靠C或者Fortran 等语言编写，然后再用R调用这些编译好的模块，达到更高的效率。

运算过程中，需要大量用到if等条件作为判别的标准。if和while都是需要数据TRUE/FALSE这样的逻辑类型变量，这就意味着，if内部，往往是对条件的判别，例如 is.na, is.matrix, is.numeric等等，或者对大小的比较，如，if(x > 0)， if(x == 1)， if(length(x)== 3)等等。if后面，如果是1行，则花括号可以省略，否则就必须要将所有的语句都放在花括号中。这和循环是一致的。

例子1

1
2
3
4
5
6
7
8
9
10
11
12
## if与条件判断

fun.test <- function(a, b, method = "add"){
    if(method == "add") { ## 如果if或者for/while；
        res <- a + b       ## 等后面的语句只有一行，则无需使用花括号。
     }
    if(method == "subtract"){
        res <- a - b
    }
    return(res)           ## 返回值

}
1
2
3
### 检验结果
fun.test(a = 10, b = 8, method = "add")
fun.test(a = 10, b = 8, method = "substract")

for循环有些时候是必须要用到的，for循环内部，往往需要用下标，访问数据内的一定元素，例如向量内的元素，这时候用方括号表示。一维的数据组合，或者数组，常常称为向量。二维的数据组合，往往称为矩阵，或者数据框。具体的访问方式主要是方括号内部有没有逗号的区别。for循环或者while循环有时候让人觉得比较困惑，可能需要专门的时间进行讲解。

例2

1
2
3
4
5
6
7
8
9
10
11
### for循环与算法

test.sum <- function(x)
{
    res <- 0               ###  设置初始值，在第一次循环的时候使用
    for(i in 1:length(x)){
        res <- res + x[i] ## 这部分是算法的核心，
##总是总右面开始计算，结果存到左边的对象
    }
    return(res)
}
1
2
3
4
### 检验函数
a <- c(1,2,1,6,1,8,9,8)
test.sum(a)
sum(a)

无论是什么样的函数，算法才是最关键的。往往需要巧妙得设计算法，让函数快捷高效。

（3）. 返回值。

返回值就是函数给出的结果。打个比方，编写一个函数，就像自己攒一个机器，例如现在攒好 一台豆浆机，该豆浆机要求输入大豆，输入的大豆就是参数， 返回的结果，就是豆浆。如果该豆浆机需要不停地输入大豆， 而不能产出豆浆，这样的机器就一定会被扔掉。函数也是一样的， 需要给出返回值。 R中默认的情况是将最后一句作为返回值。但是为了函数的可读性起见，应该尽量指名返回值。返回值用return()函数给出。 函数在内部处理过程中，一旦遇到return()，就会终止运行， 将return()内的数据作为函数处理的结果给出。

下面举例说明R函数的编写方法。

例3 计算标准差

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
sd2 <- function(x)
{
   # 异常处理，当输入的数据不是数值类型时报错   
   if(!is.numeric(x)){
      stop("the input data must be numeric!n") 
   }

   # 异常处理，当仅输入一个数据的时候，告知不能计算标准差

   if(length(x) == 1){
      stop("can not compute sd for one number,
           a numeric vector required.n")
   }
   ## 初始化一个临时向量，保存循环的结果，
   ## 求每个值与平均值的平方  
   x2 <- c()
   ## 求该向量的平均值 
   meanx <- mean(x)

   ## 循环 
   for(i in 1:length(x)){   
       xn <- x[i] - meanx  
       x2[i] <- xn^2
   }
      ## 求总平方和
   sum2 <- sum(x2)
   # 计算标准差
   sd <- sqrt(sum2/(length(x)-1))
   # 返回值
   return(sd)
}

 

## 程序的检验
## 正常的情况
sd2(c(2,6,4,9,12))

## 一个数值的情况
sd2(3)

## 输入数据不为数值类型时
sd2(c("1", "2"))

这样，一个完整的函数就编写完成了。当然，实际情况下，函数往往更为复杂，可能要上百行。但是好的编程人员往往将复杂的函数编写成小的函数。以便于程序的修改和维护，即使其中出现错误，也很好修改。

再有就是编写R函数时一定要注意缩进，编辑器用Notepad++, TinnR, Rstudio等，同时用等距字体（如Consolas, Courier new等）和语法高亮显示。这样便于快速寻找到其中的错误。

感谢 黄继红博士，杜彦君博士，毛岭峰博士，饶米德，冯刚对本文提出的意见和建议。

2011年8月19日 于中科院植物所

2011年8月22日修改

Inno setup 制作安装文件以编译FigTree为例

Windows下安装文件的制作为软件的安装和卸载提供了方便。用户不再用担心可执行文件及动态连接库等保存在固定的文件夹中，只需要进行软件的安装和卸载，就可以将全部需要的可执行文件，以及程序运行所需要的文件创建和删除。这为程序管理提供了很大的方便。安装文件的制作有很多软件，如setup factory, Install shield wizard等等，但是对于小型软件，

Inno setup，这样的小型安装文件制作工具就已经足够了。

Inno setup开始于1997年，是完全免费的，从稳定上和兼容性上，甚至超过了一些商业软件。

这里介绍一下如何用Inno setup创建一个setup文件。以FigTree为例。

FigTree是爱丁堡大学的Andrew Rambaut编写的绘制进化树的软件，用Java写成，可以在多种平台上运行。但是Andrew提供的Windows下的程序包实际上是一个包含exe文件和动态连接库的文件夹，使用起来不是特别方便，因此用Inno

Setup建立一个Windows安装文件，便于对程序的管理和使用。

FigTree可以在 http://tree.bio.ed.ac.uk/software/figtree/ 下载。下载Windows下的zip程序包，解压缩。

下载Inno setup http://www.jrsoftware.org/isinfo.php 并安装。

创建一个新的Project，按照界面给出的提示，Inno setup会自动生成一个脚本。编译该脚本，即可生成所需的exe文件。一般来说，首先要指定exe文件。其次要给出程序所要包含的文件夹，该文件夹中包含运行该exe所需的动态连接库等。

还有就是创建ico文件，ico文件是Windows的图标，如果现有jpg或者png等其他格式的文件。可以在 http://iconverticons.com/转换成ico文件，作为该setup显示的图形。

了解以上信息，相信读者都能够编译出自己的Windows setup文件了。

FigTree v1.3.1.zip

如何撰写生态学科研论文？

科研论文是科学家贡献的最直接标志。好的研究论文往往具有很强的生命力，在十几年，甚至几十年后还会有人引用。研究论文，作为当代最重要的科学研究著作的形式，理应受到特别的关注。而对于广大研究生而言，学会正确的撰写科研论文，也是学习的重要内容。科研论文最重要的目的是简明扼要得阐述科研成果。经过同行评议的科研论文意味着该成果已经正事进入学术圈，是对科学的重要贡献。因此，写好科研论文是十分必要的。

对于理论性很强的生态学科而言，更是如此。为此，笔者做了如下总结，希望和生态学同仁共勉。

数据分析好后，一般已经获得了撰写论文结果所需的图表等，因此，撰写研究论文一般按照以下顺序

  1. 结果
  2. 方法
  3. 讨论
  4. 前言
  5. 摘要

下面按照论文的一般顺序介绍各部分的主要写法

前言

这部分主要是为了引出自己的研究的问题，以及这些问题为什么重要。

  • 第一段 研究的背景

  • 第二段 研究领域的question,以及相应的假说;

  • 第三段， 关于这些假说，人们都进行了哪些研究，但是还有什么地方不清楚

  • 第四段， 提出自己要研究的问题，或者依据现有的理论，预测出的结果，经过分析之后，会给出怎样的结果，即要检验的hypothesis

  • 第五段， 简要介绍这篇文章用什么数据，用什么方法以引出方法这部分

方法

这一部分介绍研究的材料和方法，以便别人能够重复你的研究。

一般先写数据的获取，包括样地的描述，调查的方法，然后是统计分析的方法，最后介绍统计分析的软件。

结果

这一部分给出分析最重要的结果。尽量用图给出，实在不行再用表格给出。给出最主要的信息就可以，因为结果可以有很多。但是支持或者不支持假说的，只有很少的关键的一部分，重点突出这部分

讨论

这一部分是深入阐述研究结果的意义与价值，升华论文的深度。指出研究的不足与进一步的发展方向。

讨论往往需要先回顾一下讨论的问题。读者读到这里，往往忘记了作者要做什么问题，因此再重申一下检验的问题

然后分别对检验的结果，一一进行评述，包括这些结果与前人研究的异同，以及这些结果带来的启示，目的就是深化论文发现结果的意义，将论文的思想表达清楚。

之后要说一下本论文的不足之处，以及可以改进的地方，或者仍然研究不清楚的地方。甚至是研究展望。

最后，对本文的发现和意义进行总结。

摘要

摘要部分往往最后撰写，摘要往往分成两种形式，第一种，直接报道科研数据，以及结果，在摘要中往往能够发现大量的数字；第二种，文字性的报道。但是无论哪种，都要说清楚论文的目的，主要方法，结果和新发现的意义。一篇好的摘要，如同一篇好的论文，虽然不能面面俱到，也是能够吸引读者深入读下去的。

其它

​再之后，就是准备supporting information,包括论文分析数据用到的程序，一些不便放在正文中的图表等。

​经过多次认真细致的修改之后，将格式按照投稿要求改好，才能投稿。 这里祝愿大家马到成功！

PCNM动画展示：生态学与空间统计学

PCNM (Principal Coordinates of Neighbourhood Matrix)是Legendre教授等开发出的表示空间关系的方法，类似于pcoa。常用于生物多样性空间格局研究中，例如在古田山beta多样性研究中，Legendre**(2009)**就使用了这种方法。用样地的规则栅格距离可以对点两两之间的距离关系进行PCNM模拟。 这里用到了vegan程序包将距离矩阵转换成PCNM的特征向量，用imagevect函数绘制PCNM特征向量在样地中的分布。希望对有兴趣的研究人员有帮助。

将以下代码在R中运行即可。

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
#### Function imagevect
#### By Jinlong Zhang  <Jinlongzhang01@gmail.com>
#### Institute of Botany, the Chinese Academy of Sciences, Beijing ,China
#### Aug. 18, 2011
 
imagevect <- function (x, labels, contour = FALSE, gridsize = 20,
                        axes = TRUE, nlabx = 5, nlaby = 5, ...)
{
    require(fields)
   dimension <- function(x, unique = FALSE, sort = FALSE){
        ncharX <- substring(x, 2, regexpr("Y", x)-1)
        ncharY <- substring(x, nchar(ncharX)+3, nchar(x))
        if(unique){
            ncharX = unique(ncharX)
            ncharY = unique(ncharY)
        }
        if(sort){
            ncharX = sort(as.numeric(ncharX))
            ncharY = sort(as.numeric(ncharY))
        }
        res <- list(ncharX, ncharY)
        return(res)
    }
    
    formatXY <- function(x){
        ncharX <- substring(x, 2, regexpr("Y", x)-1)
        ncharY <- substring(x, nchar(ncharX)+3, nchar(x))
        resX <- c()
        for(i in 1:length(ncharX)){
            n0x <- paste( rep(rep(0, length(ncharX[i])),
                         times = (max(nchar(ncharX)) - nchar(ncharX[i]) + 1)),
                         collapse = "", sep = "")
            resX[i] <- paste("X", substring(n0x, 2, nchar(n0x)), ncharX[i], collapse = "", sep = "")
        }
        resY <- c()
        for(i in 1:length(ncharY)){
            n0x <- paste( rep(rep(0, length(ncharY[i])),
                         times = (max(nchar(ncharY)) - nchar(ncharY[i]) + 1)),
                         collapse = "", sep = "")
            resY[i] <- paste("Y", substring(n0x, 2, nchar(n0x)), ncharY[i], collapse = "", sep = "")
        }
        res <- paste(resX, resY, sep = "")
        return(res)
    }
 
    sort.x <- x[order(formatXY(labels))]
    rrr <- dimension(labels, unique = TRUE, sort = TRUE)
    dims <- c(length(rrr[[2]]), length(rrr[[1]]))
    dim(sort.x) <- dims
    par(xaxs = "i", yaxs = "i")
    image.plot(nnn <- t(sort.x), axes = FALSE, ...)
 
    if (contour) {
        contour(nnn, add = TRUE, ...)
    }
   
    if (axes) {
        points(0, 0, pch = " ", cex = 3)                      ## invoke the large plot
       
     get.axis.ticks <-
              function(nlabs = NULL, gridsize = NULL, limit_max = NULL){
    
         ngrid <- (limit_max-0)/gridsize                      ## Obtain number of labels to plot
         per_grid <- 1/(ngrid-1)                              ## Obtain length for each grid
         start <- 0 - (1/(ngrid-1))/2                         ## starting point for the ticks
         stop <- 1 + (1/(ngrid-1))/2                          ## stopping point for the ticks
        
         lab <-(0:nlabs*(limit_max/nlabs))
         at <- seq(from = start, to = stop,
                   by = ((1 + per_grid))/((length(lab)-1)))   ## Position of the ticks
        
         return(list(lab, at))
        
     }
       xaxis.position <- get.axis.ticks(nlabs = nlabx, gridsize = gridsize,
                                limit_max = gridsize * nrow(nnn))
       yaxis.position <- get.axis.ticks(nlabs = nlaby, gridsize = gridsize,
                                limit_max = gridsize * ncol(nnn))
       
        axis(1, labels = xaxis.position[[1]], at = xaxis.position[[2]])
        axis(2, labels = yaxis.position[[1]], at = yaxis.position[[2]])
    }
    invisible(nnn)
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
library(vegan)
 
## 模拟一套数据
X <- seq(10, 590, by = 20)
Y <- seq(10, 390, by = 20)
 
## Label
XY <- expand.grid(X, Y)
names <- paste("X", (XY[,1] + 10)/20 ,"Y", (XY[,2] + 10)/20 , sep = "")
rownames(XY) <- names
## 计算样方两两之间的距离
distXY <- dist(XY)
 
## 进行pcnm
gtsPCNM <- pcnm(distXY)
head(gtsPCNM$vectors)
 
## 绘制pcnm图
for(i in 1:10){
    imagevect(gtsPCNM$vectors[,i], labels = names,
               col = topo.colors(100))
    Sys.sleep(0.3)
}

img

img

img

imgimg

img

img

img

img

img

Tortoise SVN管理本地代码

什么是Tortoise SVN?

TortoiseSVN是一个windows下的文档版本管理的开源软件。用户每次对自己编写的代码进行修改，都会记录在SVN的数据库中。Tortoise SVN能够在设定好的文件夹上添加相应的“对号”，“问号”等标识，标识当前代码的编辑状态，特别是有没有在数据库中保存。 对于代码的修改，用户可以添加相应的标注。对于每一次修改，数据库都有详细的记录，从而保证所编写的文档可以回到作者保存过的任何一个版本。

这种版本控制策略在软件开发中是极为重要的，当然，在R程序包的开发中也十分重要。

用Tortoise SVN管理本地R代码的大体过程如下：

  1. 下载和安装Tortoise SVN软件，各项均选择默认即可，网址如下 http://tortoisesvn.net/downloads.html
  2. 在本地硬盘上创建一个新目录，作为数据库的保存文件夹。例如 D:/packages/phylotools
  3. 右键点击phylotools文件夹，Tortoise SVN>Create a repository here.完成后，打开phylotools文件夹，我们会发现其中新增了一些文件和文件夹。这是版本数据库相应的文件，我们暂且不管。
  4. 在本地硬盘上创建一个新文件夹，例如在C:/developing/. 点击鼠标右键，选择SVN Checkout。我们看到，developing文件夹下出现了一个phylotools空文件夹，该文件夹上有一个绿色的对号。我们发现，该文件夹是空的，绿色的对号表示，文件夹下的内容已经与数据库中的版本相同了。
  5. 在该文件夹下创建新文件，或者将之前编写好的代码拷贝到C:/developing/phylotools文件夹下。此时发现每个文件上都被加上了蓝色的问号，这表明这些文件还没有和数据库链接起来。此时我们回到上级目录，C:packages, 右键点击phylotools文件夹，点击SVN Commit这样，该文件夹下的文件就全部导入数据，并且关联起来了。
  6. 之后对其中任何代码的修改，均可以提供Comments，并且隔一段时间进行保存。
    这样以后恢复到以前的版本，就容易多了。而不用隔一段时间备份一下新文件。因为SVN已经帮你把修改信息全部存到数据库里了。

报告《R与进化生态学：系统发育比较分析概述》

随着测序和进化树推断理论和技术的不断成熟，生态学家正努力将物种系统发育关系整合到生态学研究中，并逐渐发展出一系列研究物种性状和适应性演化规律的分析方法。这些研究方法大量运用计算机模拟，包括随机生灭过程和随机化等统计学手段。目前，这些方法在R中迅速发展。本报告简要介绍了系统发育比较分析中的主要内容，包括：祖先状态重建（Ancestral reconstruction）、物种分化速率随时间的变化(Lineages through time plot)、群落系统发育与物种共存(Community phylogenetics)、生态位进化(Niche evolution),以及这些方法在R的实现。

这个报告是在中国第四届R语言会议上的报告，旨在介绍R中的系统发育比较方法。

下载本幻灯片 pdf

R与进化生态学_张金龙.pdf