Ｒを使った非線形最小二乗について

質問

はじめまして。
Ｒを使った非線形最小二乗のフィッティングについて、ご教授お願いします。
（カテゴリーが違うようでしたらご指摘お願いします）

二つのデータの組Ｙ＝（y1,y2,y3,・・・）、Ｚ＝(z1,z2,z3,・・・)について、これらをf(x) = α / { 1 + β exp(-γx) }の関数でフィッティングしたいと考えています。
ここで、f(x)の未知パラメータα、βはＹとＺで同じ値を持ち、γはＹとＺに固有のパラメータとして、フィッティングしようとしています。

以上のことをＲのnls関数を用いて行いたいが、どのようにすればよいか。

以下の流れでフィッティングしてみましたが、上手くいきませんでした。
(1)　Ｙをf(x)でフィッティングして、α、β、γを求める。
(2)　(1)のα、βを用いてＺをf(x)でフィッティングする。

そもそも、αとβを共通の値と考えること自体が間違いというご指摘があるかもしれませんが、
あくまでもαとβが共通の値と考えたとき、ＹとＺをf(x)で最も上手くフィッティング（Ｙの残差とＺの残差の和が最小になる）したいと考えています。

よろしくお願いいたします。

回答No.2

rabbit_cat

＞f(x)の未知パラメータα、βはＹとＺで同じ値を持ち、γはＹとＺに固有のパラメータとして、フィッティングしようとしています。

nls関数では、このようなことはできません。
＞(1)　Ｙをf(x)でフィッティングして、α、β、γを求める。
＞(2)　(1)のα、βを用いてＺをf(x)でフィッティングする。
お気づきだと思いますが、このやり方では、パラメータα、βの推定にYのデータしか使っていないわけで、上記のフィッテイングになってないです。
（本来は、パラメータα、βの推定にはYとZの両方のデータを使わないとおかしい）

一応、真面目に式を立てて考えれば、このような条件での、最小二乗回帰の式を導出することは可能です。（この条件の回帰ができるように、nls関数相当の関数を自前で実装することになる）

そんなことやってられん、ということであれば、階層ベイズモデルだと思って、MCMC（マルコフ連鎖モンテカルロ）してしまうのが一番簡単です。
MCMCであれば、質問の例のように観測データの一部だけが依存するパラメータがある場合とか、パラメータ間に複雑な依存関係があろうが、あるいは、欠損データが大量にあろうが、何も考えずにそのまま推定できます。
（人間が何も考えなくてよくなる代わりに、計算時間は大幅に増えてしまいますが。。）

RでMCMCというと、「WinBUGS」とか、最近だと「Stan」とかが使われているんでしょうか。
これらのキーワードで調べてみるとよいのでは。

お礼

返信遅くなりました。申し訳ありません。

やりたいことを解決する、単純な方法が無いことは理解できました。
回帰分析などやったことがないズブの素人なので、回答いただいた情報を読んでもまだピンと来ませんが、No.1様の回答と併せて勉強してみて、一番良い方法を考えてみたいと思います。

ご回答ありがとうございました。

回答No.1

kamiyasiro

質問の趣旨から外れ、回答になっていませんが、ご参考まで。

デフォルトの非線形回帰（nls）は、よくエラーを起こすので、
ライブラリーminpackのnls.lmを使って、レーベンバーグ・マルカート法で
非線形回帰することをお勧めします。

それと、最小２乗にこだわりますか。
曲線が立っているところは、ペナルティを緩めないとフィッティングしませんよ。

下は、minpackを使って、変動係数最小化で非線形回帰するスクリプトです。
データは1列目が横軸、2列目は縦軸として下さい。

# 代表的な減衰関数への近似を行う。
# power(冪乗) : 両対数グラフで直線 a*(1-x)^b
# exponetial(指数) :y片対数グラフで直線 a*exp(-b*x)
# logarithmic(対数) :x片対数グラフで直線 a*log((x+b)/(xmax+b))
# logistic(ロジスティック):1/(1+exp(a+bx))
#
# これらの関数はパラメータbに対して非線形であるため、
# levenberg-marquardt法により、データからパラメータを最小自乗推定する。
#
par(mfrow=c(2,2))
require(graphics)
library(minpack.lm)
#
# データの読み込み
#
data <- read.csv("******.csv")
x <- data[,1]
parax <- min(data[,1])
y <- data[,2]
paray <- mean(data[,2])
#
#
para <- NULL
para.c <- NULL
for (i in 1:4){
#
# モデル関数の定義
#
if(i==1){func <- function(para,xx){para$d+para$a*(xx-para$c)^(-1*para$b)}}
if(i==2){func <- function(para,xx){para$d+para$a*exp(-1*para$b*(xx-para$c))}}
if(i==3){func <- function(para,xx){para$d+para$a*log((xx+para$b)/(para$c+para$b))}}
if(i==4){func <- function(para,xx){para$d/(1+exp(para$a+para$b*(xx-para$c)))}}
#
# 残差関数を定義(変動係数の最小化)
#
res <- function(para,observed,xx){(observed-func(para,xx))/func(para,xx)}
#
# パラメータ初期値を設定
#
if(i==1){para.s <- list(a=1,b=1,c=0, d=paray)}
if(i==2){para.s <- list(a=1,b=1,c=parax,d=paray)}
if(i==3){para.s <- list(a=0,b=1,c=100, d=paray)}
if(i==4){para.s <- list(a=0,b=0,c=parax,d=paray)}
#
# レーベンバーグ・マーカート法で係数を求める
#
result <- nls.lm(par=para.s,fn=res,observed=y,xx=x,control=nls.lm.control(nprint=1,maxiter=200))
result
#
# 求められた係数で近似線を引く
#
para$a <- result$par$a
para$b <- result$par$b
para$c <- result$par$c
para$d <- result$par$d
para.c <- cbind(para.c,para)
m.corr <- round(cor(func(para,x),y),digit=4)
t <- seq(min(x),max(x),by=1)
y.pred <- func(para,t)
plot(data,pch=20)
lines(t,y.pred,col="black")
mtext(m.corr,side=3,line=0,at=NA)
}
#

お礼

返信が遅くなりました。申し訳ありません。

頂いたスクリプトを参考に、どのような方法でフィッティングするのが一番いいか、勉強して考えてみたいと思います。
ご回答ありがとうございました。