私的関数 メモ代わり

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)
}