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

张金龙

jinlongzhang01@gmail.com

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

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)