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

No.02189　Rについて　　【だー】　2007/01/09(Tue) 17:13

はじめまして．
非線形回帰式のパラメータ推定を行う必要性に迫られ，
本HP(http://aoki2.si.gunma-u.ac.jp/Hanasi/StatTalk/nls.html)を参考に
初めてRを使ってのパラメータ推定を試みております．

Rで非線形回帰式のパラメータ推定を行う場合，
nls という関数を利用するということはわかったのですが，
Y=a*X^b+cといった”べき乗”の推定を行おうとすると，
「パラメータの初期値で勾配行列が特異です」というエラーが発生してしまいます．
いろいろと初期値を変化させてみても，推定できません．
”べき乗”部分が含まれていなければ問題なく推定できるのですが．．．

Rでは”べき乗”の推定はできないでしょうか．．．

※HP上の初期値の選び方を参考にさせていただきましたが，やはり上手くいきませんでした．
http://aoki2.si.gunma-u.ac.jp/Hanasi/StatTalk/optim.html

基本的なことだと思いますがRの初心者なため，よろしくお願いします．

No.02190　Re: Rについて　　【波音】　2007/01/09(Tue) 17:17

質問へ対する直接の答えではないのですが,私も以前に同じような質問をして教えていただいたことがあるので,もしかしたら参考になるかもしれません。

1787から始まるスレッドです。

No.02193　Re: Rについて　　【青木繁伸】　2007/01/09(Tue) 18:06

> Rでは"べき乗"の推定はできないでしょうか

そんなこと，あるわけがありません。

[パラメータの初期値で勾配行列が特異]というのは，要するに初期値が悪いと言うことです。

うまく初期値を捜すかというのは，いくつか方法があります。

1. エクセルで，パラメータを入力するとそのときの曲線をデータ点と共に描くというようなワークシートを作って，どのパラメータをどのように動かしたらどのような曲線になるかを把握しながらやるとよい。でたらめに捜すよりは効率的。

2. シンプレックス法を使う

3. 変数変換をして重回帰分析として解いて，得られるパラメータを逆変換する

4. 根気よく捜す(^_^)

No.02200　Re: Rについて　　【だー】　2007/01/09(Tue) 21:01

エクセルなどを組合せて根気よく探してみます．
ご回答ありがとうございました．

No.02202　Re: Rについて　　【だー】　2007/01/09(Tue) 21:29

度々申し訳ありません．
エクセルなどにより初期値をいろいろと探しているのですが，上手くいきません．

そこで，Y=1*X＾3　という関係のデータを作成し，
Y=(a*X＾b)+c　a=1，b=3，c=0　という初期値で，パラメータ推定を行ってみました．
しかし，やはり[勾配が特異です]というエラーが発生しました．

私のやり方に何か間違いがあるのでしょうか．．．
申し訳ありませんが，よろしくお願いします．

ちなみに，作成したコマンドは次のとおりです．
> X<-c(1,  8,  27,  64,  125,  216,  343, 512, 729, 1000)
> Y<-c(1,  2,  3,  4,  5,  6,  7, 8, 9, 10)
> dat<-data.frame(X,Y)
> nls(Y~(a*X^b)+c,dat,start=list(a=1,b=3,c=0))
以下にエラーnls(Y ~ (a * X^b) + c, dat, start = list(a = 1, b = 3, c = 0)) : 
        勾配が特異です

No.02203　Re: Rについて　　【青木繁伸】　2007/01/09(Tue) 22:13

どう見ても，X と Y が逆ではないかと思いますが?
また，このデータはテストデータで，Y=X^3+0 ですよね。というか，そのように作ったんですよね。。。って，start 見ると，あきらかでですが。
ついでに言えば，nls のオンラインヘルプを見ると，誤差のないテストデータで試さないようにと書いてあったような気がします

Warning

Do not use nls on artificial "zero-residual" data.

ということで，以下のようにすれば体験できると思いますよ。
Y<-c(1, 9, 28, 65, 125, 218, 364, 458, 800, 989)
X<-c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10)
dat<-data.frame(X,Y)
nls(Y~(a*X^b)+c,dat,start=list(a=1,b=3,c=0))
実行例は以下の通り
> Y<-c(1, 9, 28, 65, 125, 218, 364, 458, 800, 989)
> X<-c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10)
> dat<-data.frame(X,Y)
> nls(Y~(a*X^b)+c,dat,start=list(a=1,b=3,c=0))
Nonlinear regression model
  model:  Y ~ (a * X^b) + c 
   data:  dat 
      a       b       c 
 1.0225  2.9948 -0.2841 
 residual sum-of-squares: 8306

Number of iterations to convergence: 2 
Achieved convergence tolerance: 2.511e-06 

No.02206　Re: Rについて　　【だー】　2007/01/10(Wed) 09:44

青木先生
すいません．
XとYが逆でした．
お恥ずかしいかぎりです．．．