「統計学関連なんでもあり」の過去ログ---027

★ Rでの非線形最小二乗法について ★

2234.　Rでの非線形最小二乗法について　furuki　2004/02/05 (木) 14:47
└2236.　Re: Rでの非線形最小二乗法について　青木繁伸　2004/02/05 (木) 16:11
　└2237.　Re^2: Rでの非線形最小二乗法について　青木繁伸　2004/02/05 (木) 16:49

2234.　Rでの非線形最小二乗法について　furuki　　2004/02/05 (木) 14:47

いつもお世話になってます。HPでの例を参考に
Rで非線形最小二乗法を行っているのですが，
エラーメッセージが例と異なっていて困っています。
どこがいけないのか教えていただけないでしょうか？
数値以外は例をコピーアンドペーストしたのですが…

x y 1 35 0 2 37 1 3 39 0 4 41 0 5 43 0 6 45 1 7 47 1 8 49 2 9 51 5 10 53 7 11 55 3 12 57 3 13 59 12 14 61 8 15 63 14 16 65 13 17 67 10 18 69 18 19 71 7 20 73 9 21 75 7 22 77 1 23 79 3 24 81 2 25 83 0 26 85 0 27 87 1 > nls(y ~ a/(1+b*exp(-c*x))+d, dat, start=list(a=1, b=1, c=1, d=1)) Error in nlsModel(formula, mf, start) : singular gradient matrix at initial parameter estimates
また，初期値は勘で入力していくしか方法がないのでしょうか？
よい方法があるのでしたらこれも教えてください。

　　　　　[このページのトップへ]

2236.　Re: Rでの非線形最小二乗法について　青木繁伸　　2004/02/05 (木) 16:11

> エラーメッセージが例と異なっていて困っています。
> どこがいけないのか教えていただけないでしょうか？

エラーメッセージの意味は，「初期値が悪い」ということです。

> 数値以外は例をコピーアンドペーストしたのですが…

全く同じような例なら，うまくいくこともあるでしょうが。
うまくいったとしても，偶然でしょう。

> また，初期値は勘で入力していくしか方法がないのでしょうか？
> よい方法があるのでしたらこれも教えてください。

うまい方法ではなくて，当たりをつける方法はあります。

まず，データのグラフを描くべきです。Excel などを使うのがいいでしょう。
あなたのデータの場合，これはあなたが当てはめようとしているデータでしょうか？
予測式が，y ~ a/(1+b*exp(-c*x))+d ということは，y は d ～ d+a の範囲にある。いわゆるs字状曲線です。x, y のデータとはずいぶん違うでしょう？

先のデータの y は，累積和を作らないといけないことになります，y は 0,1,1,1,1,2,3,5,10,17,...,122,125,127,127,127,128 となるんでしょうね。
そして，初期値のうち簡単にわかるものは，d ≒ 0, a ≒ 130。b,c を適当に決めて（変えて）さっき描いた x, y のグラフの中に予測値の曲線を加えてみましょう。b, c を式の中に直接書かずに（a, d もですが），セルに書き込んでおいてそれを用いて計算するようにします。
b, c のセルの数値を変えると予測値の曲線と y の値の一致の具合がよくわかるでしょう。

ここまでくればしめたもの，Excel のソルバーを使ってもよいし，R の nls でもよいでしょう。
ところが，R の nls はなかなか気むずかしいようで，Excel のソルバーだと多少のいい加減な初期値でも収束するようではあります。
結局のところ，試行錯誤の末
> nls(y ~ a/(1+b*exp(-c*x))+d, dat, start=list(a=130, b=360000, c=0.2, d=-0.25)) Nonlinear regression model model: y ~ a/(1 + b * exp(-c * x)) + d data: dat a b c d 1.301335e+02 5.438734e+05 2.057295e-01 4.545957e-01 residual sum-of-squares: 85.2488
となりました。（d はモデルに入れなくてもよいようですね）

　　　　　[このページのトップへ]

2237.　Re^2: Rでの非線形最小二乗法について　青木繁伸　　2004/02/05 (木) 16:49

私も R はまだまだ勉強中ですので，なかなかですが。

関数の最大化・最小化を行う optim というのがあり，初期値を探すのはこれを使うといいみたいですね。

先ほどのデータで，
> x [1] 35 37 39 41 43 45 47 49 51 53 55 57 59 61 63 65 67 69 71 73 75 77 79 81 83 [26] 85 87 > y [1] 0 1 1 1 1 2 3 5 10 17 20 23 35 43 57 70 80 98 105 [20] 114 121 122 125 127 127 127 128
が定義されていて，

最小化すべき関数が
> resid function(par) { yhat <- par[1]/(1+par[2]*exp(-par[3]*x))+par[4] sum((y-yhat)^2) }
つまり，予測誤差の二乗和を最小にする。引数は，モデル式に含まれるパラメータベクトル。x, y は大域変数として利用。

以上の用意ができたら，

> optim(c(130, 30000, 0.2, 0), resid)
とする。第1引数はパラメータの初期値ベクトル，第2引数は最小化する関数名。

うまく計算が終われば，以下のように結果が示される。

$par [1] 136.3372705 33140.5711068 0.1639743 -5.1328886 $value [1] 391.2728 $counts function gradient 373 NA $convergence [1] 0 $message NULL
これを初期値として nls を行えばいいのでしょう。

　　　　　[このページのトップへ]