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

要了解响应变量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)