私の問題を説明するために再現可能な例で更新してください
私の最初の質問は、「R での信頼領域反射最適化アルゴリズムの実装」でした。ただし、再現可能な例を作成する途中で（@Benのアドバイスに感謝します）、私の問題はMatlabの1つの関数であることに気付きましたlsqnonlin私が持っているほとんどの場合に十分です（つまり、良い開始値を選択する必要がないことを意味します）が、Rにはそのような1対1の関数はありません。さまざまな最適化アルゴリズムが、さまざまなケースで適切に機能します。アルゴリズムが異なれば、到達する解も異なります。この背後にある理由は、R の最適化アルゴリズムが Matlab の信頼領域反射アルゴリズムより劣っているということではなく、R が自動微分を処理する方法にも関連している可能性があります。この問題は、実際には 2 年前に中断された作業が原因です。当時、パッケージoptimxの作成者の 1 人である John C Nash 教授は、Matlab lsqnonlin が R の最適化関数/アルゴリズムよりも優れたパフォーマンスを発揮する理由かもしれません。
以下の例は、私が遭遇したいくつかの問題を示しています (より再現可能な例が追加されています)。サンプルを実行するには、まず を実行しますinstall_github("KineticEval","zhenglei-gao")。パッケージmkinとその依存関係をインストールする必要があり、さまざまな最適化アルゴリズム用に他のパッケージを多数インストールする必要がある場合もあります。
lsqnonlin基本的に、Matlab 関数のドキュメント ( http://www.mathworks.de/de/help/optim/ug/lsqnonlin.html )で説明されているように、非線形最小二乗曲線近似問題を解決しようとしています。私の場合の曲線は、一連の微分方程式によってモデル化されています。例を挙げてもう少し説明します。私が試した最適化アルゴリズムは次のとおりです。
nls.lmレーベンバーグ・マルカート出身のマルク- からのポート
nlm.inb - からのL-BGFS-B
optim - からのspg
optimx solnpパッケージのRsolnp
他にもいくつか試しましたが、ここには示していません。
私の質問の要約
lsqnonlin私のタイプの非線形最小二乗問題を解決できるMatlabのように、Rには信頼できる関数/アルゴリズムがありますか? （私はそれを見つけることができませんでした。）- 単純なケースで、異なる最適化が異なるソリューションに到達する理由は何ですか?
lsqnonlinR の機能の優れている点は何ですか? 信頼領域反射アルゴリズムまたはその他の理由?- 私のタイプの問題を別の方法で解決するためのより良い方法はありますか? 多分もっと簡単な解決策があるかもしれませんが、私はそれを知りません。
例 1: 単純なケース
最初に R コードを与え、後で説明します。
ex1 <- mkinmod.full(
Parent = list(type = "SFO", to = "Metab", sink = TRUE,
k = list(ini = 0.1,fixed = 0,lower = 0,upper = Inf),
M0 = list(ini = 195, fixed = 0,lower = 0,upper = Inf),
FF = list(ini = c(.1),fixed = c(0),lower = c(0),upper = c(1)),
time=c(0.0,2.8, 6.2, 12.0, 29.2, 66.8, 99.8, 127.5, 154.4, 229.9, 272.3, 288.1, 322.9),
residue = c( 157.3, 206.3, 181.4, 223.0, 163.2, 144.7,85.0, 76.5, 76.4, 51.5, 45.5, 47.3, 42.7),
weight = c( 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1)),
Metab = list(type = "SFO",
k = list(ini = 0.1,fixed = 0,lower = 0,upper = Inf),
M0 = list(ini = 0, fixed = 1,lower = 0,upper = Inf),
residue =c( 0.0, 0.0, 0.0, 1.6, 4.0, 12.3, 13.5, 12.7, 11.4, 11.6, 10.9, 9.5, 7.6),
weight = c( 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1))
)
ex1$diffs
Fit <- NULL
alglist <- c("L-BFGS-B","Marq", "Port","spg","solnp")
for(i in 1:5) {
Fit[[i]] <- mkinfit.full(ex1,plot = TRUE, quiet= TRUE,ctr = kingui.control(method = alglist[i],submethod = 'Port',maxIter = 100,tolerance = 1E-06, odesolver = 'lsoda'))
}
names(Fit) <- alglist
kinplot(Fit[[2]])
(lapply(Fit, function(x) x$par))
unlist(lapply(Fit, function(x) x$ssr))
最後の行からの出力は次のとおりです。
L-BFGS-B Marq Port spg solnp
5735.744 4714.500 5780.446 5728.361 4714.499
「Marq」と「solnp」を除いて、他のアルゴリズムは最適に達しませんでした。さらに、「spg」メソッド (「bobyqa」などの他のメソッドも) は、このような単純なケースに対してあまりにも多くの関数評価を必要とします。さらに、開始値を変更してk_Parent=0.0058、ランダムに選択された の代わりに (そのパラメーターの最適値) を作成すると0.1、「Marq」は最適値を見つけることができなくなります。(以下にコードを示します)。「solnp」が最適なものを見つけられないデータセットもありました。ただし、lsqnonlinMatlab で使用する場合、このような単純なケースでは問題は発生しませんでした。
ex1_a <- mkinmod.full(
Parent = list(type = "SFO", to = "Metab", sink = TRUE,
k = list(ini = 0.0058,fixed = 0,lower = 0,upper = Inf),
M0 = list(ini = 195, fixed = 0,lower = 0,upper = Inf),
FF = list(ini = c(.1),fixed = c(0),lower = c(0),upper = c(1)),
time=c(0.0,2.8, 6.2, 12.0, 29.2, 66.8, 99.8, 127.5, 154.4, 229.9, 272.3, 288.1, 322.9),
residue = c( 157.3, 206.3, 181.4, 223.0, 163.2, 144.7,85.0, 76.5, 76.4, 51.5, 45.5, 47.3, 42.7),
weight = c( 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1)),
Metab = list(type = "SFO",
k = list(ini = 0.1,fixed = 0,lower = 0,upper = Inf),
M0 = list(ini = 0, fixed = 1,lower = 0,upper = Inf),
residue =c( 0.0, 0.0, 0.0, 1.6, 4.0, 12.3, 13.5, 12.7, 11.4, 11.6, 10.9, 9.5, 7.6),
weight = c( 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1))
)
Fit_a <- NULL
alglist <- c("L-BFGS-B","Marq", "Port","spg","solnp")
for(i in 1:5) {
Fit_a[[i]] <- mkinfit.full(ex1_a,plot = TRUE, quiet= TRUE,ctr = kingui.control(method = alglist[i],submethod = 'Port',maxIter = 100,tolerance = 1E-06, odesolver = 'lsoda'))
}
names(Fit_a) <- alglist
lapply(Fit_a, function(x) x$par)
unlist(lapply(Fit_a, function(x) x$ssr))
最後の行からの出力は次のとおりです。
L-BFGS-B Marq Port spg solnp
5653.132 4866.961 5653.070 5635.372 4714.499
ここで最適化していることを説明します。上記のスクリプトを実行して曲線を確認した場合は、一次反応を伴う 2 コンパートメント モデルを使用して曲線を記述します。モデルを表す微分方程式は次のex1$diffsとおりです。
Parent
"d_Parent = - k_Parent * Parent"
Metab
"d_Metab = - k_Metab * Metab + k_Parent * f_Parent_to_Metab * Parent"
この単純なケースでは、微分方程式から 2 つの曲線を表す方程式を導き出すことができます。最適化されるパラメータは$M_0,k_p, k_m, c=\mbox{FF_parent_to_Met} $制約付き$M_0>0,k_p>0, k_m>0, 1> c >0$です。
$$
\begin{split}
y_{1j}&= M_0e^{-k_pt_i}+\epsilon_{1j}\\
y_{2j} &= cM_0k_p\frac{e^{-k_mt_i}-e^{-k_pt_i}}{k_p-k_m}+\epsilon_{2j}
\end{split}
$$
したがって、微分方程式を解かずに曲線を当てはめることができます。
BCS1.l <- mkin_wide_to_long(BCS1)
BCS1.l <- na.omit(BCS1.l)
indi <- c(rep(1,sum(BCS1.l$name=='Parent')),rep(0,sum(BCS1.l$name=='Metab')))
sysequ.indi <- function(t,indi,M0,kp,km,C)
{
y <- indi*M0*exp(-kp*t)+(1-indi)*C*M0*kp/(kp-km)*(exp(-km*t)-exp(-kp*t));
y
}
M00 <- 100
kp0 <- 0.1
km0 <- 0.01
C0 <- 0.1
library(nlme)
result1 <- gnls(value ~ sysequ.indi(time,indi,M0,kp,km,C),data=BCS1.l,start=list(M0=M00,kp=kp0,km=km0,C=C0),control=gnlsControl())
#result3 <- gnls(value ~ sysequ.indi(time,indi,M0,kp,km,C),data=BCS1.l,start=list(M0=M00,kp=kp0,km=km0,C=C0),weights = varIdent(form=~1|name))
## Coefficients:
## M0 kp km C
## 1.946170e+02 5.800074e-03 8.404269e-03 2.208788e-01
このようにすると、経過時間はほぼゼロになり、最適に達します。ただし、この単純なケースが常にあるとは限りません。モデルは複雑になる可能性があり、微分方程式を解く必要があります。例 2 を参照
例 2、複雑なモデル
私はずっと前にこのデータセットに取り組んでいましたが、次のスクリプトを自分で実行し終える時間がありませんでした。(実行が完了するまで数時間かかる場合があります。)
data(BCS2)
ex2 <- mkinmod.full(Parent= list(type = "SFO",to = c( "Met1", "Met2","Met4", "Met5"),
k = list(ini = 0.1,fixed = 0,lower = 0,upper = Inf),
M0 = list(ini = 100,fixed = 0,lower = 0,upper = Inf),
FF = list(ini = c(.1,.1,.1,.1),fixed = c(0,0,0,0),lower = c(0,0,0,0),upper = c(1,1,1,1))),
Met1 = list(type = "SFO",to = c("Met3", "Met4")),
Met2 = list(type = "SFO",to = c("Met3")),
Met3 = list(type = "SFO" ),
Met4 = list(type = "SFO", to = c("Met5")),
Met5 = list(type = "SFO"),
data=BCS2)
ex2$diffs
Fit2 <- NULL
alglist <- c("L-BFGS-B","Marq", "Port","spg","solnp")
for(i in 1:5) {
Fit2[[i]] <- mkinfit.full(ex2,plot = TRUE, quiet= TRUE,ctr = kingui.control(method = alglist[i],submethod = 'Port',maxIter = 100,tolerance = 1E-06, odesolver = 'lsoda'))
}
names(Fit) <- alglist
(lapply(Fit, function(x) x$par))
unlist(lapply(Fit, function(x) x$ssr))
これは、次のような警告メッセージが表示される例です。
DLSODA- At T (=R1) and step size H (=R2), the
corrector convergence failed repeatedly
or with ABS(H) = HMIN
In above message, R =
[1] 0.000000e+00 2.289412e-09
元の質問
Matlab Optimization Toolbox ソルバーで使用される方法の多くは、信頼領域に基づいています。CRAN タスク ビュー ページによると、パッケージtrust、trustOptim、minqaのみ に、信頼領域ベースのメソッドが実装されています。ただし、勾配とヘシアンが必要ですtrust。minqaでは、私が探しているものではないようです。私の個人的な経験から、Matlab の信頼領域 Reflective アルゴリズムは、R で試したアルゴリズムと比較してパフォーマンスが優れていることが多いため、R でこのアルゴリズムの同様の実装を見つけようとしました。trustOptimbobyqa
ここで関連する質問をしました: R function to search for function
Matthew Plourde が提供した回答はlsqnonlin、Matlab で同じ関数名を持つ関数を提供しますが、信頼領域反射アルゴリズムが実装されていません。Matthew Plourde の回答は一般に、関数を探している R ユーザーにとって非常に役立つと思うので、古いものを編集し、ここで新しい質問をしました。
もう一度検索しましたが、運がありません。同様の matlab 関数を実装する関数/パッケージがまだいくつかあります。そうでない場合は、Matlab 関数を直接 R に翻訳して、自分の目的に使用することが許されるかどうか疑問に思います。