R中均匀分布的最大似然估计导致荒谬的结果
Maximum likelihood estimation with uniform distribution in R leads to absurd result
我想使用 mle
函数获取 Unif(a,b)
分布中 a
和 b
的估计值。但我得到的估计值远不及 1 和 3。
library(stats4)
set.seed(20161208)
N <- 100
c <- runif(N, 1, 3)
LL <- function(min, max) {
R <- runif(100, min, max)
suppressWarnings((-sum(log(R))))
}
mle(minuslogl = LL, start = list(min = 1, max = 3), method = "BFGS",
lower = c(-Inf, 0), upper = c(Inf, Inf))
我得到了:
Call:
mle(minuslogl = LL, start = list(min = 1, max = 3), method = "BFGS")
Coefficients:
min max
150.8114 503.6586
知道发生了什么事吗?提前致谢!
我先指出你的代码哪里错了。
你需要 dunif
而不是 runif
。您可以定义:
LL <- function (a, b) -sum(dunif(x, a, b, log.p = TRUE))
我下面的代码中没有使用dunif
,因为密度只是1 / (b - a)
所以我直接写了。
- 您正在 objective 函数内生成样本。 对于
U[a,b]
这没关系,因为它的密度没有 x
。但对于其他分布,objective 函数在每次迭代时都会发生变化。
- 有了框约束,你需要
method = "L-BFGS-B"
,而不是普通的 "BFGS"
。而且您没有使用正确的约束。
现在更深入...
对于来自 U[a, b]
的长度 n
样本向量 x
,似然为 (b - a) ^ (-n)
,负对数似然为 n * log(b - a)
。显然 MLE 是 a = min(x)
和 b = max(x)
.
数值优化是完全没有必要的,实际上没有约束是不可能的。看梯度向量:
( n / (a - b), n / (b - a) )
偏导数w.r.t。 a
/ b
始终为负数/正数且不能为 0。
当我们施加框约束时,数值方法变得可行:-Inf < a <= min(x)
和 max(x) <= b < Inf
。我们肯定知道迭代在边界处终止。
我下面的代码同时使用了 optim
和 mle
。注意 mle
在反转 Hessian 矩阵时会失败,因为它是奇异的:
-(b - a) ^ 2 (b - a) ^ 2
(b - a) ^ 2 -(b - a) ^ 2
代码:
## 100 samples
set.seed(20161208); x <- runif(100, 1, 3)
# range(x)
# [1] 1.026776 2.984544
## using `optim`
nll <- function (par) log(par[2] - par[1]) ## objective function
gr_nll <- function (par) c(-1, 1) / diff(par) ## gradient function
optim(par = c(0,4), fn = nll, gr = gr_nll, method = "L-BFGS-B",
lower = c(-Inf, max(x)), upper = c(min(x), Inf), hessian = TRUE)
#$par
#[1] 1.026776 2.984544 ## <- reaches boundary!
#
# ...
#
#$hessian ## <- indeed singular!!
# [,1] [,2]
#[1,] -0.2609022 0.2609022
#[2,] 0.2609022 -0.2609022
## using `stats4::mle`
library(stats4)
nll. <- function (a, b) log(b - a)
mle(minuslogl = nll., start = list(a = 0, b = 4), method = "L-BFGS-B",
lower = c(-Inf, max(x)), upper = c(min(x), Inf))
#Error in solve.default(oout$hessian) :
# Lapack routine dgesv: system is exactly singular: U[2,2] = 0
我想使用 mle
函数获取 Unif(a,b)
分布中 a
和 b
的估计值。但我得到的估计值远不及 1 和 3。
library(stats4)
set.seed(20161208)
N <- 100
c <- runif(N, 1, 3)
LL <- function(min, max) {
R <- runif(100, min, max)
suppressWarnings((-sum(log(R))))
}
mle(minuslogl = LL, start = list(min = 1, max = 3), method = "BFGS",
lower = c(-Inf, 0), upper = c(Inf, Inf))
我得到了:
Call:
mle(minuslogl = LL, start = list(min = 1, max = 3), method = "BFGS")
Coefficients:
min max
150.8114 503.6586
知道发生了什么事吗?提前致谢!
我先指出你的代码哪里错了。
你需要
dunif
而不是runif
。您可以定义:LL <- function (a, b) -sum(dunif(x, a, b, log.p = TRUE))
我下面的代码中没有使用
dunif
,因为密度只是1 / (b - a)
所以我直接写了。- 您正在 objective 函数内生成样本。 对于
U[a,b]
这没关系,因为它的密度没有x
。但对于其他分布,objective 函数在每次迭代时都会发生变化。 - 有了框约束,你需要
method = "L-BFGS-B"
,而不是普通的"BFGS"
。而且您没有使用正确的约束。
现在更深入...
对于来自 U[a, b]
的长度 n
样本向量 x
,似然为 (b - a) ^ (-n)
,负对数似然为 n * log(b - a)
。显然 MLE 是 a = min(x)
和 b = max(x)
.
数值优化是完全没有必要的,实际上没有约束是不可能的。看梯度向量:
( n / (a - b), n / (b - a) )
偏导数w.r.t。 a
/ b
始终为负数/正数且不能为 0。
当我们施加框约束时,数值方法变得可行:-Inf < a <= min(x)
和 max(x) <= b < Inf
。我们肯定知道迭代在边界处终止。
我下面的代码同时使用了 optim
和 mle
。注意 mle
在反转 Hessian 矩阵时会失败,因为它是奇异的:
-(b - a) ^ 2 (b - a) ^ 2
(b - a) ^ 2 -(b - a) ^ 2
代码:
## 100 samples
set.seed(20161208); x <- runif(100, 1, 3)
# range(x)
# [1] 1.026776 2.984544
## using `optim`
nll <- function (par) log(par[2] - par[1]) ## objective function
gr_nll <- function (par) c(-1, 1) / diff(par) ## gradient function
optim(par = c(0,4), fn = nll, gr = gr_nll, method = "L-BFGS-B",
lower = c(-Inf, max(x)), upper = c(min(x), Inf), hessian = TRUE)
#$par
#[1] 1.026776 2.984544 ## <- reaches boundary!
#
# ...
#
#$hessian ## <- indeed singular!!
# [,1] [,2]
#[1,] -0.2609022 0.2609022
#[2,] 0.2609022 -0.2609022
## using `stats4::mle`
library(stats4)
nll. <- function (a, b) log(b - a)
mle(minuslogl = nll., start = list(a = 0, b = 4), method = "L-BFGS-B",
lower = c(-Inf, max(x)), upper = c(min(x), Inf))
#Error in solve.default(oout$hessian) :
# Lapack routine dgesv: system is exactly singular: U[2,2] = 0