r - 为什么R中的optimx无法为这种简单的非参数似然最大化提供正确的解决方案？

optimx()提供了错误的解决方案还是我遗漏了一个简单的观点？谢谢!

我正在尝试最大化一个非常简单的可能性。这是非参数可能性
从某种意义上说，F的分布不是通过参数指定的。相当，
对于每个观察到的xi，f(xi)=pi以及log(Likelihood)=Sum(log(f(xi)))=Sum(log(pi))。

我试图最大化的功能是:sum(log(pi))+lamda(sum(pi-1))其中sum(pi)=1(即这是一个受约束的最大化问题，可以使用拉格朗日乘数来解决)。

这个问题的答案是pi=1/n，其中n是数据点的数量。但是，optimx似乎没有提供这种解决方案。有人有什么主意吗？如果是n=2，我要最大化的函数是log(p1)+log(p2)+lamda(p1+p2-1)。

这是我的代码和R的输出:

n=2
log.like=function(p)
{
  lamda=p[n+1]
  ll=0
  for(i in 1:n){
    temp = log(p[i])+lamda*p[i]-lamda/(n)
    ll=ll+temp
  }
  return(-ll)
}


mle = optimx(c(0.48,.52,-1.5),
             log.like,
             lower=c(rep(0.1,2),-3),
             upper=c(rep(.9,2),-1),
             method = "L-BFGS-B")

> mle
             par  fvalues   method fns grs itns conv  KKT1 KKT2 xtimes
1 0.9, 0.9, -1.0 1.010721 L-BFGS-B   8   8 NULL    0 FALSE   NA      0

n=2为p1=p2=1/2和lamda=-2时方程的解。但是，使用optimx时我没有得到这个。任何想法？

最佳答案

optimx没问题。退后一步，查看您要最大化的函数:log(p1) + log(p2) + lamda*(p1+p2-1)。直觉上说，最佳解决方案是使所有变量尽可能大，不是吗？看到optimx正确返回了您指定的上限。

那么，您的方法有什么问题呢？使用拉格朗日乘数时，临界点是上面函数的鞍点，而不是像optimx这样的局部最小值不能帮助您找到。因此，您需要以使这些鞍点成为局部最小值的方式来修改您的问题。这可以通过优化梯度的范数来完成，这很容易针对您的问题进行分析计算。这里有一个很好的例子(带有图片):

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 [...]

关于r - 为什么R中的optimx无法为这种简单的非参数似然最大化提供正确的解决方案？，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/11876291/