汎用最適化関数 nlm



最尤推定や最小自乗法等、ある目的関数の最適化に用いる汎用関数

  • RWikiの該当ページはこちら
  • Rの Htmlヘルプ"nlm{stats}"も参考のこと
  • ごくごく簡単に
    • 今、y=ax^2に近似できるデータを持っているとする。aの値がわからないので、それを推定したいとする。データは datay,dataxというペアデータとしてもっているとする。
    • 今、このdatay,dataxのデータセットについて、datay - (a ¥times datax^2)が、yの値の差であり、その自乗(datay - (a ¥times datax^2))^2をすべてのデータペアについて足し合わせた値がもっとも小さくなるような a の値を求めたいものとする
    • 最小にしたい値を計算する関数 f を次のように書く

f<-function(x,y,a) {
sum((y-a*x^2)^2)
}

    • aの推定にあたっては、この式と、推定を開始するためのaの初期値と、データセットをnlm関数に次のように渡す

nlm(f,3,x=datax,y=datay)

    • 関数fが受け取るべき変数x,y,aのうち、xとyとは、nlm関数に与えており、値の与えていないaについて、nlm関数の第2引数3が初期値として与えられたものと判断されて処理される
    • では、datax,datayにy=5 ¥times x^2なるデータを仮データとして作成して、上記を実行する

datax<-c(1:100)
datay<-3*datax^2
f<-function(x,y,a){
sum((y-a*x^2)^2)
}
nlm(f,3,x=datax,y=datay)

    • 結果の表示は以下の通り(aの推定値として3が得られている $estimate [1] 3)

$minimum
[1] 0

$estimate
[1] 3

$gradient
[1] 0

$code
[1] 1

$iterations
[1] 1


    • 計算の仮定を詳しく表示させるには、print.levelオプションで指定する。aの初期値に100を与えて、print.levelに2を指定した場合は次のようになる

> nlm(f,100,x=datax,y=datay,print.level=2)
iteration = 0
Step:
[1] 0
Parameter:
[1] 100
Function Value
[1] 1.929159e+13
Gradient:
[1] 397764871055

iteration = 1
Step:
[1] -100
Parameter:
[1] -2.376571e-14
Function Value
[1] 1.8453e+10
Gradient:
[1] -12301997932

iteration = 2
Step:
[1] 2.999998
Parameter:
[1] 2.999998
Function Value
[1] 0.00807352
Gradient:
[1] -1986.180

iteration = 3
Parameter:
[1] 2.999999
Function Value
[1] 0.004613245
Gradient:
[1] 1.005457e-05

連続した繰り返しが許容範囲内です
現在の繰り返しでおそらく解が得られたでしょう

$minimum
[1] 0.004613245

$estimate
[1] 2.999999

$gradient
[1] 1.005457e-05

$code
[1] 2

$iterations
[1] 3