2006-02-18

汎用最適化関数 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関数の第２引数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に２を指定した場合は次のようになる



> 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

2006-02-18

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

2006-02-18

私的関数　メモ代わり

R 最適化

nls関数を用いて、テキストファイルからデータを読み取り、非線形最適化をし、その推定パラメタをファイル出力する関数(ごくごく私的な用途、メモ代わり)



est<-function(file){

data<-read.table(file,TRUE,"\t")

datax<-data$x

datay<-data$y
datax<-datax/100000

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

summary(result)

pars<-result$m$getPars()

pars

yset<-(1-result$m$getPars()["c"])*(1-result$m$getPars()["m"])+result$m$getPars()["c"]
pars["ysep"]<-yset

pars

outfile <-paste(file,"out.txt")
write(pars,outfile)

return(result)
}

さらに、返った結果を用いて、プロットする。重ねプロットにpar(new=T)を使うことなどの備忘録に。plotの引数cesは点のサイズ、xlim,ylimを共通にすることで、描図の範囲をそろえる。



plotest<-function(result){

estx<-c(-1000:2500)/1000

esty<-(1-result$m$getPars()["c"])*(1-result$m$getPars()["m"])*exp(-result$m$getPars()["k"]*estx)+result$m$getPars()["c"]

plot(estx,esty,cex=.05,xlim=c(-0.2,2.5),ylim=c(-0.1,1.3))

par(new=T)

x1<-c(rep(1,length(estx)))

plot(estx,x1,cex=.01,xlim=c(-0.2,2.5),ylim=c(-0.1,1.3))

par(new=T)

y1<-c(rep(0,length(esty)))

plot(y1,esty,cex=.01,xlim=c(-0.2,2.5),ylim=c(-0.1,1.3))

}

残差もプロットさせてみる



estLDwithPlotRes<-function(file){
plotestLD3<-function(result){

estx<-c(-1000:2500)/1000

esty<-(1-result$m$getPars()["c"])*(1-result$m$getPars()["m"])*exp(-result$m$getPars()["k"]*estx)+result$m$getPars()["c"]

plot(estx,esty,cex=.05,xlim=c(-0.2,2.5),ylim=c(-0.1,1.3))

par(new=T)

x1<-c(rep(1,length(estx)))

plot(estx,x1,cex=.01,xlim=c(-0.2,2.5),ylim=c(-0.1,1.3))

par(new=T)

y1<-c(rep(0,length(esty)))

plot(y1,esty,cex=.01,xlim=c(-0.2,2.5),ylim=c(-0.1,1.3))

}
plotestRes<-function(result,datax,datay){

estx<-datax

esty<-datay-((1-result$m$getPars()["c"])*(1-result$m$getPars()["m"])*exp(-result$m$getPars()["k"]*estx)+result$m$getPars()["c"])


plot(estx,esty,cex=.05,xlim=c(-0.2,2.5),ylim=c(-0.1,1.3))

par(new=T)

x1<-c(rep(0,length(estx)))

plot(estx,x1,cex=.01,xlim=c(-0.2,2.5),ylim=c(-0.1,1.3))

par(new=T)

y1<-c(rep(0,length(esty)))

plot(y1,esty,cex=.01,xlim=c(-0.2,2.5),ylim=c(-0.1,1.3))

}
data<-read.table(file,TRUE,"\t")

datax<-data$x

datay<-data$y
datax<-datax/100000

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

summary(result)

pars<-result$m$getPars()

pars

yset<-(1-result$m$getPars()["c"])*(1-result$m$getPars()["m"])+result$m$getPars()["c"]
pars["ysep"]<-yset

pars

outfile <-paste(file,"out.txt")
write(pars,outfile)
plot(datax,datay,cex=.01,xlim=c(-0.2,2.5),ylim=c(-0.1,1.3),col="red")

par(new=T)

plotestLD3(result)

par(new=T)

plotestRes(result,datax,datay)
return(result)
}