次数が未知の最小二乗法について

質問

最小二乗法により実験データを解析したいと考えています。
n回の測定でyとxを求め、
y＝ax^2+bx+c
という二次関数でフィッティングする場合、未知パラメータabcを擬似逆行列で求めるというところまでは分かったのですが、
y=ax^d+bx+c
というように、未知パラメータがabcに加えて、次数のdも未知の場合にはどのように解いたら良いのでしょうか。
どなたか分かる方、教えてください。

回答No.2

masudaya · ベストアンサー

Excelを用いるのであれば,a,b,c,dをパラメータとして
ソルバーを用いるのが良いのではないでしょうか.
まずは,ツール＝＞アドインでソルバーが入っているか確認して
入っていなければ,ソルバー アドインのところをチェック下さい.
シート上にパラメータa,b,c,d,をきまったセルに入れます.
xの値とyの値とax^d+bx+cの値を横に並べてy-(ax^d+bx+c)も並べます.
取得したx,yの分,同じ用に縦に並べます.
y-ax^d+bx+cのセルをsumsq()で平方和します.

この状態でツール=>ソルバーを起動させると,
目的セルは平方和を入れたセルで,目標値は最小値とします.
変化させるセルはパラメータa,b,c,dになります.
(各パラメータについて正のみなどの制約条件も入れられます.)
うまくいかない場合は,オプション設定で制限時間や反復回数を変化させると良いと思います.

補足

回答有難う御座います。
早速ソルバーで計算させてみました。(詳しい説明ありがとう御座いました。)
残差平方和がより小さくなるように、オプションで
制限時間や探索法など色々変えて見たのですが、
測定値となかなか一致しません。(全体的に右に少しシフトした感じに
なっています。)
現在、残差平方和が3.47719E-05ほどなのですが、
通常はどのくらいの残差平方和を一致の目安としているのでしょうか。
度々申し訳御座いませんが、回答のほど宜しくお願い致します。

回答No.4

masudaya

もともとの値がいくらか分かりませんが,
1程度の数値だとすると,5桁落ちはたいしたものと思います.
最小二乗法などの回帰曲線は近似曲線なので,測定点を通ることは
ほとんどありえません.

全体的に右に少しシフトした感じ=>これはx軸上に平行移動
したようなイメージなのでしょうか.
そうであれば,y=ax^d+bx+cを
y=a(x-e)^d+b(x-e)+c=a(x-e)^d+bx+c'
と平行移動させれば良いかと思います.
パラメータを1個増やすことになりますが...
ちなみに,最小二乗法ではパラメータ数は
データ数の半分程度が限界と何かで読んだことがあります.

お礼

回答有難う御座います。

>最小二乗法などの回帰曲線は近似曲線なので,測定点を通ることは
ほとんどありえません.

そうだったんですね。勉強不足でした。有難う御座います。

>全体的に右に少しシフトした感じ=>これはx軸上に平行移動
したようなイメージなのでしょうか.

そうです。傾きも少し大きくなっていますが・・。

>そうであれば,y=ax^d+bx+cを
y=a(x-e)^d+b(x-e)+c=a(x-e)^d+bx+c'
と平行移動させれば良いかと思います.
パラメータを1個増やすことになりますが...

有難う御座います。そうやればいいんですね。
しかし、そうすると求めたい式の形が変わってしまう
気がするのですが・・・。

>ちなみに,最小二乗法ではパラメータ数は
データ数の半分程度が限界と何かで読んだことがあります.

それは初耳でした。大変勉強になりました。
フィッティングするように色々と試して見たいと思います。
有難う御座いました。

回答No.3

TACOSS

最も大きな次数の項の寄与が大きい(i.e. y≒ax^d)のなら、Logを取って
　　log y ＝ d・log x + a
として最小二乗フィットさせれば、最大次数ｄを推定できると思います。
普通、実験結果からザっと次数を決定するときは、結果をLog-Logプロットして
その勾配を見る、ってしませんか？
最大次数が推定できたら、後は次数に応じて普通に最小二乗法を適用すれば
良いのでは。

補足

回答有難う御座います。
教えて頂いた通り
log y ＝ d・log x + a
として最小二乗フィットさせ、
dを求めましたところ、d=0.00679となったので、この次数を用いて
y = ax^d + bx + c
に最小二乗フィットさせてみました。
その結果、測定値と大きくずれてしまったので、探索方法などを変えて
フィットさせて行くうちに、dはd=2に近づいて行きました。
また、未知パラメータの初期値は測定値のグラフにおいて、excelの
近似曲線の追加で出た値を入力しています。
間違えている点などありましたら、ご指摘頂けると有難いです。

回答No.1

usokoku

dの値で大きくグラフが変化するでしょう。このような場合には、行列ではちょっと無理です。
ただし、
Y=aiX^i + a(i-1)X^(i-1) +... + a(1)x + a(0)
の形に直せれば、可能です。

なおせないとして
作図法で、エイヤーと、a,b,c,d を決定してください。
a,b,c,d のいずれかひとつを少し(普通は1/100くらい)変化させて、計算してみて、残差が減ったらば、変化後の値が正しい
という作業を飽きれるまで(山さが十分少なくなったと判断できるまで)つづけてくだ最。
今のパソコンですと早いですから、2－10秒でできるでしょう。
8086無wait , 無コプロで10分ぐらいの計算ですから。

お礼

早速の回答、有難うございます。
excelのソルバーで図を描かせて、そのような作業をしてみたいと思います。
有難う御座いました。