博士论文致谢

时光荏苒,三年不经意间就过去了。这三年中有无数憧憬和企盼,也有数不清的酸甜苦辣。百味之中,才足见这三年的厚重。

首先,要感谢恩师马克平研究员,能够在马老师的指导下从事研究工作,是我的荣幸。马老师渊博的学识,谦逊的品格,严谨的治学态度,深厚的文史修养,都是我一生的学习榜样。

感谢中国科学院生态环境中心陈圣宾博士提供部分保护区名录数据,以及在数据库整理中的繁重工作,并提出建设性的意见和建议。感谢中国科学院华南植物园葛学军研究员为DNA条形码实验的顺利进行给予诸多的便利。华南植物园裴男才博士、颜海飞博士,浙江古田山自然保护区陈声文工程师,中国科学院植物研究所冯刚、浙江师范大学饶米德参与了古田山DNA条形码样品的采集工作。特别是裴男才博士完成了多半部分DNA样品的采集整理DNA提取等,并提供了构建进化树的详细指导,才能使实验及后续分析得以顺利完成。

感谢米湘成副研究员、美国密歇根州大大学的Nathan Swenson博士对论文的部分内容提出修改意见。感谢黄继红博士、阳文静女士、李宗善博士、刘勇波博士、陈国科博士参与部分章节的讨论,并提出修改意见。感谢北京林业大学博士生范娟对论文提出的宝贵建议。

感谢杨永老师、梁宇老师、陈铁梅老师、杜晓军老师、裴克全老师、张守仁老师、赖江山老师、付增娟老师,以及厦门大学李振基教授,中南大学生命学院王丽副教授,湖南商学院李陈华教授,植物所研究生处赵剑锋老师、吴海燕老师和张文娟老师给予的帮助。树木年轮组的张齐兵研究员、邱红岩老师、刘彩云老师,在很多方面给予我亲切的关怀和热情的帮助,在此表示衷心的感谢。朱丽博士、陈彬博士、金冬梅博士、邵长亮博士、丁琼博士、孙秀峰博士、于胜祥博士、黄菊莹博士、余华博士、王建秀博士、毛建丰博士、祝燕博士、吕立新、杜彦君、张潮,以及戴文燕、张传岭、毛岭峰、刘晓娟、刘冰、冯刚、宋凯、韩娟娟、张朝、张灿娟、曼兴兴、程佳佳、郭翠燕,以及华中农业大学武欣同学都给予我很多帮助,在此一并表示感谢。

最后感谢我的父母、妹妹张玉翠和姑父、姑母对我的理解和关爱。

张金龙

于北京香山中国科学院植物研究所地主大院

2011年6月19日

两棵进化树一致性的校对

首先, 要保证获得的进化树是有根树, 设定好外类群,并且,分支出现的次序经过排序reorder,这可以在Figtree软件中, 对newick格式的进化树进行操作, export as currently displayed 即可。

其次, 用R软件ape程序包的 cophyloplot函数, 可以实现, 将两棵进化树面对面绘制, 并将名称相同的分类单元用线段联系起来, 此时, 如果两个进化树结构有差别, 就一目了然了。

举例:

1
2
3
4
library(ape)
labs1 <- paste("t", 1:10, sep = "")
labs2 <- paste("t", 1:10, sep = "")
cophyloplot(x = rtree(20), y = rtree(10), space = 100, length.line = 0.1, gap = 10,assoc = as.matrix(data.frame(labs1, labs2)))

R绘图一则 (2017年1月9日修订)

2017年1月9日修订

这里给出R常见绘图用到的参数。供生成出版质量的TIFF图。

调整R绘图的字体,一页多图, 分辨率600dpi, tiff, 压缩方式lzw

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
tiff(filename = "test.tif",
width = 3600, height = 3600, units = "px", pointsize = 8,
compression = "lzw", res = 600) ## 打开绘图器
par(mfrow = c(2,2), font.lab = 6, font.main = 6,font.axis = 6, font = 6) ## 分割成2行2列,共四个图

## 图1
## 标题部分斜体
hist(rnorm(200), main = expression(paste("Testing ", italic(Hist1))), col = "gray")
# 添加边框
box()
## 图中添加图例
legend.label <- c("Type I", "Type II", "Type III")
legend(1.5, 30, legend = legend.label, pch = c(19, 20, 21), col = c(1, 2, 3))
legend(-3.0, 30, legend = legend.label, lty = c(1, 2, 3), col = c(1, 2, 3))

## 图2
## 标题中显示度等特殊字符,标题中显示希腊字母
hist(rnorm(200), main = expression(paste(plain(sin) * phi, " and ", plain(cos) * phi)), xlab = expression(paste("Latitude ", degree)), col = "gray")
## 图中添加公式
text(-2, 30, expression(bar(x) == sum(frac(x[i], n), i==1, n)), cex = 1.2)

## 图3
x <- rnorm(100)
m=mean(x)
sd=sd(x)
hist(x)
x00 <- seq(-10, 10, 0.01)
lines(x00, m*dnorm(x00,m,sd), lty=1, col="red")

## 图4
y = 1:50 + (5*rnorm(50) + 20)
x = 1:50
plot(y ~ x, main = "Scatter plot", col = "gray")
abline(lm(y ~ x), col = 2)
dev.off() ## 关闭绘图器
getwd() ## 查看图保存的位置

LyX生成中文pdf文档

摘要 本文概述了用LyX编译中文文档的过程, 可以作为LyX中文入门的指南。

引言

LyX软件可以看做是LaTeX的前端。LaTeX则是学术领域应用很广的排版软件。计算机领域有大量的专业人士以及爱好者醉心于LaTeX。LaTeX不仅可以排版论文、图书、还可以用来写简历,制作做幻灯片,生成的文件简洁、 清晰,十分美观。同时,LaTeX最优越之处在于能够十分方便得输入数学公式。

但是用LaTeX编写文档,用户常常需要面对大量的代码。此时不得不记忆大量的命令。有人形容,用LaTeX写文档, 就相当于写程序。幸运的是,有一批软件实现了对LaTeX代码进行图形化的管理。LyX就是其中的佼佼者。LyX软件实现了LaTeX代码管理的图形化,对于入门级用户来说, 即便没有学过LaTeX,也能快速做出漂亮的LaTeX文档。

1 下载安装程序

安装好后,如图1[图-1-LyX2.0.4版本界面(2012)]所示
img

图 1 LyX2.0.4版本界面(2012)

2 设定文档的类型

每一次新建文档,首先要选择文档的类型,这里以article的类型举例。

点击 File > New ,新建一个空白文档, 点击 Document> Settings

分别点击左侧的导航栏,在Document class中选择 CTeX,这是由于我们的文档中包含中文,参考文献, 目录的格式,CTeX类型的文档能给予很好的支持。

Language中选 Chinese (simplefied),Encoding 选 Unicode (XeTeX) (utf8) 这是由于我们在后续编译LyX文档的时候, 要用到XeTeX程序, 而要在正文中使用中文,xeCJK宏包的效果很好。

在output设置Default output format 为 PDF (XeTeX),我们用XeTeX编译LyX文档, 生成PDF文档。

在LaTeX Preamble中输入如下LaTeX 代码

usepackage{xeCJK}

setCJKmainfont{SimSun}

img

图 2 增加两行命令,作为LaTeX开始时导入的命令

LaTeX Preamble保证在LyX生成LaTeX源代码的时候,将这部分加入到LaTeX源代码中。usepackage{xeCJK} 是导入xeCJK宏包。由于xeCJK包只有设定了字体之后才能生效,因此需要输入setCJKmainfont{SimSun}。 如果不输入setCJKmainfont{SimSun}, 将出现PDF能够编译, 但是中文不能显示出来的问题。

这样, article的中文环境格式就设定好了。

3 调整文本格式

默认条件下,文本部分都是作为Standard的格式,这也是正文中最多的格式。但是一篇论文,需要有标题(Title) 作者, 日期,摘要,以及前言,材料与方法,结果和讨论等,有时候还要有插图和参考文献的引用。这些内容,都需要逐项设定。设定的方法,就是在左侧的下拉菜单中,选择对应的格式。例如LyX生成pdf中文文档,应设置为Title ,作者Author。

前言, 材料与方法, 结果, 讨论,都通过类似的方法设定。

3.1 公式的输入

有时候, 在材料和方法中要输入公式

例如

$$NRI=({MPD}{Observed}-mean({MPD}{Null}))/(sd({MPD}_{Null}))$$

则点击Insert菜单,选择要输入的公式类型, 则屏幕下方将出现公式管理的按钮,点击即可调整。

Insert菜单,不仅包括包括公式,还包括输入标签,以及插入参考文献,插入目录,图表等功能。公式、图表、 参考文献等物件, 在LyX中会自动编号,动态调整,非常方便。

3.2 参考文献的输入

设定参考文献的方法, 在文章的最后, 输入要引用的参考文献的纯文本。

例如 S.W. Kembel, P.D. Cowan, M.R. Helmus, W.K. Cornwell, H. Morlon, D.D. Ackerly, S.P. Blomberg, and C.O. Webb. 2010. Picante: R tools for integrating phylogenies and ecology. Bioinformatics 26:1463-1464.

选中该部分文字,在设定格式的菜单栏中,选定Biblography,将该部分作为参考文献。当然, 也可以引用用BibTeX,方法是,先在计算机上创建一个BibTeX纯文本文档, 注意以utf8格式编码,然后在正文中引用即可。很多软件或网站,如Mendeley,R,及google scholar都提供BibTeX格式的文献导出。而参考文献的引用,点击Insert>Citation即可。如果使用的是BibTeX,直接输入Insert>List/TOC>BibTeX Biblography,选择扩展名为.bib的文件即可。

.bib文件内的格式为

1
2
3
4
5
6
7
8
@Article{, 
title = {Picante: {R} tools for integrating phylogenies and ecology},
author = {S.W. Kembel and P.D. Cowan and M.R. Helmus and W.K. Cornwell and H. Morlon and D.D. Ackerly and S.P. Blomberg and C.O. Webb},
journal = {Bioinformatics},
volume = {26},
pages = {1463--1464},
year = {2010},
}

假设,撰写论文中有一句话,“NRI[equation 1]表示群落内物种系统发育关系的远近, 通过picante程序包可以计算SES.MPD等价于NRI[equation 1]” 需要引用公式以及计算NRI用到的软件picante

则先创建一个标签, 方法为 Insert>Label, 随后输入Insert cross reference ,在出现的对话框中选择引用该标签即可。

3.3 插入图片

通常情况下,文章中出现图片会更吸引人。点击Insert>Graphics即可输入图片。推荐使用pdf格式的图片。 因为pdf格式的图片在整个文档编译成pdf后, 仍然是矢量的, 不会有分辨率的损失。jpg,png等格式也都是可以的。

4 文档的编译

前面的内容都准备好后,只需点击View>View\[pdf(XeTeX)]即可完成编译,此时, 你的pdf阅览器应该能够自动打开编译好的pdf文档。看,一篇用LyX编译的文档诞生了。其它类型的文档, 还包括图书,幻灯片(Beamer)等,方法大同小异,只不过设定的文档结构和格式不同罢了。例如,报告或者图书, 需要目录,图书在最后往往还要设置索引等内容。

世界时与力学时之差:Delta T的计算

为了便于天文计算,需要使用均匀流逝的时间标尺, 这种时间系统最早称为历书时, 后来,由于天文计算均以力学规律为基础, 又改称力学时。 人们发现原子内电子的能级跃迁振荡频率非常均匀, 如铯133的共振频率为每秒9192631770周。遵循这一原理, 科学家制作出了原子钟,作为均匀的时间标尺。因此, 原子时又成为力学时的等价时间。

通常, 人们用到的时间为世界时(UT)。 世界时是完全根据地球自转而确定的时间参考系统。 但是地球自转是不均匀的(章动), 同时也包括一些未知的因素, 造成了世界时与力学时之间的差异。

星历表、日月食、行星动态等, 都是先通过力学时为基准计算的, 再通过世界时与力学时之差转换为世界时。因此,

在天象计算中, Delta T是决定计算精度的重要因素之一。

但是Delta T只能通过实测来获得。 对远古时代的Delta T, 只能通过古代的可靠的天象记录进行推测, 例如, 2500年前, Delta T 可能达到4到5小时。

为此, 若需要某一特定时刻的Delta T, 则需要使用插值法。 插值法包括线性插值, 二次插值, 三次插值, 拉格朗日插值法等多种, 这里提供三次插值法的R函数, 以及相应的插值表, 以备参考。

详细参考 天文算法 http://www.fjptsz.com/xxjs/xjw/rj/117/03.htm

这里给出R函数, 计算Delta T

数据和函数改编自 许剑伟 老师的 寿星万年历 javascript 源代码。

int2 <- function (v) { floor**(v)**; }

dt_ext <- function**(y, jsd){** #二次曲线外推,用于数值外插

​ dy = **(y - 1820)/**100;

​ return (-20 + jsd ***** dy ***** dy);

}

dt_calc <- function**(y){** #传入年, 返回世界时UT与原子时(力学时 TD)之差, ΔT = TD - UT

​ y0 = dt_at**[length(dt_at)** - 1**]**; #表中最后一年

​ t0 = dt_at**[length(dt_at)]**; #表中最后一年的 deltatT

if**(y >= y0){**

​ jsd = 31; # sjd是y1年之后的加速度估计。

​ # 瑞士星历表jsd=31, NASA网站jsd=32, skmap的jsd=29

if**(y > y0 + 100){**

​ return**(dt_ext(y, jsd))**

};

​ v = dt_ext**(y, jsd)**; #二次曲线外推

​ dv = dt_ext**(y0, jsd)** - t0; # ye年的二次外推与te的差

​ return (v - dv(y0 + 100 - y)/100)*;

}

​ d = dt_at;

for**(i in seq(1, length(d), by = 5)** ) {

if**(y < d[i + 5])** break; # 判断年所在的区间

}

​ t1 = (y - d[i]) / (d[i + 5] - d**[i])** ***** 10 #### 三次插值, 保证精确性

​ t2 = t1 ***** t1

​ t3 = t2 ***** t1;

​ res <- d**[i + 1]** +d[i + 2] ***** t1 +d[i + 3] ***** t2 +d[i + 4] ***** t3

​ return**(** res );

}

# TD - UT1 插值表

dt_at = c**(**

**-**4000, 108371.7, **-**13036.80, 392.000, 0.0000,

**-**500, 17201.0, **-**627.82, 16.170, **-**0.3413,

**-**150, 12200.6, **-**346.41, 5.403, **-**0.1593,

150, 9113.8, **-**328.13, **-**1.647, 0.0377,

500, 5707.5, **-**391.41, 0.915, 0.3145,

900, 2203.4, **-**283.45, 13.034, **-**0.1778,

1300, 490.1, **-**57.35, 2.085, **-**0.0072,

1600, 120.0, **-**9.81, **-**1.532, 0.1403,

1700, 10.2, **-**0.91, 0.510, **-**0.0370,

1800, 13.4, **-**0.72, 0.202, **-**0.0193,

1830, 7.8, **-**1.81, 0.416, **-**0.0247,

1860, 8.3, **-**0.13, **-**0.406, 0.0292,

1880, **-**5.4, 0.32, **-**0.183, 0.0173,

1900, **-**2.3, 2.06, 0.169, **-**0.0135,

1920, 21.2, 1.69, **-**0.304, 0.0167,

1940, 24.2, 1.22, **-**0.064, 0.0031,

1960, 33.2, 0.51, 0.231, **-**0.0109,

1980, 51.0, 1.29, **-**0.026, 0.0032,

2000, 63.87, 0.1, 0, 0,

2005, 64.7, 0.4, 0, 0, #一次项记为x, 则 10x=0.4秒/年*(2015-2005), 解得x=0.4

2015, 69

);

## 输入年, 计算delta T

dt_calc**(1800)**

dt_calc**(1960)**

dt_calc**(1965)**

dt_calc**(2012)**

瑞典:美丽的西博滕

7月22日到28日, 我在瑞典北部的西博滕省西部的Asele进行植物调查。

Asele地区是拉普兰的一部分, 湖泊众多, 景色怡人。拉普兰是传说中圣诞老人居住的地方。 这里一年大部分时间被积雪覆盖, 夏季阳光充足, 但是物种比较单一。 分布着广泛的欧洲赤松(Pinus sylvestris)、欧洲刺柏(Juniperus communis)、欧洲山杨(Populus tremula)、 欧洲白桦(Betula pendula)林 。

植物调查的内容, 是在划分好的边长5km的方格内,尽可能详尽得记录出现的植物,并编写植物名录。这是为了编写西博滕省植物志的需要。西博腾省共包括2700多个方格,此项调查已经进行了三十年,现在即将完成。三十年来,参加调查的,是几百位志愿者。 此次参加调查的志愿者,是来自瑞典各地,有些驱车十几小时,从乌普萨拉甚至是更遥远的地方驱车而来,其中有大学教授, 也有林业工人, 学生, 都算是对认植物有一定兴趣的爱好者。 有的爱好者对植物分类学有很深的造诣, 和专业人士比起来毫不逊色。 在瑞典, 这样的调查在很多省份都已经完成, 这样的调查为更好的认识物种分布、 物种入侵、全球变化等, 提供了最为详实的基础资料。

调查虽然只有五天, 我却深深折服于拉普兰的美, 那里的蓝天, 白云, 绿草, 那里清澈见底的河流, 广阔的沼泽, 深深映入人心。 除此之外, 也不得不折服这些瑞典人对工作的敬业与坚持: 一个基本植被调查项目, 在没有收益的情况下进行了三十年, 实属不易。再看看他们每天工作的那种执着精神,每人胸前挂着一个高倍的放大镜,讨论一个物种最细致入微的区别,一丝一毫也不放过,工作到这种状态,让人觉得很可爱。那种细致入微的工作方式,是一份纯真,也是一种美好,这些也许更耐人寻味。

img

调查用的植物志

img

调查用的地图和表格

img

租用的工作场地: 人民活动中心

img

img

img

驯鹿

img

林间小路

img

湖泊

img

水毛茛, 毛茛科毛茛属植物,随水漂浮, 沉水叶裂片线型, 细长, 挺水叶裂片掌状分裂。花白色

img

湖泊

img

img

img

img

矢车菊

img

img

img

img

img

img

img

img

欧洲赤松

img

img

欧洲赤松

img

img

沼生马先蒿, 玄参科马先蒿属植物, 花粉红色,常分布在河岸或溪流旁。

img

欧洲赤松

img

img

睡菜,为龙胆科睡菜属植物, 三出复叶,椭圆形, 喜欢生长在浅水中, 在沼泽中也有分布。

img

img

img

欧洲赤松

img

狐尾藻, 小二仙草科植物, 完全沉在水下,随水流摇摆, 仅开花时, 部分出露到水面。

img

睡莲

img

球尾花,报春花科珍珠菜属植物, 喜欢生长在河边,溪流旁。

img

img

发草

img

img

驯鹿

img

img

睡莲

img

img

睡莲

img

萍蓬草

img

萍蓬草

img

img

阵雨袭来

img

雨过天晴

img

img

img

img

北极花,又名林耐草,是忍冬科小灌木, 在我国的大兴安岭, 长白山等地有分布。 北极花纯白淡雅, 小巧精致,惹人喜爱。分类学之父,瑞典人卡尔-林耐用自己的姓氏为这个属命名。

img

麦瓶草

img

img

img

img

img

img

img

img

img

img

The bumblebee beetle

img

img

img

雏菊

img

img

狸藻,是著名的食虫植物, 喜欢生长在有一定流水的沼泽中,沉水叶中有捕虫囊, 能够捕捉水蚤等多种生物

img

img

小河截弯取直形成的牛轭湖

img

长叶茅膏菜, 只生长在贫瘠的沼泽中

img

溪流中的苔藓

img

夏季小屋,瑞典人一般在城镇拥有一套房屋, 而在乡间也有专门用于度假的房屋,仅供夏季休闲娱乐

img

沼泽

img

欧石南

img

img

羊群

img

img

img

img

img

img

奶牛

img

瑞典大量使用机械, 这是牧场上打包的青贮饲料。

How to compute phylobeta diversity indices?

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
## Files required

## phylocom.exe
## phylomatic.exe
## davies_dated.new
## species.table.csv Format shown below:
Family Genus Species
Actinidiaceae Actinidia 1
Actinidiaceae Actinidia 2
Apiaceae Aegopodium 6
...

## sample.csv
01 1 sp1
01 1 sp2
01 1 sp6
01 1 sp10
01 1 sp12
01 1 sp13
...



setwd("C:/Jinlong/phylobeta")
library(picante)
library(simba)
library(hydroTSM)

taxatab <- read.csv("species.table.csv", header = TRUE)

taxatab$SPECIES <- paste("sp", taxatab$SPECIES, sep = "")
taxatab$GENUS <- tolower(taxatab$GENUS)
taxatab$DAVIES_Family <- tolower(taxatab$DAVIES_Family)
write.table(unique(taxatab), "taxatable.txt", row.names = FALSE, col.names = FALSE, quote = FALSE, sep = "/")

shell("phylomatic -f davies_dated.new -t taxatable.txt>phylo.to.modify")
## Delete the outer most edge
dat <- read.csv("sample.csv", header = TRUE)
if(any(nchar(dat$PLOT) == 1)){
dat$PLOT[nchar(dat$PLOT) == 1] <- paste("0", dat$PLOT[nchar(dat$PLOT) == 1], sep = "")
}
dat$PLOT <- paste(dat$PLOT)
dat$SPECIES <- paste("sp", dat$SPECIES, sep = "")
dat$Abund <- rep(1, nrow(dat))
write.table(dat[,c(1, 3, 2)], "sample", quote = FALSE, row.names = FALSE, col.names = FALSE)


## Calculating Diversity metrics using picante
test.tree <- read.tree("phylo")
test.sample <- sample2matrix(dat[,c(1, 3, 2)])

## PhyloSor ## Not available from Phylocom
test.phylosor <- phylosor(samp =test.sample, tree =test.tree)

## UniFrac ## Not available from Phylocom
test.unifrac <- unifrac(comm =test.sample, tree = test.tree)

## Dnn COMDIST
test.comdist <- comdist(comm =test.sample, dis = cophenetictest.tree))

## Dpw comdistnt
test.comdistnt <- comdistnt(comm =test.sample, dis = cophenetictest.tree))

## Rao D and Rao H
test.rao <- raoD(comm =test.sample, phy =test.tree)

## raoD.Dkl ("US.species.raoD.Dkl.pairlist.ver2.csv ") is rao's D
test.raoD <- as.dist(test.rao$Dkl)

# and "US.species.raoD.H.pairlist.ver2.csv" is Rao's H.
test.raoH <- as.dist(test.rao$H)

test.phylosor.list <- liste(test.phylosor)
test.unifrac.list <- liste(test.unifrac)
test.comdist.list <- liste(test.comdist)
test.comdistnt.list <- liste(test.comdistnt)
test.rao.D.list <- liste(test.raoD )
test.rao.H.list <- liste(test.raoH)


result <- data.frame(test.phylosor.list[,3], test.unifrac.list[,3],
test.comdist.list[,3], test.comdistnt.list[,3],
test.rao.D.list[,3], test.rao.H.list[,3])

colnames(result) <- c("phylosor", "unifrac", "Dpw", "Dnn", "raoD", "raoH")

pdf("phylo.beta.pairs.pdf", width = 10, height = 10)
hydropairs(result)
dev.off()

result2 <- data.frame(test.phylosor.list, test.unifrac.list[,3],
test.comdist.list[,3], test.comdistnt.list[,3],
test.rao.D.list[,3], test.rao.H.list[,3])
colnames(result2) <- c("PlotI","PlotII", "phylosor", "unifrac", "comdist", "comdistnt", "raoD", "raoH")

write.csv(result2, "phylobeta.result.csv")


write.csv(test.phylosor.list ,"test.phylosor.list.csv" , row.names = FALSE)
write.csv(test.unifrac.list ,"test.unifrac.list.csv" , row.names = FALSE)
write.csv(test.comdist.list ,"test.Dpw.list.csv" , row.names = FALSE)
write.csv(test.comdistnt.list ,"test.Dnn.list.csv" , row.names = FALSE)
write.csv(test.rao.D.list ,"test.rao.D.list.csv" , row.names = FALSE)
write.csv(test.rao.H.list ,"test.rao.H.list.csv" , row.names = FALSE)

论文投稿用什么软件修改图片?

学术期刊要求在投稿时, 一般要求矢量图(eps, pdf, svg等)或者高分辨率(600DPI以上)的栅格图(一般为TIFF)。

这里简要总结如下:

1 生成图片

绘制用散点图、柱状图、 折线图等,一般用R (免费)、GNUplot (免费)、Python(免费)、Matlab (需购买)、 SigmaPlot(需购买)、Excel等;绘制进化树,使用FigTree (免费);绘制地图等用ArcGIS (需购买); 绘制流程图等,用OpenOffice的Draw(免费)。处理公式, 常用MathType(需购买), 也可以用EqualX Equation Editor (免费)。

2 不同格式图的修改

修改矢量图,一般保存为pdf格式或eps格式,以pdf最方便处理和浏览。

eps、 pdf等格式的矢量图, 可以用inkscape (免费)、Adobe Illustrator (需购买) 修改。

tiff等栅格图(包括照片), 一般需要在600dpi以上,可以用paint.net (免费)、 GIMP (免费)或Adobe Photoshop (需购买)修改,

总之, R、 OpenOffice、 Inkscape、GIMP、paint.net、EqualX Equation Editor 等免费软件, 完全可以保证对图片作出符合投稿要求的修改。 科研工作者几乎不需要购买昂贵的软件或使用盗版软件, 即可获得高质量的图片, 更好地展示自己的科研成果。

[转载]被子植物大于500个种的属

摘编自 维基百科 http://en.wikipedia.org/wiki/List_of_the_largest_genera_of_flowering_plants

序号 Genus 中文名 Species 中文科名 Family
1 Astragalus 黄芪属 3,270 豆科 Fabaceae / Leguminosae
2 Bulbophyllum 石豆兰属 2,032 兰科 Orchidaceae
3 Psychotria 九节属 1,951 茜草科 Rubiaceae
4 Euphorbia 大戟属 1,836 大戟科 Euphorbiaceae
5 Carex 薹草属 1,795 莎草科 Cyperaceae
6 Begonia 秋海棠属 1,484 秋海棠科 Begoniaceae
7 Dendrobium 石斛属 1,371 兰科 Orchidaceae
8 Acacia 金合欢属 c. 1,353 豆科 Fabaceae / Leguminosae
9 Solanum 茄属 c. 1,250 茄科 Solanaceae
10 Senecio 千里光属 c. 1,250 菊科 Asteraceae / Compositae
11 Croton 巴豆属 1,223 大戟科 Euphorbiaceae
12 Pleurothallis 叶上花属 1,120+ 兰科 Orchidaceae
13 Eugenia 番樱桃属 1,113 紫金牛科 Myrtaceae
14 Piper 胡椒属 1,055 胡椒科 Piperaceae
15 Ardisia 紫金牛属 1,046 紫金牛科 Myrsinaceae
16 Syzygium 蒲桃属 1,041 桃金娘科 Myrtaceae
17 Rhododendron 杜鹃花属 c. 1,000 杜鹃花科 Ericaceae
18 Miconia 大叶野牡丹属 1,000 野牡丹科 Melastomataceae
19 Peperomia 椒草属 1,000 胡椒科 Piperaceae
20 Salvia 鼠尾草属 945 唇形科 Lamiaceae / Labiatae
21 Erica 欧石楠属 860 杜鹃花科 Ericaceae
22 Impatiens 凤仙花属 850 凤仙花科 Balsaminaceae
23 Cyperus 莎草属 839 莎草科 Cyperaceae
24 Phyllanthus 叶下珠属 833 大戟科 Phyllanthaceae
25 Allium 葱属 815 葱科 Amaryllidaceae
26 Epidendrum 树兰属 800 兰科 Orchidaceae
27 Vernonia 斑鸠菊属 800–1,000 菊科 Asteraceae / Compositae
28 Lepanthes 细叶兰属 c. 800 兰科 Orchidaceae
29 Anthurium 安祖花属 789 天南星科 Araceae
30 Diospyros 柿属 767 柿树科 Ebenaceae
31 Ficus 榕属 750 桑科 Moraceae
33 Indigofera 胡枝子属 700+ 豆科 Fabaceae / Leguminosae
32 Silene 蝇子草属 700 石竹科 Caryophyllaceae
34 Oxalis 酢浆草属 700 酢浆草科 Oxalidaceae
35 Crotalaria 猪屎豆属 699 豆科 Fabaceae / Leguminosae
36 Centaurea 矢车菊属 695 菊科 Asteraceae / Compositae
37 Cassia 决明属 692 豆科 Fabaceae / Leguminosae
38 Eucalyptus 桉属 681 桃金娘科 Myrtaceae
39 Oncidium 文心兰属 680 兰科 Orchidaceae
40 Galium 拉拉藤属 661 茜草科 Rubiaceae
41 Cousinia 刺头菊属 655 菊科 Asteraceae / Compositae
42 Ipomoea 番薯属 650 旋花科 Convolvulaceae
43 Dioscorea 薯蓣属 631 薯蓣科 Dioscoreaceae
44 Cyrtandra 浆果苣苔属 622 苦苣苔科 Gesneriaceae
45 Helichrysum 蜡菊属 600 菊科 Asteraceae / Compositae
46 Ranunculus 毛茛属 600 毛茛科 Ranunculaceae
47 Habenaria 玉凤花属 600 兰科 Orchidaceae
48 Justicia 爵床属 600 爵床科 Acanthaceae
49 Schefflera 鹅掌柴属 584 五加科 Araliaceae
50 Ixora 龙船花属 561 茜草科 Rubiaceae
51 Berberis 小檗属 556 小檗科 Berberidaceae
52 Quercus 栎属 531 壳斗科 Fagaceae
53 Pandanus 露兜树属 c. 520 露兜树科 Pandanaceae
54 Panicum 黍属 500+ 禾本科 Poaceae / Gramineae
55 Eria 毛兰属 500 兰科 Orchidaceae
56 Polygala 远志属 500 远志科 Polygalaceae
57 Potentilla 委陵菜属 500 蔷薇科 Rosaceae

R中的极大似然估计

人们经常用样本(sample)来估计总体(population)的参数, 假定一个概率密度函数,如正态分布,并给出所有可能的参数组合,那么每一种参数组合所生成当前样本的可能性分别有多大?将最可能产生当前样本的参数,作为总体参数的估计,这就是极大似然(Maximum Likelihood)思想的体现。

在实际的计算中, 人们往往需要先假定某一概率密度, 如正态分布,

$p.d.f. (y) = \frac{1} {\sqrt{2 \pi} \sigma} e^{-1/2 \frac{(y- \mu)^2} {\sigma^2}}$

1
pdfi <- 1/(sqrt(2*pi)*sigma) * exp(-1/2 * (obs - mu)^2/(sigma^2))

其中,

  • obs为某一次观测值
  • mu 为 总体平均值的任意估计值
  • sigma 为 总体标准差的任意估计值,
    则可以求得观测到该次观测值的概率 pdfi

对所有的样本数据的概率密度相成, 即获得联合概率密度, 即似然函数。
似然函数 L = pdf1 * pdf2 * pdf3 * ... * pdfn

$L = \prod_{i=1}^n \frac{1} {\sqrt{2 \pi} \sigma} e^{-1/2 \frac{(y- \mu)^2} {\sigma^2}}$

那么,什么样的参数, 才能保证Likelihood函数取最大值?
由于联合概率密度的值往往非常小,在计算机中难以处理,一般情况下, 将似然函数取对数, 便于数值操作。 此时,做一系列的变换。 :
如:

$LogL = \sum_{i=1}^n \log\left ( \frac{1} {\sqrt{2 \pi} \sigma} e^{-1/2 \frac{(y- \mu)^2} {\sigma^2}} \right )$

$LogL= \sum_{i=1}^n \left [ \log\left ( \frac{1} {\sqrt{2 \pi} \sigma} \right ) + \log\left ( e^{-1/2 \frac{(y- \mu)^2} {\sigma^2}} \right ) \right ]$

$LogL= \sum_{i=1}^n \log(2 \pi)^{-1/2} + \sum_{i=1}^n \log\left (\frac{1} {\sigma}\right ) + \sum_{i=1}^n \left ( -1/2 \frac{(y- \mu)^2} {\sigma^2} \right )$

$LogL= -\frac{n} {2} \log(2 \pi) - n \log(\sigma) - \frac{1} {2 \sigma^2} \sum_{i=1}^n ( y - \mu)^2$

若认为LogL为$\mu$的函数, 若保证LogL最大值,令 $\sum_{i=1}^n ( y - \mu)^2$ 取最大值,则

$\mu = \frac{\sum_{i=1}^n y_i}{n}$
另外一方面, 若认为LogL为$\sigma$的函数, 将LogL函数对$\sigma$求导, 令导数为0, 可求得极值。

$\frac{\delta LogL}{\delta \sigma} = -\frac{n}{\sigma} - (-2) \frac{1}{2} \sum_{i=1}^n (y - \mu)^2 \sigma^{-3}$

$\frac{\delta LogL}{\delta \sigma} = -\frac{n}{\sigma} + \frac{1}{\sigma^2} \sum_{i=1}^n (y - \mu)^2$

令导数为0, 得
$\frac{n}{\sigma} = \frac{1}{\sigma^3} \sum_{i=1}^n (y - \mu)^2$

$\sigma^2 = \frac{\sum_{i=1}^n (y - \mu)^2}{n}$

此时, 带入之前$\mu$,求得$\sigma$的估计

$\hat{\sigma}^2 = \frac{\sum_{i=1}^n (y - \bar{y})^2}{n}$

用极大似然估计求得的$\sigma$值为有偏估计, 即样本的最可能取值。

以上是解析解。 在实际情况下, 人们常用多种优化方法, 求得极大似然估计的数值解。 这些方法包括

  • "Nelder-Mead" (Nelder and Mead ,1965),
  • "BFGS" (Broyden, Fletcher, Goldfarb and Shanno, 1970),
  • "CG" (Fletcher and Reeves,1964),
  • "L-BFGS-B" (Byrd et. al. ,1995) ,
  • "SANN" (Belisle 1992), 等。

由此看来, 极大似然估计, 首先是按照假定的概率密度, 求得联合概率密度, 即似然函数,进一步写成对数似然函数, 再用数值优化的方法, 求得对数似然函数所对应的点(各参数)。即获得参数的极大然估计。

下面提供一个实例:
问题:
假设一批树木中抽取若干个体,胸径分别为 150, 124, 100, 107, 170, 144, 113, 108, 92, 129, 123, 118,试求这批树木胸径的平均值和标准差。

分析:
如果用这批样本数据直接计算平均值和方差, 可以得到总体平均值和标准差的估计。

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
dat <- c(150, 124, 100, 107, 170, 144, 113, 108, 92, 129, 123, 118)

> mean(dat)
[1] 123.1667

### 求样本平均值
mean.sample <- sum(dat)/length(dat)

### 设置空向量
var.sample <- c()
for(i in 1:length(dat)){
var.sample[i] <- (dat[i] - mean.sample)^2
}
dd <- sum(var.sample)

## 求样本标准差
sqrt(dd/(length(dat)))

## 求总体标准差
sqrt(dd/(length(dat) - 1))
以上是普通的估计方法


下面用极大似然方法估计
###########################################################################
### 联合概率密度

like = function(data, mu, sigma) {
like = 1
for(i in 1:length(data)){
# like = like *( 1/(sqrt(2*pi)*sigma) * exp(-1/2 * (data[i] - mu)^2/(sigma^2)))
like = like * dnorm(x = dat[i], mean = mu, sd = sigma)
}
return(like)
}

like(dat, 120, 20)

### 对数似然函数
loglike = function(data, mu, sigma) {
loglike = 0
for(i in 1:length(data)){
loglike = loglike + log(dnorm(x = dat[i], mean = mu, sd = sigma))
}
return(loglike)
}
loglike(dat, 122, 20)

### 扩展参数组合
params = expand.grid(mu = seq(50, 200, 1),
sigma = seq(10, 40, 1))
LogLik.surf = with(params, loglike(dat, mu, sigma))

dim(LogLik.surf) <- c(length(seq(50, 200, 1)), length(seq(10, 40, 1)))

library(lattice)
contourplot(LogLik.surf ~ mu*sigma, data = params, cuts = 20)

# Zooming in
params = expand.grid(mu = seq(120, 126, 0.01), sigma = seq(18, 25, 0.1))
LogLik.surf.zoom = with(params, loglike(dat, mu, sigma))

dim(LogLik.surf.zoom) <- c(length(seq(120, 126, 0.01)), length(seq(18, 25, 0.1)))

library(lattice)
contourplot(LogLik.surf.zoom ~ mu*sigma, data = params, cuts = 20, region = TRUE)

### nlm进行极大似然优化, 需要提供初始值
### 对数似然函数
loglike000 = function(p) {
loglike = 0
for(i in 1:length(dat)){
loglike = loglike + log(dnorm(x = dat[i], mean = p[1], sd = p[2]))
}
return(-loglike)
}

(MLEst <- nlm(f = loglike000, p = c(100, 10)))

参考

http://www.quantumforest.com/2011/10/maximum-likelihood/