R中的最大可能性

Question

我对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 

代替〜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值不会有帮助。

Answer 1

这里有几件事要注意。为了澄清，我们首先将错误项的分布从统一分布runif(x, -3, 3)更改为std。正态分布：rnorm(x)。现在，我们可以很容易地模仿你的数据，然后设置（减）数似然最大化（minizime）由：

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太远离真正的值，你会得到Inf。如果你靠拢的话，你会得到相同的可能性（例如许多可能的最小值。）：

> 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)，它可以通过一个简单的搜索算法找到。