线性模型参数的极大似然估计

要了解响应变量y与因变量x之间的关系， 人们常用到线性回归。 常用 y = a*x + b表示， 然而更一般格式为

y = beta1*x + beta0 + error

其中beta0， beta1是参数,分别为截距和因变量的系数，error为残差。按照线性模型的假设， 残差应符合正态分布，一般取均值mean为0， 只需要估计标准差sd即可。 。

本文的目的是用R程序的实例说明如何用极大似然估计进行简单的线性回归的参数。 按照以下步骤完成：

因本文需要用到bbmle程序包， 没有安装的读者， 需要用 install.packages("bbmle")安装bbmle程序包。

步骤：

  1. 生成要拟合的数据， 并绘图
  2. 通过R函数的内置lm方法建立线性模型， 并计算AIC以及LogLikelihood以便进行比较
  3. 建立对数似然函数 Log Likelihood
  4. 利用bbmle程序包的 mle2函数， 用数值优化算法(Optimization)求似然函数的最大值，从而给出所求参数的估计值。
  5. 基于最大对数似然值 Maximum Log Likelihood求该模型的AIC值（Akaike Information Criterion）
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
library(bbmle)  ### 导入bbmle程序包
set.seed(12345) ### 设定随机数种子
x = 1:60        ### 生成x数据， 共60个值
y = x*5 + 20 + rnorm(60, 0, 20) ## 生成y数据， 这里假设已知 beta1 = 5， beta0 = 20， mu = 0, sd = 20
plot(y ~ x, main = "Fitting of a Linear Model ") ## 观察数据
linmod <- lm(y ~ x) ## 建立线性模型
lines(x, fitted(linmod) , lty = 1, col = "red", lwd = 2)  ## 添加线段
segments(x0 = x, y0 = y, x1 = x, y1 = fitted(linmod))  ## 添加残差线
summary(linmod)  ## 查看模型的拟合参数
(AIC1 <- AIC(linmod))  ## Akaike Information Criterion
(logLik1 <- logLik(linmod)) ## Loglikelihood

##################################################
### Log Likelihood Function
### 函数的一般形式
### y = beta1*x + beta0 + error
### 而error一项， 服从正态分布
### 所以 error = y - beta1*x - beta0
### 所以Loglikelihood函数可以写成如下形式：
LL <- function(beta0, beta1, mu, sigma){
   R = y - x*beta1 - beta0  
   -sum(dnorm(R, mu, sigma, log = TRUE))
}

### mle2函数， 可以用来估计似然函数去极值时，各参数的取值， 调用方式如下：
res <- mle2(LL, start = list(beta0 = 1, beta1 = 2, mu = 1, sigma = 1), fixed = list(mu = 0))
summary(res)

### 因为mle2的结果是S4类型， 所以这里用@提取
### 基于参数估计，计算y2的取值
y2 = x*(res@coef[2]) + res@coef[1]
### 添加极大似然估计值
points(y2 ~ x, col= "blue", type = "b", pch = 19)
### 用mle2得出来的Maximum Log Likelihood
(logLik2 <- logLik(res))
### 计算AIC, 的公式为 AIC = -2*logLik + 2*p, 其中logLike为 Maximum Log Likelihood， p为参数的个数。
### 因为在参数估计中， mu为一个常数0，我们只估计了sigma， 所以估计的参数是3个
(AIC2 = -2*logLik(res) + 2*3)

R拟合统计分布并绘制曲线

R拟合统计分布并绘制曲线

在统计分析中， 经常会用到统计分布拟合， 这里介绍MASS程序包的 fitdistr(), 该函数可以用来拟合 "beta", "cauchy", "chi-squared", "exponential", "f", "gamma", "geometric", "log-normal", "lognormal", "logistic", "negative binomial", "normal", "Poisson", "t" 以及 "weibull" 的概率密度分布， 现举例如下：

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
115
116
117
118
119
120
121
122
123
124
125
126
127
128
library(MASS)

#### Gammadistribution
x <-rgamma(100, shape =5, rate =0.1)
res <- fitdistr(x, "gamma")
h <-hist(x)
xfit <-seq(min(x), max(x), by=(max(x)-min(x))/100)
yfit <-dgamma(xfit,res[[1]][1],res[[1]][2])
yfit <- yfit*diff(h$mids[1:2])*length(xfit)
lines(xfit, yfit, col="blue", lwd=2)

####"beta"                
x <-rbeta(100, shape1 =5, shape2 =2)
res <- fitdistr(x, "beta", start=list(shape1 =4, shape2 =1))
h <-hist(x)
xfit <-seq(min(x), max(x), by=(max(x)-min(x))/100)
yfit <-dbeta(xfit,res[[1]][1],res[[1]][2])
yfit <- yfit*diff(h$mids[1:2])*length(xfit)
lines(xfit, yfit, col="blue", lwd=2)

####"cauchy",            
"cauchy",
x <-rcauchy(100, location =200, scale=300)
res <- fitdistr(x, "cauchy")
h <-hist(x)
xfit <-seq(min(x), max(x), by=(max(x)-min(x))/100)
yfit <-dcauchy(xfit,res[[1]][1],res[[1]][2])
yfit <- yfit*diff(h$mids[1:2])*length(xfit)
lines(xfit, yfit, col="blue", lwd=2)

####"chi-squared",        
"chi-squared",
x <-rchisq(100, df=5, ncp =2)
res <- fitdistr(x, "chi-squared", start=list(df=2, ncp =1))
h <-hist(x)
xfit <-seq(min(x), max(x), by=(max(x)-min(x))/100)
yfit <-dchisq(xfit,res[[1]][1],res[[1]][2])
yfit <- yfit*diff(h$mids[1:2])*length(xfit)
lines(xfit, yfit, col="blue", lwd=2)

####"exponential",        
"exponential"
x <-rexp(100, rate =10)
res <- fitdistr(x, "exponential")
h <-hist(x)
xfit <-seq(min(x), max(x), by=(max(x)-min(x))/100)
yfit <-dexp(xfit,res[[1]][1])
yfit <- yfit*diff(h$mids[1:2])*length(xfit)
lines(xfit, yfit, col="blue", lwd=2)

####"f",                  
"f"
x <-rf(100, df1 =60, df2 =50)
res <- fitdistr(x, "f", start=list(df1 =10, df2 =2))
h <-hist(x)
xfit <-seq(min(x), max(x), by=(max(x)-min(x))/100)
yfit <-df(xfit,res[[1]][1], res[[1]][2])
yfit <- yfit*diff(h$mids[1:2])*length(xfit)
lines(xfit, yfit, col="blue", lwd=2)

####"lognormal",          
"lognormal"
x <-rlnorm(100, meanlog =5, sdlog =1.5)
res <- fitdistr(x, "lognormal")
h <-hist(x)
xfit <-seq(min(x), max(x), by=(max(x)-min(x))/100)
yfit <-dlnorm(xfit,res[[1]][1], res[[1]][2])
yfit <- yfit*diff(h$mids[1:2])*length(xfit)
lines(xfit, yfit, col="blue", lwd=2)

####"logistic",          
"logistic",
x <-rlogis(100, location =5, scale=2)
res <- fitdistr(x, "logistic", start=list(location=5, scale=2))
h <-hist(x)
xfit <-seq(min(x), max(x), by=(max(x)-min(x))/100)
yfit <-dlogis(xfit,res[[1]][1], res[[1]][2])
yfit <- yfit*diff(h$mids[1:2])*length(xfit)
lines(xfit, yfit, col="blue", lwd=2)

#### "negativebinomial",  
"negative binomial",
x <-rnbinom(100, size =300, prob =.3)
res <- fitdistr(x, "negative binomial")
h <-hist(x)
xfit <-floor(seq(min(x), max(x), by=(max(x)-min(x))/100))
yfit <-dnbinom(xfit,size = res[[1]][1], mu = res[[1]][2])
yfit <- yfit*diff(h$mids[1:2])*length(xfit)
lines(xfit, yfit, col="blue", lwd=2)
####"normal",            
"normal"
x <-rnorm(100, mean=15, sd=2)
res <- fitdistr(x, "normal")
h <-hist(x)
xfit <-seq(min(x), max(x), by=(max(x)-min(x))/100)
yfit <-dnorm(xfit, mean= res[[1]][1], sd= res[[1]][2])
yfit <- yfit*diff(h$mids[1:2])*length(xfit)
lines(xfit, yfit, col="blue", lwd=2)

####"Poisson",            
"Poisson",
x <-rpois(100, lambda =400)
res <- fitdistr(x, "Poisson")
h <-hist(x)
xfit <-floor(seq(min(x), max(x), by=(max(x)-min(x))/100))
yfit <-dpois(xfit,res[[1]][1])
yfit <- yfit*diff(h$mids[1:2])*length(xfit)
lines(xfit, yfit, col="blue", lwd=2)

####"t"                  
"t"
x <-rt(100, df=5)
res <- fitdistr(x, "t")
h <-hist(x)
xfit <-seq(min(x), max(x), by=(max(x)-min(x))/100)
yfit <-dt(xfit,res[[1]][3])  
yfit <- yfit*diff(h$mids[1:2])*length(xfit)
lines(xfit, yfit, col="blue", lwd=2)

####"weibull"            
"weibull"
x <-rweibull(100, shape =5, scale=9)
res <- fitdistr(x, "weibull", start=list(shape =1, scale=1))
h <-hist(x)
xfit <-seq(min(x), max(x), by=(max(x)-min(x))/100)
yfit <-dweibull(xfit,res[[1]][1], res[[1]][2])
yfit <- yfit*diff(h$mids[1:2])*length(xfit)
lines(xfit, yfit, col="blue", lwd=2)

R程序包整合C++：Rcpp使用简明指南

R语言轻巧简洁，便于开发新算法。但是R函数运行速度较慢， 在进行迭代计算或者跑蒙特卡罗马尔科夫链时，常常不能满足要求。 此外， 由于用C或者C++已经有很多算法库，人们希望能直接调用这些库中的函数。Rcpp正是为了让R开发者更方便地运用C++函数生的。 本文简要介绍如何编写一个含有C++函数的R程序包。

要在Windows下创建R程序包，首先要安装：

然后，还需要安装 Rcpp （http://cran.r-project.org/web/packages/Rcpp/index.html ）程序包。

并配置好启动 路径， 即将所需要的任何exe文件的路径， 拷贝到 电脑>属性>高级系统设置>高级>环境变量>系统变量>路径 中 如:

此外， 这里假设读者使用文本编辑器 Notepad++ 。

Rcpp程序包已经提供了多种C++的类模板， 可以让C++函数返回多种R对象。假设现有一个C++函数， 名叫Eccentricity，该函数用于计算某天体的轨道偏心率随时间的变化，输入参数为儒略日（日期的某种记录方法）, 在正常的C++语法中， 该函数定义如下：

1
2
3
4
5
6
7
8
double Eccentricity(double JD)

{
double JD3  = JD;
double T = (JD3 - 2451545.0) / 36525.0;
double Tsquared = T*T;
return ((1 - 0.002516*T - 0.0000074*Tsquared));
}

现在我们需要有一个Rcpp形成的包骨架(package skeleton)，假设名叫 skycalc。在R中运行如下命令生成基于Rcpp的R包骨架:

1
2
3
library(Rcpp)
Rcpp.package.skeleton("skycalc")
getwd()

skycalc保存在getwd()所示的文件夹下， 有 man、R、src三个文件夹以及DESCRIPTION、NAMESPACE 以及 Read-and-delete-me三个文本文件，这个三个文件没有任何扩展名。

  • man文件夹保存的是所有R函数以及R程序包的帮助文件文件，遵循是Latex格式，需要逐项填写。
  • R文件夹保存的是R函数的文件，每个R函数单独一个文件
  • src文件夹放的是cpp文件， 即c++的源文件

默认的情况下，各文件夹下有以下示例文件：

  • man文件夹中包含 rcpp_hello_world.Rd 文件，
  • R文件夹中包含RcppExports.R文件
  • src 文件夹中包含rcpp_hello_world.cpp 和 RcppExports.cpp 两个C++源文件

这些示例文件在本例中均可以删除。

我们的目标是：将Eccentricity函数，放在C++源代码中。假设函数保存在名为Eccentricity.cpp的纯文本文件， 为了能够在R中直接调用该函数， 需要将其参数数据类型以及返回值的数据类型更改为R能够识别的格式。因为R只能识别一种称为SEXP的数据类型， 所以这里做如下更改：

1
2
3
4
5
6
7
RcppExport SEXP Eccentricity(SEXP JD)
{
double JD3  = Rcpp::as<double>(JD);
double T = (JD3 - 2451545.0) / 36525.0;
double Tsquared = T*T;
return (wrap(1 - 0.002516*T - 0.0000074*Tsquared));
}

在类型前，更加上RcppExport，返回值return里面，先套用wrap函数，以便返回为R能识别的SEXP数据类型。

为在Eccentricity.cpp文件中加上引用的头文件：

1
2
3
#include <Rcpp.h>
using namespace Rcpp;
using namespace std;

此时， Eccentricity.cpp文件内容如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
///////////////cpp文件开始//////////////////////

#include <Rcpp.h>
using namespace Rcpp;
using namespace std;
RcppExport SEXP Eccentricity(SEXP JD)
{
double JD3  = Rcpp::as<double>(JD);
double T = (JD3 - 2451545.0) / 36525.0;
double Tsquared = T*T;
return (wrap(1 - 0.002516*T - 0.0000074*Tsquared));
}

///////////////cpp文件结束//////////////////////

这样 Eccentricity.cpp 文件就准备好了，将其保存在src文件夹中。

为了能够在R中调用编译后的C++函数，可使用R的.Call()函数调用dll库（cpp文件经过编译后， 在Windows中编译为dll文件）。但为了更为方便地使用这些函数，一般为编写一个与C++函数形式上类似的R函数，该R函数的名称和各参数最好能与C++函数对应。为此，我们可以写一个R函数EccentricityR， 该R函数可定义如下：

1
2
3
EccentricityR <- function(JD) {
   .Call('Eccentricity', JD, PACKAGE = 'skycalc')
}

其中JD为该R函数的参数， 直接通过.Call传递给DLL中的Eccentricity。

在Notepad++中，拷贝以上R函数， 并且另存为 EccentricityR.R 文件， 放在skycalc程序包下的R文件夹下。

为了生成该EccentricityR函数的帮助文件，需要将该函数粘贴到Rconsole中， 之后运行R命令 prompt(EccentricityR)，所形成的EccentricityR.Rd 文件需要进一步根据参数的意义填写和编辑。之后放到man文件夹下。

编辑和填写DESCRIPTION以及NAMESPACE两个文本文件 (参考 http://blog.sciencenet.cn/blog-255662-247614.html )。删除 Read-and-delete-me文件。这样，程序包的源文件就做好了。

下面准备编译或安装R程序包的Window批处理文件，bat文件

R包的检查的bat文件

新建一个纯文本文件，并将扩展名txt更改为.bat，将以下内容拷贝到该文件夹中：

1
2
Rcmd check skycalc
PAUSE

将该文件命名为 check skycalc package.bat， 双击该bat文件， 可以对skycalc程序包中的错误进行检查。 如果要提交程序包到CRAN，检查过程中不能有任何错误或者warning。

建立Windows版本的安装包的bat文件

(2)新建一个纯文本文件，并将扩展名txt更改为.bat

将以下内容拷贝到该文件夹中

1
2
Rcmd INSTALL --build skycalc
PAUSE

将该文件命名为 build skycalc package Windows Binary.bat， 双击该文件， 可以建立Windows系统下的R程序包。

创建源代码包的bat文件

新建一个纯文本文件，并将扩展名txt更改为.bat，将以下内容拷贝到该文件夹中

1
2
Rcmd build skycalc
pause

将该文件命名为 build skycalc package Linux Source Code.bat， 双击该文件，可以建立Linux系统下的安装包，其实是源文件。

直接安装程序包的bat文件

新建一个纯文本文件，并将扩展名txt更改为.bat，将以下内容拷贝到该文件夹中：

1
2
Rcmd INSTALL skycalc
PAUSE

将该文件命名为 install skycalc package.bat， 双击该文件， 可以安装skycalc到当前的R中。

运行批处理命令

将以上四个.bat文件， 放置到skycalc文件夹所在的文件夹（例如， skycalc所 在的文件夹为study， 则bat文件需要放置在study文件夹下）

双击文件， 即可完成程序包的编译和检查等。

注意： 如果使用的64bit的Windows，源代码 skycalc文件夹下， 每次都会生成 src-x64 和 src-i386两个文件夹， 在编译linux源代码的时候， 需要删除。

进一步阅读

在频率直方图上添加拟合的曲线－以对数正态分布为例

以下为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
\##### 第一个方案
\##### 拟合正态分布
x <- log(rlnorm(1000))*500 ### 提取数据框
h <- hist(x, ylim = c(0, 300))  ## 绘制频率直方图
xfit <- seq(min(x),max(x),length=500) ### 取要计算的x值
yfit <- dnorm(xfit, mean=mean(x), sd=sd(x)) ### 假设是正态分布， 则对应的概率密度为 yfit
yfit <- yfit*diff(h$mids[1:2])*length(x) ### 求x向量内每一个值与前面数值的差, 以确定在y方向上， 应该放大多少倍
lines(xfit, yfit, col="blue", lwd=2)  ### 绘制曲线


\##### 第二个方案
\##### 拟合对数正态分布
x <- rlnorm(1000)*1000  ### 提取第一列
h <- hist(x, ylim = c(0, 1000), breaks = 25)  ## 绘制频率直方图
xfit <- seq(min(x),max(x),length=500) ### 取要绘图的x值取值范围
yfit <- dlnorm(xfit, meanlog=log(mean(x)), sdlog=sd(log(x))) ### 假设是对数正态分布， 则对应的概率密度为 yfit
yfit <- yfit*diff(h$mids[1:2])*length(x) ### 求x向量内每一个值与前面数值的差, 以确定在y方向上， 应该放大多少倍
lines(xfit, yfit, col="blue", lwd=2) ### 绘制曲线


\##### 第三个方案
\##### 基于核密度估计
x <- rlnorm(1000)*1000   ### 提取第一列
h <- hist(x, breaks = 25)  ## 绘制频率直方图
xfit <- seq(min(x),max(x),length=500) ### 取要绘图的x值取值范围
yfit <- density(x, n = length(xfit))$y ### 假设是对数正态分布， 则对应的概率密度为 yfit
yfit <- yfit*diff(h$mids[1:2])*length(xfit) ### 求x向量内每一个值与前面数值的差, 以确定在y方向上， 应该放大多少倍
lines(xfit, yfit, col="blue", lwd=2) ### 绘制曲线

PyEphem基本功能介绍

PyEphem为Python下的一个程序包， 用来进行天文历算， 虽然是爱好者编写的， 但是由于使用VOS87行星运动数据， 计算精度达到了很高的精度， 足以满足一般的观测需要。 详情参见 http://rhodesmill.org/pyephem/ 本文是学习该程序包的读书笔记。

火星在2018年的星座

1
2
3
4
import ephem
m = ephem.Mars('2018')
print(ephem.constellation(m))
('Lib', 'Libra')

计算不同天体在1970年1月1日0时的位置等信息

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
from ephem import *
mars = Mars('1970') # 火星
print('%s, %s, %.10f' % (mars.name, mars.elong, mars.size))

sun = Sun('1970') # 太阳
print('%s, %s, %.10f' % (sun.name, sun.elong, sun.size))

moon = Moon('1970') # 月亮
print('%s, %s, %.10f' % (moon.name, moon.elong, moon.size))

mercury = Mercury('1970') # 水星
print('%s, %s, %.10f' % (mercury.name, mercury.elong, mercury.size))

venus = Venus('1970') # 金星
print('%s, %s, %.10f' % (venus.name, venus.elong, venus.size))

jupiter = Jupiter('1970') # 木星
print('%s, %s, %.10f' % (jupiter.name, jupiter.elong, jupiter.size))

saturn = Saturn('1970') # 土星
print('%s, %s, %.10f' % (saturn.name, saturn.elong, saturn.size))

uranus = Uranus('1970') # 天王星
print('%s, %s, %.10f' % (uranus.name, uranus.elong, uranus.size))

neptune = Neptune('1970') # 海王星
print('%s, %s, %.10f' % (neptune.name, neptune.elong, neptune.size))

pluto = Pluto('1970') # 冥王星
print('%s, %s , %.10f' % (pluto.name, pluto.elong, pluto.size))
Mars, 62:04:48.0, 5.9295825958
Sun, 0:00:00.0, 1951.8359375000
Moon, -89:27:48.3, 1831.0117187500
Mercury, 18:52:17.4, 7.5804867744
Venus, -5:42:32.4, 9.9605798721
Jupiter, -67:50:24.1, 34.2474060059
Saturn, 111:52:32.6, 18.8122367859
Uranus, -91:26:29.1, 3.8586533070
Neptune, -40:17:59.5, 2.1982953548
Pluto, -102:16:48.8 , 0.2606950402
### a_ra — Astrometric geocentric right ascension for the epoch specified
### a_dec — Astrometric geocentric declination for the epoch specified
### g_ra and ra — Apparent geocentric right ascension for the epoch-of-date
### g_dec and dec — Apparent geocentric declination for the epoch-of-date
### elong — Elongation (angle to sun)
### mag — Magnitude
### size — Size (diameter in arcseconds)
### radius — Size (radius as an angle)
### circumpolar — whether it stays above the horizon
### neverup — whether is stays below the horizon
############################################ 
hlon — Heliocentric longitude (see next paragraph)
### hlat — Heliocentric latitude (see next paragraph)
### sun_distance — Distance to Sun (AU)
### earth_distance — Distance to Earth (AU)
### phase — Percent of surface illuminated

已知观测者的位置， 包括 经纬度和海拔， 观测日期和时刻， 求天体的高度和方位

1
2
3
4
5
6
7
8
gatech = Observer()
gatech.lon = '114.169576'
gatech.lat = '22.451939'
gatech.elevation = 5
gatech.date = '2014/9/22 22:30:00'
v = Venus(gatech)
print('%s, %s' % (v.alt, v.az))
11:47:49.3, 89:41:14.6

已知彗星的根数， 计算彗星的亮度和位置， 读取的格式为 xephem格式

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
line = "C/2002 Y1 (Juels-Holvorcem),e,103.7816,166.2194,128.8232,242.5695,
0.0002609,0.99705756,0.0000,04/13.2508/2003,2000,g  6.5,4.0"
yh = ephem.readdb(line)
yh.compute('2007/10/1')
print('%.10f' % yh.earth_distance)
print(yh.mag)
14.8046731949
23.96
## 已知人造卫星的轨道根数， 计算其在任意时刻的高度和方位
### 读取 xephem格式用 readdb 方法
### 写出 xepheme格式， 用 writedb 方法
######### 人造卫星为 tle格式, 这里用readtle方法来阅读
line1 = "ISS (ZARYA)"
line2 = "1 25544U 98067A   03097.78853147  .00021906  00000-0  28403-3 0  8652"
line3 = "2 25544  51.6361  13.7980 0004256  35.6671  59.2566 15.58778559250029"
iss = ephem.readtle(line1, line2, line3)
iss.compute('2003/3/23')
print('%s %s' % (iss.sublong, iss.sublat))
-76:24:18.3 13:05:31.1
('Sgr', 'Sagittarius')
50.54
V2 - P274
V2 - P135
V1 - P273

计算月亮在某时刻所在星座

1
2
3
m = Moon('1980/6/1')
print(constellation(m))
('Sgr', 'Sagittarius')

查询某一赤经赤纬在不同星图中的页码

1
2
3
4
5
6
7
8
9
####### Uranometria by Johannes Bayer.
####### Uranometria 2000.0 edited by Wil Tirion.
#######  Millennium Star Atlas by Roger W. Sinnott and Michael A. C. Perryman.

ra, dec = '7:16:00', '-6:17:00'print(ephem.uranometria(ra, dec))

print(ephem.uranometria2000(ra, dec))

print(ephem.millennium_atlas(ra, dec))

转换儒略日

1
ephem.julian_date('2000/1/1')

计算 Delta T

1
print(ephem.delta_t('1980'))

计算两个物件的角距

1
2
3
4
5
m1 = ephem.Moon('1970/1/16')
m2 = ephem.Moon('1970/1/17')
s = ephem.separation(m1, m2)
print("In one day the Moon moved %s" % s)
In one day the Moon moved 12:33:28.5

坐标系转换

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
np = Equatorial('0', '90', epoch='2000')
g = Galactic(np)
print('%s %s' % (g.lon, g.lat))
####   Equatorial
####   ra — right ascension
####   dec — declination
####   epoch — epoch of the coordinate
####   Ecliptic
####   lon — ecliptic longitude (+E)
####   lat — ecliptic latitude (+N)
####   epoch — epoch of the coordinate
####   Galactic
####   lon — galactic longitude (+E)
####   lat — galactic latitude (+N)
####   epoch — epoch of the coordinate
122:55:54.9 27:07:41.7

计算某地点的天象

天体的高度和方位

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
#### Describes a position on Earth’s surface.lowell = ephem.Observer()
lowell.lon = '-111:32.1'
lowell.lat = '35:05.8'
lowell.elevation = 2198
lowell.date = '1986/3/13'
j = ephem.Jupiter()
j.compute(lowell)
print('%s %s' % (j.alt, j.az))
##  Attributes can be set:
##      date — Date and time
##      epoch — Epoch for astrometric RA/dec
##      lat — Latitude (+N)
##      lon — Longitude (+E)
##      elevation — Elevation (m)
##      temp — Temperature (°C)
##      pressure — Atmospheric pressure (mBar)
##  ##      The date defaults to now().
##      The epoch defaults to '2000'.
##      The temp defaults to 25°C.
##      The pressure defaults to 1010mBar.
##      Other attributes default to zero.
0:57:44.7 256:41:01.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
############## 内置的城市位置数据
city = city('Beijing')
print('%s %s' % (city.lat, city.lon))
############# 天体的升降和中天
sitka = Observer()
sitka.date = '1999/6/27'
sitka.lat = '57:10'
sitka.lon = '-135:15'm = Mars()
print(sitka.next_transit(m))
print('%s %s' % (m.alt, m.az))
print(sitka.next_rising(m, start='1999/6/28'))
print('%s %s' % (m.alt, m.az))
#### 给定观测者之后， 每个天体可以使用方法 （人造卫星除外）
### previous_transit()
### next_transit()
### previous_antitransit()
### next_antitransit()
### previous_rising()
### next_rising()
### previous_setting()
### next_setting()
39:54:15.2 116:24:26.7
1999/6/27 04:22:45
21:18:33.6 180:00:00.0
1999/6/28 23:28:25
-0:00:05.8 111:10:41.6

人造卫星升落和中天时刻等

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
##### 某年某月某日， 某地， 计算人造卫星掠过的时刻
line1 = "IRIDIUM 80 [+]"
line2 = "1 25469U 98051C   09119.61415140 -.00000218  00000-0 -84793-4 0  4781"
line3 = "2 25469  86.4029 183.4052 0002522  86.7221 273.4294 14.34215064557061"
iridium_80 = ephem.readtle(line1, line2, line3)

city.date = '2009/5/1'
info = city.next_pass(iridium_80)
print("Rise time: %s azimuth: %s" % (info[0], info[1]))
####  返回值的时刻####  
0  Rise time####  
1  Rise azimuth####  
2  Transit time####  
3  Transit altitude####  
4  Set time####  
5  Set azimuth
##############################################
Rise time: 2009/5/1 00:45:05 azimuth: 6:13:37.5

计算日出时刻

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
sun = Sun()
greenwich = Observer()
greenwich.lat = '51:28:38'
print(greenwich.horizon)
greenwich.date = '2007/10/1'
r1 = greenwich.next_rising(sun)
greenwich.pressure = 0
greenwich.horizon = '-0:34'
greenwich.date = '2007/10/1'
r2 = greenwich.next_rising(sun)
print('Visual sunrise: %s' % r1)
print('Naval Observatory sunrise: %s' % r2)
0:00:00.0
Visual sunrise: 2007/10/1 05:59:30
Naval Observatory sunrise: 2007/10/1 05:59:50

恒星时

1
2
3
4
5
6
7
madrid = ephem.city('Madrid')
madrid.date = '1978/10/3 11:32'
print(madrid.sidereal_time())
12:04:28.09
2000/3/20 07:35:17
2000/6/21 01:47:51
Spring lasted 92.8 days

春分， 秋分， 冬至， 夏至的时刻

1
2
3
4
5
6
7
8
9
10
11
12
13
d1 = ephem.next_equinox('2000')
print(d1)
d2 = ephem.next_solstice(d1)
print(d2)
t = d2 - d1
print("Spring lasted %.1f days" % t)
#### 相应的方法
#### previous_solstice()
#### next_solstice()
######## previous_equinox()
#### next_equinox()
######## previous_vernal_equinox()
#### next_vernal_equinox()

月相

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
d1 = ephem.next_full_moon('1984')
print(d1)

d2 = ephem.next_new_moon(d1)
print(d2)
### 朔
### previous_new_moon()
### next_new_moon()
### 上弦月
### previous_first_quarter_moon()
### next_first_quarter_moon()
### 望
### previous_full_moon()
### next_full_moon()
### 下弦月
### previous_last_quarter_moon()
### next_last_quarter_moon()
1984/1/18 14:05:10
1984/2/1 23:46:25

角度转换

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
a = ephem.degrees('180:00:00')
print(a)
a
print("180° is %f radians" % a)
h = ephem.hours('1:00:00')
deg = ephem.degrees(h)
print("1h right ascension = %s°" % deg)
## 时间和角度输入
ephem.degrees(ephem.pi / 32)
ephem.degrees('5.625')
ephem.degrees('5:37.5')
ephem.degrees('5:37:30')
ephem.degrees('5:37:30.0')
ephem.hours('0.375')
ephem.hours('0:22.5')
ephem.hours('0:22:30')
ephem.hours('0:22:30.0')
180:00:00.0
180° is 3.141593 radians
1h right ascension = 15:00:00.0°
#### 日期时间的输入d = ephem.Date('1997/3/9 5:13')
print(d)
d
d.triple()
d.tuple()
d + ephem.hour
print(ephem.date(d + ephem.hour))
print(ephem.date(d + 1))
### 输入日期的格式ephem.Date(35497.7197916667)
ephem.Date('1997/3/10.2197916667')
ephem.Date('1997/3/10 05.275')
ephem.Date('1997/3/10 05:16.5')
ephem.Date('1997/3/10 05:16:30')
ephem.Date('1997/3/10 05:16:30.0')
ephem.Date((1997, 3, 10.2197916667))
ephem.Date((1997, 3, 10, 5, 16, 30.0))
1997/3/9 05:13:00
1997/3/9 06:13:00
1997/3/10 05:13:00
35497.71979166667

亮星列表

1
2
3
rigel = ephem.star('Rigel')
print('%s %s' % (rigel._ra, rigel._dec))
5:14:32.30 -8:12:06.0

城市列表

1
2
3
4
5
6
7
8
9
10
11
12
stuttgart = ephem.city('Stuttgart')
print(stuttgart)
print(stuttgart.lon)
<ephem.Observer date='2018/8/27 05:17:01' 
epoch='2000/1/1 12:00:00' 
lon='9:10:50.8' 
lat='48:46:37.6' 
elevation=249.205185m 
horizon=0:00:00.0 
temp=15.0C 
pressure=983.6685758505818mBar>
9:10:50.8

天文常量

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
#### PyEphem provides constants for the dates of a few major star-atlas epochs:
#### 不同历元#  B1900# B1950# J2000
### PyEphem provides, for reference, the length of four distances, all in meters:

### 天文常量print(ephem.meters_per_au)
print(ephem.earth_radius)
print(ephem.moon_radius)
print(ephem.sun_radius)
#### PyEphem provides the speed of light in meters per second:
## 光速 m/sprint(ephem.c)
149597870000.0
6378160.0
1740000.0
695000000.0
299792458.0

《黑龙江常见野生植物图鉴》出版

《黑龙江常见野生植物图鉴》由李晶、苍晶、张金龙编著，本书以野外拍摄的精美照片展示了黑龙江省常见野生植物资源，对每种植物的形态特征、生境和用途进行了必要的描述。

全书共收录黑龙江省常见的种子植物449种(83科，276属)，学名以全球植物名录数据库中的最新修订为准，植物编排参照恩格勒系统；对被子植物APGⅢ分类系统进行了简单介绍，并附每种植物在APGⅢ系统中的分类概况。

为了方便非专业的植物分类学爱好者阅读，大部分植物附加了当地广泛接受的俗名，并在书后附有常用植物形态学术语图示和常见科植物的野外识别要点。

本书配套了数字课程，内容包括与书中记录物种相关的大量照片等参考资料，并会不断补充和更新。本书可作为高等院校植物学野外实习用书，也是相关学科教师、学生、科研工作者查阅和鉴定黑龙江省常见野生植物的参考书，以及植物分类学爱好者的实用工具书。

http://product.dangdang.com/1050991625.html

一个简单的线性插值问题

举例： 一个均匀的坐标轴上，x取0.5时，y为0.7， x取11.5时，y为13.9，问x取6.1时， y的值。

该题为插值问题，可以用线性插值公式解决：
参见： http://zh.wikipedia.org/wiki/%E7%BA%BF%E6%80%A7%E6%8F%92%E5%80%BC
该公式其实也可以通过解一元二次方程来推导。

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
## 线性方程一般的形式如下：
y = a*X + b

所以题目中的关系也可用以下方程组来表示： 
## 0.7  = a*0.5 + b
## 13.9 = a*11.5 + b

即： 
a*0.5 + b = 0.7
a*11.5 + b ＝ 13.9

##写成矩阵形式如下
Ax =b

###在R中，可用solve 函数对线性方程组求解
### 示例R code 如下：
A <- matrix(c(0.5, 1, 11.5, 1), nrow = 2, byrow = TRUE)
b <- c(0.7, 13.9)
res <- solve(A, b)

### 写成函数
fx <- function(x){   
   y <- res[1]*x + res[2]
   return(y)
}

fx(6.1)
### 结果为7.42
### 所以当x取值为6.1时， y为7.42

R packages for Phylogenetic Comparative Methods

adephylo:
exploratory analyses for the phylogenetic comparative method

ape:
Analyses of Phylogenetics and Evolution

apTreeshape:
Analyses of Phylogenetic Treeshape

BAMMtools:
Analysis and visualization of macroevolutionary dynamics on phylogenetic trees

caper:
Comparative Analyses of Phylogenetics and Evolution in R

cladoRcpp:
C++ implementations of phylogenetic calculations

corHMM:
Analysis of binary character evolution

DDD:
Diversity-dependent diversification

distory:
Distance Between Phylogenetic Histories

diversitree:
comparative phylogenetic analyses of diversification

evobiR:
evolutionary biology in R

expoTree:
Calculate density dependent likelihood of a phylogenetic tree

geomorph:
Geometric morphometric analyses of 2d/3d landmark data

GUniFrac:
Generalized UniFrac distances

HMPTrees:
Statistical Object Oriented Data Analysis of RDP-based Taxonomic trees from Human Microbiome Data
Modeling, Visualization, and Two-Group Comparison

HyPhy:
Macroevolutionary phylogentic analysis of species trees and gene trees

iteRates:
Parametric rate comparison

laser:
Likelihood Analysis of Speciation/Extinction Rates from Phylogenies

MPSEM:
Modeling Phylogenetic Signals using Eigenvector Maps

mvSLOUCH:
MultiVariate Stochastic Linear Ornstein-Uhlenbeck models for phylogenetic Comparative hypotHeses

ouch:
Ornstein-Uhlenbeck models for phylogenetic comparative hypotheses

paleotree:
Paleontological and Phylogenetic Analyses of Evolution

pastis:
Phylogenetic Assembly with Soft Taxonomic Inferences

PCPS:
PCPS - Principal Coordinates of Phylogenetic Structure

pGLS:
Generalized Least Square in comparative Phylogenetics

phangorn:
Phylogenetic analysis in R

phyclust:
Phylogenetic Clustering (Phyloclustering)

phyext:
An extension of some of the classes in phylobase. Tree objects now support subnodes on branches

phylobase:
Base package for phylogenetic structures and comparative data

phyloclim:
Integrating Phylogenetics and Climatic Niche Modeling

phyloland:
Modelling Competitive Exclusion and Limited Dispersal in a Statistical Phylogeographic Framework

phylolm:
Phylogenetic Linear Regression

phylotools:
Phylogenetic tools for Eco-phylogenetics

phyreg:
Implements the Phylogenetic Regression of Grafen (1989)

phytools:
Phylogenetic Tools for comparative biology (and other things)

picante:
R tools for integrating phylogenies and ecology

PVR:
Computes phylogenetic eigenvectors regression (PVR) and phylogenetic signal-representation curve (PSR) (with null and Brownian expectations)

RADami:
R Package for Phylogenetic Analysis of RADseq Data

SigTree:
Determine significantly responsive branches in phylogenetic trees

spacodiR:
Spatial and Phylogenetic Analysis of Community Diversity

spider:
Species Identity and Evolution in R

SYNCSA:
Analysis of functional and phylogenetic patterns in metacommunities

TESS:
Fast simulation of reconstructed phylogenetic trees under time-dependent birth-death processes

treebase:
An R package for discovery, access and manipulation of online phylogenies

TreePar:
Estimating birth and death rates based on phylogenies

TreeSimGM:
Simulating Phylogenetic Trees under a General Model