9

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=2isp1=p2=1/2と。のときの方程式の解lamda=-2。ただし、optimxを使用すると、これは取得されません。何か案が？

4

1 に答える 1

23

何も問題はありません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 [...]
于 2012-08-09T05:12:25.767 に答える