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

Candice
Candice
Candice
96
文章
0
评论
2020-04-1709:05:00 评论 2,333 1987字
摘要

特定形式的非线性模型可以通过广义线性模型 (glm()) 拟合。但是许多时候,我们必须把非线性拟合的问题 作为一个非线性优化的问题解决。 R的非线性优化程序是 optim() 和 nlm()。

我们通过搜寻 参数值使得缺乏度(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.

来源:数据分析网

  • 我的微信公众号
  • 微信扫一扫
  • weinxin
  • 我的微信公众号
  • 微信扫一扫
  • weinxin
匿名

发表评论

匿名网友 填写信息

:?: :razz: :sad: :evil: :!: :smile: :oops: :grin: :eek: :shock: :???: :cool: :lol: :mad: :twisted: :roll: :wink: :idea: :arrow: :neutral: :cry: :mrgreen: