最適化関数 nls



  • Rの Htmlヘルプ"nls{stats}","nlsModel{stats}"も参考のこと
  • ごくごく簡単に
    • 今、y=a¥times e^{k¥times x}+bに近似できるデータを持っているとする。a,b,kの値がわからないので、それを推定したいとする。データは datay,dataxというペアデータとしてもっているとする。
    • このようなとき、近似式と、推定の初期値を次のように与える。trace=TRUEは実行時の標準出力のオプションである

result<-nls(datay ~ a * exp(-k*datax)+b,start = list(a=1,k=1,b=1),trace=TRUE)

    • 結果の表示には、次のようにする

summary(result)

    • また、推定変数を取り出すには、次のようにする

result$m$getPars()

    • 今、datax,datayにペアデータが入っているものとする(関数の解説文書にあるように、式の通りに算出したデータセットを用いて実行すると、推定論理の都合でうまく進まないので、用いるデータセットは、「理論値」からある程度のずれがあるようにする)
    • 今、2カラムのタブ区切りテキストファイル"hoge.txt"の第1カラムにxの値が、第2カラムにyの値が入っているものとして、そこからの処理を以下に書く。1行目はヘッダ行で、それぞれx,yと記されているものとする

> data<-read.table("hoge.txt",TRUE,"\t")
> datax<-data$x
> datay<-data$y
> result<-nls(data$y ~ a*exp(-k*datax) + b,start = list(a=1,b=1,k=1),trace=TRUE)
202.2528 : 1 1 1
3.777557 : 0.3029115 0.4300912 1.5836115
0.1444613 : 0.4318135 0.3402358 3.2362493
0.0291222 : 0.5121943 0.3164979 3.5311545
0.02900385 : 0.5158914 0.3169041 3.5554486
0.02900286 : 0.5161296 0.3169499 3.5598961
0.02900283 : 0.5161722 0.3169585 3.5607202
0.02900283 : 0.5161801 0.3169601 3.5608725
0.02900283 : 0.5161816 0.3169604 3.5609006
> summary(result)

Formula: data$y ~ a * exp(-k * datax) + b

Parameters:
Estimate Std. Error t value Pr(>|t|)
a 0.516182 0.005123 100.75 <2e-16 ***
b 0.316960 0.001155 274.43 <2e-16 ***
k 3.560901 0.059238 60.11 <2e-16 ***

    • -

Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.01213 on 197 degrees of freedom

Correlation of Parameter Estimates:
a b
b -0.0006447
k 0.5964732 0.537

> result$m$getPars()
a b k
0.5161816 0.3169604 3.5609006
> result$m%getPars()["a"]
a
0.5161816

実行条件nls.controlへの値の与え方


> result<-nls(y~a*exp(b*x)+c, start = list(a=0,b=0.5,c=0.5),trace=TRUE,control=nls.control(maxiter = 100,tol = 1e-05, minFactor = 1/2096))