r - R の最適化の問題

Question

Rの最適関数に関する質問

これまでのところ、次のコードがあり、X と T の値を入力することを知る必要があります。X は 10 個の値のベクトルであり、T は平均と分散に関連する 10*2 値のベクトルです。出力を、alpha、mean1、mean2、var1、および var2 の 1 つの新しい値の形式にしたいと考えています。入力データを正しく取得する方法がわかりません。

この関数で X のすべての値を実行したいのですが、T の最初の行 (10 個の値) だけを実行したいのですが、T に対してこれを行う方法がわかりません。2 番目の行には別の関数があります。

R <-runif(10, 0, 1)
S <-1-R
T <-t(cbind(R,S))


X <- runif(10, 25, 35)


Data1 <- function(xy) { 
  alpha <- xy[1]
  mean1 <- xy[2]
  mean2 <- xy[3]
  var1 <- xy[4]
  var2 <- xy[5] 

  -sum(0.5*(((X)-mean1)/var1)^2+alpha*mean1+log(2.5*var1)+log(exp(-alpha*mean1)+exp(-alpha*mean2))*(T))
}
starting_values <- c(0.3, 28, 38, 4, 3)
optim(starting_values, Data1, lower=c(0, 0, 0, 0, 0), method='L-BFGS-B')

次のエラーコードも取得します。

Error in optim(starting_values, Data1, lower = c(0, 0, 0, 0, 0), method = "L-BFGS-B") : 
  L-BFGS-B needs finite values of 'fn'

どんな助けにも乾杯。

編集

含めるための 2 番目の関数

0.5*((y1-mean2)/var2)^2+alpha*mean2+log(2.5*var2)+ log(exp(-alpha*mean1)+exp(-alpha*mean2)))*T

わかりましたので、私がやりたいことをできるだけ明確に説明してください。上記の元の投稿の最初の関数は、X の 10 個の値すべてを一度に 1 つずつ取得し、T データの最初の行 (ここでは R とラベル付けされています) を取得する必要がありますが、これを行う方法がわかりません。

上記の 2 番目の関数は、X の 10 個の値すべてを連続して取得し、次に T から 2 番目のデータ行を取得する必要があります (以下では S とラベル付けされています)。

これらすべてが合計されます。このようにして、5 つの未知のパラメータが推定されます。

T

       [,1]      [,2]      [,3]       [,4]      [,5]        [,6]      [,7]      [,8]      [,9]     [,10]
R 0.1477715 0.3055021 0.2963543 0.04149945 0.8342484 0.996865333 0.1592568 0.4623762 0.8000778 0.6979342
S 0.8522285 0.6944979 0.7036457 0.95850055 0.1657516 0.003134667 0.8407432 0.5376238 0.1999222 0.3020658

編集2

同じシードを実行しても、ベンと同じ値が得られません。すべてのパッケージがインストールされていることを確認しましたが、インストールされているように見えます。同じ最終的な答えが得られず、opt2$par の個々のアイテムを呼び出すこともできません。一連の出力を提供する代わりに、最初の数行と最後の数行を提供します。

0.3 28 38 4 3 -74.97014 -120.7212 
Loading required package: BB
Loading required package: quadprog
Loading required package: ucminf
Loading required package: Rcgmin
Loading required package: Rvmmin

Attaching package: ‘Rvmmin’

The following object(s) are masked from ‘package:optimx’:

    optansout

Loading required package: minqa
Loading required package: Rcpp
0.3 28 38 4 3 -74.97014 -120.7212 
0.9501 28 38 4 3 -176.3368 -265.9074 
1.9001 28 38 4 3 -324.7782 -478.4652 
0.9501 28.95 38 4 3 -179.9994 -260.8711 
0.9501 28 38.95 4 3 -176.3366 -283.0445 
0.9501 28 38 4.95 3 -176.7836 -265.9074 
0.9501 28 38 4 3.95 -176.3368 -254.6188

...................

16.32409 27.86113 38.54337 3.940143 2.504167 -2566.194 -3826.233 
16.32409 27.86113 38.54337 3.940044 2.504167 -2566.194 -3826.233 
16.32409 27.86113 38.54337 3.940093 2.504199 -2566.194 -3826.232 
16.32409 27.86113 38.54337 3.940093 2.504136 -2566.194 -3826.234 
> opt2$par
$par
[1] 16.324085 27.861134 38.543373  3.940093  2.504167

> opt2$par["mean1"]
$<NA>
NULL

score 7 · Accepted Answer

最初のクラック: 上記のコードを使用しました。set.seed(101)再現性のために最初に追加しました。

明確にするために関数をわずかに再フォーマットしましたが、重要なことは何も変更せず、cat()デバッグ目的でステートメントを追加しました。

Data1 <- function(xy) {
    alpha <- xy[1]; mean1 <- xy[2]; mean2 <- xy[3]
    var1 <- xy[4]; var2 <- xy[5]

    r1 <- -sum(0.5*((X-mean1)/var1)^2+
           alpha*mean1+
           log(2.5*var1)+
           log(exp(-alpha*mean1)+
               exp(-alpha*mean2))*T[1,])
    r2 <- -sum(0.5*((X-mean2)/var2)^2+
           alpha*mean2+
           log(2.5*var2)+
           log(exp(-alpha*mean1)+exp(-alpha*mean2))*T[2,])

    cat(xy,r1,r2,"\n")
   r1+r2
}

わずかに圧縮されたバージョンで、(1)with関数をよりクリーンにするために利用します。(2) R の複製およびベクター再利用機能を使用する

Data2 <- function(xy) {
    with(as.list(xy),
    {
       mmat <- rep(c(mean1,mean2),each=ncol(T))
       vmat <- rep(c(var1,var2),each=ncol(T))
       -sum(0.5*((X-mmat)/vmat)^2+
          alpha*mmat+
          log(2.5*vmat)+
          log(exp(-alpha*mean1)+exp(-alpha*mean2))*T)
    })
}

ここで、開始値の名前付きベクトルが必要です。

 starting_values <- c(alpha=0.3, mean1=28, mean2=38, var1=4, var2=3)

結果が一致することを確認します。

 Data1(starting_values) ##  [1] -195.6913
 Data2(starting_values) ##  [1] -195.6913

これは失敗します (ただし、どのように失敗するかについての情報が得られます):

 optim(par=starting_values, Data1, lower=rep(1e-4,5), method='L-BFGS-B',
     control=list(trace=6))

多くの出力が生成され、次のようになります。

##  21.29998 27.97361 37.98915 4.011199 3.001 -6014.225 
## 21.29998 27.97361 37.98915 4.011199 2.999 -6014.225 
## 85.29991 27.89318 37.95606 4.04533 3 Inf 
## Error in optim(par = starting_values, Data1, lower = rep(1e-04, 5), 
##    method = "L-BFGS-B",  : 
##     L-BFGS-B needs finite values of 'fn'

これにより、少なくともどこで問題が発生したかがわかります。ここで、エクスプレッションを 1 つずつ評価して、どのビットがオーバーフローしたかを確認します。

チャットルームのコメンター（ジャスティン）が言ったように、

alpha、mean1、および mean2 は無制限であるため、第 3 項 log(exp(...) + exp(...)) は非常に迅速に -Inf になります。exp(-大きい数 * 大きい数) ~ 0

さらにデバッグするには、次のことができます。

オーバーフローを避けるために、関数の評価を再配置してみてください
オーバーフローを避けるために、一部のパラメータに上限を設定します
関数 test があり、適切な場合に Inf ではなく非常に大きな値を返す

残念ながら、L-BFGS-B他のいくつかのオプティマイザーよりも壊れやすく、非有限値を許可しません。

次に、境界を許可し、非有限値を処理するパッケージのオプティマイザーを試しました (また、導関数を使用しない方法であり、一般に導関数ベースの方法よりもわずかに遅くなる傾向がありますが、より堅牢です):動作するbobyqaようですわかりましたが、答えが賢明かどうかはわかりません。optimx

library(optimx)
opt2 <- optimx(par=starting_values, 
      Data1, lower=rep(1e-4,5), method='bobyqa')
opt3 <- optimx(par=starting_values, 
      Data2, lower=rep(1e-4,5), method='bobyqa')

問題ないように見えます（これが賢明な答えである場合、私にはわかりません）。

> opt2$par
$par
    alpha     mean1     mean2      var1      var2 
16.330752 27.815324 38.497483  3.894179  2.447219 

> opt3$par
$par
    alpha     mean1     mean2      var1      var2 
16.330900 27.820813 38.491290  3.887975  2.456052

答えがわずかに異なることに注意してください (var2 の場合は約 0.5%)。これは、適合がわずかに不安定である/表面が非常に平らである可能性があることを示唆しています。(そしてData2は同じData1答えを与えるはずであり、開始値についてもそうしますが、操作の順序により、一部の入力に対して非常にわずかに異なる答えが得られると思います-またはどこかで失敗しました...)

この近似から個々のコンポーネントを抽出するには、たとえばmean1、ベクトルインデックスを使用します。

opt3$par["mean1"]  ## 27.820813

r - R の最適化の問題

1 に答える 1

Related

Reference