用R语言进行数据分析：最小二乘法和最大似然法

我们通过搜寻 参数值使得缺乏度（lack-of-fit）最低，如 nlm() 就是通过循环调试各种参数值得到最优值。 和线性回归不同，程序不一定会 收敛到一个稳定值。 nlm() 需要设定参数搜索的起始值， 而参数估计是否收敛在很大程度上依赖于 起始值设置的质量。

一、最小二乘法

拟合非线性模型的一种办法就是使误差平方和（SSE）或残差平方和最小。 如果观测到的误差极似正态分布， 这种方法是非常有效的。

下面是例子来自 Bates & Watts (1988)，51页。具体数据是：

> x <- c(0.02, 0.02, 0.06, 0.06, 0.11, 0.11, 0.22, 0.22, 0.56, 0.56,

1.10, 1.10)

> y <- c(76, 47, 97, 107, 123, 139, 159, 152, 191, 201, 207, 200)

被拟合的模型是：

> fn <- function(p) sum((y - (p[1] * x)/(p[2] + x))^2)

为了进行拟合，我们需要估计参数拟合起始值。一种判断 合理起始值的办法把数据图形化，然后估计一些参数值， 并且利用这些值初步绘制模型曲线。

> plot(x, y)

> xfit <- seq(.02, 1.1, .05)

> yfit <- 200 * xfit/(0.1 + xfit)

> lines(spline(xfit, yfit))

当然，我们可以做的更好，但是初始值 200 和 0.1 已经足够了。 现在做拟合：

> out <- nlm(fn, p = c(200, 0.1), hessian = TRUE)

拟合后，out$minimum 是误差的平方和（SSE）， out$estimate 是参数的最小二乘估计值。 为了得到参数估计过程中近似的标准误(SE)，我们可以：

> sqrt(diag(2*out$minimum/(length(y) - 2) * solve(out$hessian)))

上述命令中的2表示参数的个数。一个95% 的信度区间可以通过 +/- 1.96 SE 计算得到。我们可以把这个最小二乘拟合曲线画在一个图上：

> plot(x, y)

> xfit <- seq(.02, 1.1, .05)

> yfit <- 212.68384222 * xfit/(0.06412146 + xfit)

> lines(spline(xfit, yfit))

标准包 stats 提供了许多用最小二乘法 拟合非线性模型的扩充工具。我们刚刚拟合过的模型是 Michaelis-Menten 模型，因此可以利用下面的命令得到类似的结论。

> df <- data.frame(x=x, y=y)

> fit <- nls(y ~ SSmicmen(x, Vm, K), df)

> fit

Nonlinear regression model

model: y ~ SSmicmen(x, Vm, K)

data: df

Vm K

212.68370711 0.06412123

residual sum-of-squares: 1195.449

> summary(fit)

Formula: y ~ SSmicmen(x, Vm, K)

Parameters:

Estimate Std. Error t value Pr(>|t|)

Vm 2.127e+02 6.947e+00 30.615 3.24e-11

K 6.412e-02 8.281e-03 7.743 1.57e-05

Residual standard error: 10.93 on 10 degrees of freedom

Correlation of Parameter Estimates:

Vm

K 0.7651

二、最大似然法

最大似然法（Maximum likelihood）也是一种非线性拟合办法。它甚至可以用于 误差非正态的数据。这种方法估计的参数 将会使得对数似然值最大或者负的对数似然值 最小。下面的例子来自 Dobson (1990), pp. 108–111。这个例子对剂量－响应数据拟合 logistic模型 （当然也可以用 glm() 拟合）。数据是：

> x <- c(1.6907, 1.7242, 1.7552, 1.7842, 1.8113,

1.8369, 1.8610, 1.8839)

> y <- c( 6, 13, 18, 28, 52, 53, 61, 60)

> n <- c(59, 60, 62, 56, 63, 59, 62, 60)

要使负对数似然值最小，则：

> fn <- function(p)

sum( - (y*(p[1]+p[2]*x) - n*log(1+exp(p[1]+p[2]*x))

+ log(choose(n, y)) ))

我们选择一个适当的起始值，开始拟合：

> out <- nlm(fn, p = c(-50,20), hessian = TRUE)

拟合后，out$minimum 就是负对数似然值， out$estimate 就是最大似然拟合的参数值。 为了得到拟合过程近似的标准误，我们可以：

> sqrt(diag(solve(out$hessian)))

一个 95% 信度期间可根据公式 +/- 1.96 SE 计算得到。

End.

来源：数据分析网