R语言中的最大似然估计

Question

R语言中的最大似然估计

4

我是一个对R和统计学都很陌生的新手。我正在尝试使用最大似然估计方法，并且得到了一些错误的结果。我想要用一个简单的线性函数来建模x：

x<-apply(matrix(seq(1,10,1), nrow=1), 1, function(x) 10*x+runif(10,-3,3))
LL<-function(a,b){
    R=apply(x,1,function(y) a*y+b)
    -sum(log(R))
    }
mle(LL, start=list(a=10, b=0))

我得到了以下结果:

Coefficients:
    a         b 
43571.957  1338.345

将a替换为~10，b替换为~0。

我根据Spacedman的建议修改了代码:

set.seed(99)
x<-apply(matrix(seq(1,10,1), nrow=1), 1, function(x) 10*x+runif(10,-3,3))
LL<-function(a,b){
R = x[,1] - a*(1:10) + b
-sum(R^2)
}
library(stats4)
mle(LL, start=list(a=11, b=0.3))

Error in solve.default(oout$hessian) : 
Lapack routine dgesv: system is exactly singular: U[1,1] = 0

我不知道如何解决这个错误。更改查找并重新生成x值也无济于事。

- robert

2

mle函数来自哪个包？在生成数据之前，请使用set.seed(99)，以便我们都使用相同的随机数。 - Spacedman

@Spacedman mle 是来自 stats4 包。 - robert

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- J.R. · Accepted Answer

这里有几件事情需要注意。为了澄清，我们首先将误差项的分布从均匀分布runif(x, -3, 3)更改为标准正态分布：rnorm(x)。现在我们可以轻松地模拟您的数据，然后设置您的（负）对数似然并通过以下方式最大化（最小化）：

a <- 10 
b <- 0
set.seed(99)
x <- apply(matrix(seq(1, 10, 1), nrow=1), 1, function(x) b + a * x + rnorm(10))
minuslogL <- function(a, b) -sum(dnorm(x[, 1] - (b + a * 1:10), log = TRUE))
library(stats4)
mle(minuslogL, start = list(a = 11, b = 0.3))

Call:
mle(minuslogl = minuslogL, start = list(a = 11, b = 0.3))

Coefficients:
        a         b 
9.8732793 0.5922192

请注意，这很有效，因为可能性是光滑的，mle()在优化时使用了“BFGS”，例如拟牛顿梯度方法。让我们尝试一下相同的内容，但加入均匀误差：

set.seed(99)
x <- apply(matrix(seq(1, 10, 1), nrow=1), 1, function(x) b + a * x + runif(10, -3, 3))
minuslogL2 <- function(a,b) -sum(dunif(x[, 1] -(a * 1:10 + b), -3, 3, log = TRUE))
mle(minuslogL2, start = list(a = 11, b = 0.3))

Error in optim(start, f, method = method, hessian = TRUE, ...) : 
  initial value in 'vmmin' is not finite

这个失败了！为什么？因为均匀误差限制了参数空间，你将得不到平滑的似然函数。如果你把参数a、b移动得离真实值太远，就会出现无穷大的情况。如果你移动足够接近，你将得到相同的似然函数（例如：许多可能的最小值）：

> minuslogL2(11, 0.3)
[1] Inf
> minuslogL2(10, 0)
[1] 17.91759
> minuslogL2(10.02, 0.06)
[1] 17.91759

最大化这个可能性相当于找到集合：{a，b}：-logL（a，b）== -logL（10,0），它可以通过普通的搜索算法找到。