没有什么问题optimx
。退一步看看你想要最大化的功能:log(p1) + log(p2) + lamda*(p1+p2-1)
。很直观的是,最佳解决方案是使所有变量尽可能大,不是吗?看到那个optimx
正确返回您指定的上限。
那么你的方法有什么问题呢?使用拉格朗日乘子时,临界点是上面函数的鞍点,而不是像局部最小值optimx
会帮助你找到。因此,您需要修改问题,使这些鞍点成为局部最小值。这可以通过优化梯度范数来完成,这很容易针对您的问题进行分析计算。这里有一个很好的例子(带图片):
http://en.wikipedia.org/wiki/Lagrange_multiplier#Example:_numerical_optimization http://en.wikipedia.org/wiki/Lagrange_multiplier#Example:_numerical_optimization.
对于您的问题:
grad.norm <- function(x) {
lambda <- tail(x, 1)
p <- head(x, -1)
h2 <- sum((1/p + lambda)^2) + (sum(p) - 1)^2
}
optimx(c(.48, .52, -1.5),
grad.norm,
lower = c(rep(.1, 2), -3),
upper = c(rep(.9, 2), -1),
method = "L-BFGS-B")
# par fvalues method fns grs [...]
# 1 0.5000161, 0.5000161, -1.9999356 1.038786e-09 L-BFGS-B 13 13 [...]
跟进:如果您不想或无法自己计算梯度,您可以让 R 计算数值近似值,例如:
log.like <- function(x) {
lambda <- tail(x, 1)
p <- head(x, -1)
return(sum(log(p)) + lambda*(sum(p) - 1))
}
grad.norm <- function(x) {
require(numDeriv)
return(sum(grad(log.like, x)^2))
}
optimx(c(.48, .52, -1.5),
grad.norm,
lower = c(rep(.1, 2), -3),
upper = c(rep(.9, 2), -1),
method = "L-BFGS-B")
# par fvalues method fns grs [...]
# 1 0.5000161, 0.5000161, -1.9999356 1.038784e-09 L-BFGS-B 13 13 [...]