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

No.21004　nlsで推定できるモデル式について　　【晴】　2014/04/15(Tue) 19:59

交通量推計業務に携わる者です。
初歩的な質問なのでひょっとしたらどこかで回答がなされているかもしれませんが，ネット上をいろいろ検索しても，結局わからなかったのでお願いいたします。

質問は，「非線形最小二乗法によるパラメータ推定」のための関数nlsは，どのようなモデル式に対してもパラメータ推定が可能なのか？ということです。

私が推定しようとしているのは，交通工学では代表的なモデル式です。
(ロジットモデルのパラメータを距離Disに依存させた形のモデル式になっています。)

nls(P~1/(exp(-a*(Dis^b)*DT+c*log(Dis)+d)+1),Data,start=list(a=2.2,b=-1.0,c=0.5,d=0.5))

P(転換率(高速道路利用率)),Dis(距離),DT(高速ルートと一般道ルートの一般化時間の差)一般化時間は，所要時間と時間換算された通行料金の和です。

Rが利用できるかどうか知りたくて，当てはまりが良いはずのサンプルを作成して推定を試みました。
サンプルは，a=2.2,b=-0.964,c=0.442,d=0.552を上のモデル式に代入して得ました。このパラメータは先行研究で求められていた値です。

業務では，実績データを用いて新たにパラメータ推定することになります。

さて，実際に上のサンプルでパラメータ推定を試みたのですが，動きがありません。
はじめ，a=1.0,b=1.0,c=1.0,d=1.0を初期値としたときには，「勾配が特異です」とのエラーメッセージが出ました。
そこで，代入したパラメータに近い値を初期値として設定してみたのです。
すると，エラーメッセージは出なくなりましたが，反応がないのです。
もし，推定できない形のモデル式であれば，何らかのエラーメッセージが表示されるのでしょうか？

Rを始めたばかりでよくわからない状態で，計算が終わっているかどうかすら判断できないくらいのレベルで困惑しております。何卒よろしくお願い申し上げます。

参考までに，作成したプログラムをお示しします。

### Multinominal Logit estimation program (parameters depending on distance)

Data<-read.csv("D:/R/tenkanrei3.csv",header=TRUE)

nls(P~1/(exp(-a*(Dis^b)*DT+c*log(Dis)+d)+1),Data,start=list(a=2.2,b=-1.0,c=0.5,d=0.5))

result<-nls(P~1/(exp(-a*(Dis^b)*DT+c*log(Dis)+d)+1),Data,start=list(a=2.2,b=-1.0,c=0.5,d=0.5))

summary(result)

もしRでは推定できないのであれば，fortranプログラムでやってみるつもりです。

No.21005　Re: nlsで推定できるモデル式について　　【青木繁伸】　2014/04/15(Tue) 20:49

プログラムを見せられただけでは，データが提示されないと，追試できません。
一般的に言えば，データに欠陥がない限り，何でもあてはめできるでしょう。
ただ，アルゴリズムにより初期値に敏感なものもあるので，ある程度の試行錯誤，よりよい初期値の推定値を用いること，初期値に敏感でないアルゴリズムによるなど，工夫が必要なこともあるでしょう。

No.21006　Re: nlsで推定できるモデル式について　　【晴】　2014/04/16(Wed) 10:55

ありがとうございます。用いたデータは以下の通りです。

P Dis DT
0.407 10 5
0.694 10 10
　　中略
0.1 100 15
0.112 100 20
0.126 100 25

説明変数の値は，10<=Dis<=100,5<=DT<=25の範囲で5刻みであり，Pは，モデル式にa=2.2,b=-0.964,c=0.442,d=0.552を代入しただけの値(小数第4位で四捨五入)です。
よろしくお願いいたします。

No.21007　Re: nlsで推定できるモデル式について　　【青木繁伸】　2014/04/16(Wed) 11:18

21006 の P は「モデル式にa=2.2,b=-0.964,c=0.442,d=0.552を代入しただけの値(小数第4位で四捨五入)」なんでしょ？

summary によって，a, b, c, d にきわめて近い推定値が表示されているでしょう（Estimate の列）。

> summary(result)

Formula: P ~ 1/(exp(-a * (Dis^b) * DT + c * log(Dis) + d) + 1)

Parameters:
Estimate Std. Error t value Pr(>|t|)
a 2.1975635 0.0031932 688.2 <2e-16
b -0.9637574 0.0004493 -2144.8 <2e-16
c 0.4418875 0.0005869 753.0 <2e-16
d 0.5523523 0.0022497 245.5 <2e-16

Residual standard error: 0.0002683 on 91 degrees of freedom

Number of iterations to convergence: 4
Achieved convergence tolerance: 6.474e-07

plot(result$m$fitted(), Data$P) を表示すると，ほぼ一直線状に点が表示されるでしょう。
もし，「小数第4位で四捨五入」なんてけちなことをやらずにできうる限り正確な P を使えば，a, b, c, d は，P を作るときに使った値になるでしょう。もっとも，そのようなことはしてはいけないと，オンラインヘルプに書いてあります。そんなことをする意味がないから。
Warning: Do not use nls on artificial "zero-residual" data.

あてはめはできると思いますから，「うまくいくかどうか試す」なんてことしなくて，実際の P とDis, DT を使って nls で係数の推定をすればよいでしょう。それでうまくいかないときは，その理由，対応策を考えればよいだけの話です。

No.21008　Re: nlsで推定できるモデル式について　　【晴】　2014/04/16(Wed) 12:53

ご迷惑をおかけしましたが，助かりました。

Rの操作の基本を知らなかったために，上のプログラムをsourceで読み込めば，summaryの最後まで自動的に流れると思っていました。

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