はじめまして。
Rを使った非線形最小二乗のフィッティングについて、ご教授お願いします。
(カテゴリーが違うようでしたらご指摘お願いします)
二つのデータの組Y=(y1,y2,y3,・・・)、Z=(z1,z2,z3,・・・)について、これらをf(x) = α / { 1 + β exp(-γx) }の関数でフィッティングしたいと考えています。
ここで、f(x)の未知パラメータα、βはYとZで同じ値を持ち、γはYとZに固有のパラメータとして、フィッティングしようとしています。
以上のことをRのnls関数を用いて行いたいが、どのようにすればよいか。
以下の流れでフィッティングしてみましたが、上手くいきませんでした。
(1) Yをf(x)でフィッティングして、α、β、γを求める。
(2) (1)のα、βを用いてZをf(x)でフィッティングする。
そもそも、αとβを共通の値と考えること自体が間違いというご指摘があるかもしれませんが、
あくまでもαとβが共通の値と考えたとき、YとZをf(x)で最も上手くフィッティング(Yの残差とZの残差の和が最小になる)したいと考えています。
よろしくお願いいたします。
A 回答 (3件)
- 最新から表示
- 回答順に表示
No.3
- 回答日時:
少し考えたら、EM法みたいにして、α/βと、γを交互に推定すれば、nls関数でもいけそうな気がする。
http://ja.wikipedia.org/wiki/EMアルゴリズム
1. まず、Yをf(x)でフィッティングして、α、β、γを求める(αとβの初期値を決める)
2. α,βを固定して、YとZそれぞれにnls関数を使って、γYとγZを求める
3. γYとγZを固定して、YとZ両方を使って、αとβを求める(nls関数で)
で、2と3を、収束するまで繰り返す、とすれば良いでしょう。
No.2
- 回答日時:
>f(x)の未知パラメータα、βはYとZで同じ値を持ち、γはYとZに固有のパラメータとして、フィッティングしようとしています。
nls関数では、このようなことはできません。
>(1) Yをf(x)でフィッティングして、α、β、γを求める。
>(2) (1)のα、βを用いてZをf(x)でフィッティングする。
お気づきだと思いますが、このやり方では、パラメータα、βの推定にYのデータしか使っていないわけで、上記のフィッテイングになってないです。
(本来は、パラメータα、βの推定にはYとZの両方のデータを使わないとおかしい)
一応、真面目に式を立てて考えれば、このような条件での、最小二乗回帰の式を導出することは可能です。(この条件の回帰ができるように、nls関数相当の関数を自前で実装することになる)
そんなことやってられん、ということであれば、階層ベイズモデルだと思って、MCMC(マルコフ連鎖モンテカルロ)してしまうのが一番簡単です。
MCMCであれば、質問の例のように観測データの一部だけが依存するパラメータがある場合とか、パラメータ間に複雑な依存関係があろうが、あるいは、欠損データが大量にあろうが、何も考えずにそのまま推定できます。
(人間が何も考えなくてよくなる代わりに、計算時間は大幅に増えてしまいますが。。)
RでMCMCというと、「WinBUGS」とか、最近だと「Stan」とかが使われているんでしょうか。
これらのキーワードで調べてみるとよいのでは。
返信遅くなりました。申し訳ありません。
やりたいことを解決する、単純な方法が無いことは理解できました。
回帰分析などやったことがないズブの素人なので、回答いただいた情報を読んでもまだピンと来ませんが、No.1様の回答と併せて勉強してみて、一番良い方法を考えてみたいと思います。
ご回答ありがとうございました。
No.1
- 回答日時:
質問の趣旨から外れ、回答になっていませんが、ご参考まで。
デフォルトの非線形回帰(nls)は、よくエラーを起こすので、
ライブラリーminpackのnls.lmを使って、レーベンバーグ・マルカート法で
非線形回帰することをお勧めします。
それと、最小2乗にこだわりますか。
曲線が立っているところは、ペナルティを緩めないとフィッティングしませんよ。
下は、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)
}
#
返信が遅くなりました。申し訳ありません。
頂いたスクリプトを参考に、どのような方法でフィッティングするのが一番いいか、勉強して考えてみたいと思います。
ご回答ありがとうございました。
お探しのQ&Aが見つからない時は、教えて!gooで質問しましょう!
似たような質問が見つかりました
- 数学 1変数関数に陰関数ってあるんですか? 1変数関数は f(x)=xの式 f(x)はxの値で決まるもの( 4 2023/05/08 18:47
- 数学 2013 慶応(らしいです) 1 2022/06/14 21:15
- 数学 ラグランジュの未定乗数法を用いる問題 3 2023/05/15 14:48
- 数学 【高1 数学Ⅰ 二次関数】 二次関数 f(x)=x^2-4ax+8a がある。ただし、aは正の定数と 3 2022/07/23 15:46
- 統計学 統計学の問題です よろしくお願いします 回帰直線 次のデータから集計表を作成し,以下の問いに答えよ。 2 2023/01/31 23:36
- 統計学 統計学の問題です よろしくお願いします 回帰直線 次のデータから集計表を作成し,以下の問いに答えよ。 1 2023/01/31 18:55
- 高校 合成関数の定義域につきまして 1 2022/05/18 17:26
- 物理学 物理の惑星の問題 2 2023/03/21 18:51
- 数学 写真の数学の質問です。 「最小値がf(x)=x^2++x+a>g(x)=x^2+x+2aになる理由」 4 2023/01/04 11:47
- 数学 条件付き極値問題といわれる問題です。ラグランジュの乗数法 について、質問したいことがあります。 条件 3 2023/05/15 21:38
おすすめ情報
デイリーランキングこのカテゴリの人気デイリーQ&Aランキング
-
決定係数がマイナスになる例っ...
-
切片あり回帰と切片なし回帰
-
Yハットの出し方やミュートと...
-
修正済み決定係数(R2乗)がマ...
-
線形相関係数
-
ある1点で傾きが急激に変化する...
-
回帰式と近似式について
-
Rを使った非線形最小二乗について
-
統計学について質問です!
-
残差について
-
EXCELで両対数を取った重回帰分...
-
回帰関係の有意性と回帰係数の...
-
非線形回帰分析の定義が分かり...
-
重回帰分析で偏回帰係数を全て...
-
相関分析の相関係数と重回帰分...
-
二つのデータの波形が似てるか...
-
有意に大きい(小さい)とは?
-
帳票出力に関して
-
x^2+y^2-x-y=0 の実すうかいを...
-
2つのDataTableをJoin
マンスリーランキングこのカテゴリの人気マンスリーQ&Aランキング
-
決定係数がマイナスになる例っ...
-
切片あり回帰と切片なし回帰
-
回帰式と近似式について
-
ある1点で傾きが急激に変化する...
-
修正済み決定係数(R2乗)がマ...
-
回帰水を売ってる会社大丈夫か
-
Excel分析ツールでのポアソン回...
-
原点強制通過させたときの相関係数
-
Yハットの出し方やミュートと...
-
最小二乗法の傾きと切片について
-
統計用語?
-
残差について
-
重回帰分析で偏回帰係数を全て...
-
実験データの分析について
-
numbersで重回帰分析をしたい
-
ロジスティック回帰分析におけ...
-
エクセル 重回帰 グラフ
-
相関分析の相関係数と重回帰分...
-
重回帰分析の定数は、どっちの...
-
確率統計です。
おすすめ情報