非線形関数に対する最小値問題の解法として，Levenberg-Marquerdt法は非常に良く見かける方法だと思う。概念としては解から遠いところでは再急降下法を使い，解に近づいたらガウスーニュートン法に切り替えて，繰り返し計算により目的関数を小さくしていく方法だ。繰り返し1回毎に最急降下法とガウス―ニュートン法のどちらを使うかという度合いは，重み関数λを導入し，それを適宜更新することで切り替えていく。
確かにLevenberg-Marquerdt法は優秀だが，Hessianを計算しなくてはならず，俺のような面倒くさがりにとっては身構えてしまう面がある。ちょっと試しにコードを書いてみるにしても正直言って面倒くさいのだ。さらにもっと言えばλの導入によりガウス-ニュートン法よりなお複雑とも言える。また各パラメータによる勾配を計算するのにも計算間違いの危険性を孕んでいると，なかなかに心理的ハードルが高い。
こんな時，Rがインストールされていると心強い。パッケージ"minpack.lm"に入っている"nls.lm"を使えばさしたる面倒もなくLevenberg-Marquerdt法を使える。関数を設定し，パラメータのlist，誤差関数，引数に代入するデータ列などを指定するだけで繰り返し計算を行ってくれるのである。nls.lmが優れているのは，各パラメータによる関数の微分表式を打ち込む必要が無いことだと思う。
Let's try
自分の備忘録的な意味も込め，例を出して使い方をメモしておく。ノイズが乗った測定データをうまくフィッティングして，パラメータの真値を知りたいような状況としよう。Rでminpack.lmを読みこんだ上で例えば次のように打ち込む。
# まず，データ点のx座標を準備
x <- seq(-1.0, 1.0, length = 30)
# 真のモデルの表式はLorentzian
getPred <- function(parS, xx) parS$a * 1/(1 + (parS$b * xx^2))
# パラメータの真の値はこんな感じ
pp <- list(a=1.0, b=2.0)
おなじみのろーれんちあんが現れる。
実際に測定されるデータを再現するために標準偏差0.01で正規分布するノイズを与えてデータ列をsimDに格納してみよう。
simD <- getPred(pp, x) + rnorm(length(x), sd = 0.1)
測定データとしてはかなりデタラメな感じになってしまった。本当なら山のテッペンあたりで密にデータを取るべきだが，まあここでは気にせず進めよう。
# 残差関数
residFun <- function(parS, observed, xx) observed - getPred(parS, xx)
# 繰り返し計算の初期値。適当に与える。
parStart <- list(a = 0.5, b = 1.0)
# 呪文実行
nls.out <- nls.lm(par = parStart, fn = residFun, observed = simD, xx = x, control = nls.lm.control(nprint = 1))
基本的にはこれだけで操作は完了。フィッティングされたパラメータをもとに曲線を描き加えてみる。
だいたい良い感じである。
# サマリーを表示
summary(nls.out)
すると
Parameters:
Estimate Std. Error t value Pr(>|t|)
a 0.99331 0.01685 58.94 < 2e-16 ***
b 2.04580 0.12535 16.32 7.77e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.04465 on 28 degrees of freedom
Number of iterations to termination: 5
Reason for termination: Relative error in the sum of squares is at most `ftol'.
とまあ，めでたく簡潔に結果が表示される。大量の試行データをフィットするときには工夫がいるかもしれないが，とりあえずちょっとフィットしてみるときには便利なツールだろう。
