一般化最小二乗法の定式化

質問

現在，一自由度系の振動問題において(1)式の微分方程式の定式化を行い最小二乗法を用いたM(質量),B(粘性),K(剛性)のパラメータ同定を行っています。
F=M(dx^2/dy^2)+B(dx/dy)+Kx・・・(1)
F:力(測定値)，x:変位(測定値)
dx/dy：速度（差分による計算値）
dx^2/dy^2：加速度（差分による計算値）
ここで，式誤差の影響を減らすために，一般化最小二乗法を用いたM,B,Kの算出を行おうと考えていますが，文献を調べても簡単に記載されているものが見つかりません。もしよろしければ，わかりやすいサイト，文献やC言語などのソースが手に入る場所などがあれば教えて頂けますでしょうか。
よろしくお願いします。

回答No.3

stomachman · ベストアンサー

実務でお困りなのだろうと思います。
　対象は減衰項を持つ振動子の系ですね。Fは強制的に与える外力でしょう。こんな間接的なやりかたで質量を精度よく測りたい。しかもばね定数や粘性も未知だという状況ですから、例えば、ナノサイズのカンチレバーにくっついたDNAの質量を測ろう、なんて話かも知れません。だとすると、測定値に大きな誤差が入るのはしょうがないことです。

●測定誤差が気になるような場合には、数値微分はひじょーにアブナイ。ってか大抵使い物になりません。(1)式右辺の３つの項の測定誤差が互いに強い相関を持つ訳ですから、理論上も無茶です。
　こういうときは積分形を考える方がずっと良い。積分によってノイズが平滑化されるからです。(1)式をyで積分(0～yまでの定積分）すれば、
∫Fdy = M∫x''dy+B∫x'dy+K∫xdy
= Mx'(y)+Bx(y)+K∫xdy-(Mx'(0)+Bx(0))
もう一回yで積分(0～yまでの定積分）すれば、
∫∫F(dy)^2 =M∫x'dy+B∫xdy+K∫∫x(dy)^2-(Mx'(0)+Bx(0))∫dy
=Mx(y)+B∫xdy+K∫∫x(dy)^2-(Mx'(0)+Bx(0))y-Mx(0)
ですから、
P=-(Mx'(0)+Bx(0))
Q=-Mx(0)
とおくと
∫∫F(dy)^2=Mx(y)+B∫xdy+K∫∫x(dy)^2+Py+Q
であり、線形モデルになっています。そして、測定した(y, F(y), x(y))から数値積分で∫∫F(dy)^2 , ∫xdy, ∫∫x(dy)^2を計算すると（これらの項に関する誤差がyについて独立でなくなってしまうけれども）、線形最小二乗法で近似解(M0, B0, K0, P0, Q0)が得られるでしょう。測定誤差がごく小さい場合、あるいは測定回数が多い場合には、これで十分です。

●測定回数が制限されていて測定誤差が大きいときには、もっとまじめに、測定誤差の独立性を壊さないように扱うべきです。そのためには微分方程式を陽に解いてモデルを作る必要があります。幸い、(1)式のような線形微分方程式は(たとえば演算子法やラプラス変換を使うと)すんなり解ける。話を簡単にするために、y<0のときはF(y)=x(y)=0である、とすると、(1)の解は
x(y) = m(M, B, K, y)*F(y)
という形に書けます。ここに、"*"は畳み込み積分(convolution)です。
　（ちなみに制御工学の用語で言うと、mは系の「インパルス応答（F(y)=δ(y)のときのxの振る舞い）」あるいは「伝達関数(modulation transfer function）」、そして、この系は「二次遅れ系」と呼ばれます。）

　さて、無視できない測定誤差がx(y)とF(y)の両方に含まれている場合には、話がかなり難しくなります。yの測定誤差まで考えるともっと大変なことに。

　なのでとりあえず、ここではyとF(y)の測定値は信用できるものとしましょう。すなわち、xの測定値をXと書くことにして、Xに入っている誤差εだけを考慮しようという訳です。
ε(y)=X(y)-x(y)
とすると、
X(y) = m(M, B, K, y)*F(y)+ε(y)
である。もしεは仰る所の「ホワイトノイズ」（つまり、ガウス分布に従う、測定ごとに独立な確率変数）だと考えてよいのなら、非線形最小二乗法で
E=Σ(|ε(y)|^2)
を最小化する(M, B, K）を探索すれば良い。（これをfittingと言います。また、探索の出発値として、上記の近似解(M0, B0, K0)が利用できます。）
　ここで忘れちゃいけないのは、fittingの結果得られたモデルと測定値とずれ、すなわち残差
ε(y)=X(y)-m(M, B, K, y)*F(y)
が実際に「ホワイトノイズ」らしい性質を持つ事を検証する、ということです。（単にプロットしてみるだけでも意味があります。）

　この検証をやってみると、例えば「|ε(y)|が(x(y)+C)に比例して激しく変動してた」なんてこともしばしばあります。つまり、
γ(y)=ε(y)/(x(y)+C)
とすると、γ(y)こそが「ホワイトノイズ」である。そういう時には、
X(y)+C = (m(M, B, K, y)*F(y)+C)(1+γ(y))
従って
ln(X(y)+C) = ln((m(M, B, K, y)*F(y))+C) + ln(1+γ(y))
ここで、|γ(y)|<<1であれば
ln(1+γ(y))≒γ(y)
なので、
ln(X(y)+C) = ln((m(M, B, K, y)*F(y))+C) + γ(y)
というモデルが出来ますから、これを使って非線形最小二乗法をやれば良い。
　その他いろんな場合があり得ますけど、誤差が「ホワイトノイズ」になるようにモデルを変形してやる訳です。で、これが仰る所の「一般化最小二乗法」「残差を白色化」ということでしょう。
　なお、最小二乗法の教科書としては、いつも「中川,小柳：最小二乗法による実験データ解析,東京大学出版会」をお薦めしています。

●ところでこの話は、単に最小二乗法の問題というよりも、機械工学（制御工学、振動工学）の問題と捉えるのが適切でしょう。その観点からは、モデルが実際の系をきちんと表しているかどうかについて、常に留意すべきです。たとえば、外力Fや変位xを測る手段自体が遅延などの応答特性を含んだ系である、ということがしばしばあるし、考慮していない要因に由来する系統誤差が入ることも多い。残差を調べてみて、モデルと測定結果が系統的にずれるようなら、モデルが実際の系に合っていない可能性があります。
　その場合、より正しいモデルを構築しなくちゃいけません。どう直すべきか心当たりがない時には、とりあえず系が線形であると仮定する。と、モデルは未知の線形伝達関数m(y)によって
x(y) = m(y)*F(y)
と表せます。で、既知の入力Fに対する出力xを測定して、系の伝達関数mをまるごと推定することを考える。古典制御理論の逆問題(inversion。この場合はdeconvolution)であり、イチオウ数値的に解くことができます。（これも線形最小二乗法の一種の応用です。しかし逆問題はノイズに極めて敏感なので、測定回数をうんと増やして情報量を増やす必要があります。）得られたmを物理的に解釈し直します。その結果、これまで無視していた影響要因が発見され、M, B, K, その他の要因に関するパラメータを含むモデルが構成できたとすると、xの測定値からこれらのパラメータを知る手段が得られるかも知れません。（得られないかも知れません。）
　系に非線形性がある場合にも、ある条件下（たとえば、xの振幅がある範囲内に入るという条件下）でなら線形で近似できることも多いので、条件を絞って解析する手があります。

● いずれにせよ、他人が作ったプログラムを拾って来てなんとかする、というような話ではありませんぜよ。

お礼

stomachmanさん
非常に丁寧に回答して頂きありがとうございました．知らなかったことばかりで私が不勉強であることを認識致しました．

これらを受け，まず提案頂いた積分形の解法について行いたいと考えています．

そこで基本的な質問になりますが，
∫Fdy = M∫x''dy+B∫x'dy+K∫xdy
0～yまでの定積分をすると教えいただいていますが，yを変数としたまま定積分を行う方法は一般的なのかについてです．

これまでの私の理解では定積分を行う際の積分範囲は0～Y(Yは定数)として考えるものだと考えていました．そうすると∫x''dy⇒x'(Y)，∫x'dy⇒x(Y)と∫xdyは定数になり，もう一度0～yまでの定積分する際には教え頂いた∫∫F(dy)^2 =Mx(y)+B∫xdy+K∫∫x(dy)^2-(Mx'(0)+Bx(0))y-Mx(0)
の式が算出できなくなります。

tを変数として微分方程式を定積分していく方法に関して，参考になる文献などありましたら教えて頂ければと思います．

次に，∫∫F(dy)^2=Mx(y)+B∫xdy+K∫∫x(dy)^2+Py+Q
において，測定値より∫∫F(dy)^2と∫∫x(dy)^2を計算する際はまず，∫Fdyと∫xdyつまり，0～y（時間）におけるF(力)とx（変位）との面積を定数として算出し，もう一度0～yの定積分を行う際は∫Fdyと∫xdyは定数となっているためにyを掛けるということになると考えます．
もし考え方に間違いがあるようであればご指摘頂ければうれしいです。

以上，よろしくお願い致します．

回答No.5

stomachman

ANo.4のコメントについてです。

　まず訂正。ANo.4の最後の式の左辺は
えふ(わい)=…
じゃなくて、正しくは
げふ(わい)=…
でげす。

> 「X(Y）」と 「えくす(わい）」の関係が独立ではなくなってしまう

　ご指摘の通りです。ANo.3で申し上げた通り、真面目にやるのなら、「測定誤差の独立性を壊さないように扱うべきです。そのためには微分方程式を陽に解いてモデルを作る必要があります。」でもANo.3前半にある手抜きのやり方をしても、「測定誤差がごく小さい場合、あるいは測定回数が多い場合には、これで十分」である。なぜなら、「積分によってノイズが平滑化されるから」です。

　つまり、この手抜きが通用するのは、＜測定データが、モデルと良く合っている場合＞に限られる。なので、fitting自体はANo.3の前半の手抜き法を使って行うにしても、微分方程式を陽に解いてモデルを作る、というところまではまじめにやって、fittingの結果得られたパラーメタをモデルに代入したものと、測定データとのずれ（残差）をプロットしてみるべきです。

お礼

ご回答ありがとうございました。
「微分方程式を陽に解いてモデルを作る」というところからはじめるが必要ということがわかりました。
再検討してみたいと思います。

回答No.4

stomachman

ANo.3のコメントについてです。

> yを変数としたまま定積分を行う方法は一般的なのか

　これはごく普通の定積分です。積分変数のyと積分範囲のyは別物なんですが、stomachmanが両者を同じ記号で書いたから混乱なさったのでしょう。手抜きですいませんね。

>これまでの私の理解では定積分を行う際の積分範囲は0～Y(Yは定数)として考えるものだと考えていました．

その通りです。Yは定数だと思って計算し、それからYを変数yに置き換えれば良いんです。

> tを変数として微分方程式を定積分していく方法に関して，参考になる文献

微分方程式を解くわけじゃなく、単に積分するだけなんですから、高校の教科書で充分でしょう。

> もう一度0～yの定積分を行う際は∫Fdyと∫xdyは定数となっているためにyを掛けるということになると考えます．

　違います。
　積分変数と積分範囲の変数を変えて書くなら、∫Fdy（積分範囲は0～Y)はYの関数です。これを
えふ(Y)=∫Fdy（積分範囲は0～Y)
と書くことにする（勿論、数値積分で、Yがyの各サンプル点の場合について計算します）。これをさらにYで積分すると、
げふ(わい）=∫えふ(Y)dY（積分範囲は0～わい）
これはわいの関数です（勿論、数値積分で、わいがyの各サンプル点の場合について計算します）。
　同様に、
X(Y）=∫xdy（積分範囲は0～Y）
えくす(わい）=∫X(Y)dy（積分範囲は0～わい）
も計算して、モデル
えふ(わい)=Mx(わい)+BX(わい)+Kえくす(わい)+Pわい+Q+誤差(わい)
を作るんです。（わいがyの各サンプル点である場合について。つまり、yのサンプル点の個数だけ、式ができることになります。）

　定積分の範囲yを積分変数yと同じ記号で書いた理由がお分かりになったでしょ？

補足

回答ありがとうございます。確かに式上では定積分の範囲yを積分変数yと同じ記号で書いた理由が理解できました。

何度も申し訳ありませんが，離散化した場合について一つ確認事項があります。

X(Y）=∫xdy（積分範囲は0～Y）
えくす(わい）=∫X(Y)dy（積分範囲は0～わい）
に関して，台形則の数値積分を行うと
X(Y）={(x1+2*x2+x3)*dy/4}+{(x1+2*x2+x3)*dy/4}+・・・・・・
えくす(わい）={(x1+2*x2+x3)*dy^2/4}+{(x2+2*x3+x4)*dy^2/4}+・・・・・・
となります。
つまり，X(Y）= えくす(わい）*dy
が成り立つ式になり，「X(Y）」と 「えくす(わい）」の関係が独立ではなくなってしまうと思われますが，問題ないでしょうか。
どうぞよろしくお願い致します。

回答No.2

usokoku

さいしょう二乗法全般として
最適解の精度を上げるには、実験値の精度を上げること。ホワイトノイズを可能な限り除いてください。ただ、１桁精度を上げるのに最低で実験回数が10倍になります。nを無限大にすれば、積分されて誤差が減りますので。

真空中で空気抵抗を除くとか、振動物と押さえているものの間の接触抵抗を減らしたりする(振り子のアルコール脱脂による手垢の除去等。手垢がついていると蒸発して変な力がかかることがある)ことは、やっていますね。

単に計算精度を上げる、ざん差を減らす、というのであれば、発散しにくい式を使ってください。二乗法ではなく絶対値法をつかうとか、一次式ならばOR法を使うとか。
計算機Eとオーバーフローとけたおちに注意してください。収束が悪くておかしいな、と調べたら(１週間かけて手計算してみたら)桁おちになっていた事が何階かあります。単精度と倍精度を比較してみる・倍精度と４倍制度を比較してみるなどを行ってみてください。極端な場合には、自分で多倍精度塩山ルーチンを用意する必要があります(昔のFM7で24bit整数演算に変換して計算したことがあります。整数演算でしたらばメモリー制限以外の制限無しの可変長のライブラリ(一万桁１０進演算とか)を自作しています)。
1000回の繰り返しがあるなんて場合には、２桁の測定精度でも、計算に必要な有効桁は6桁になります。単精度演算ですと５－６桁目がおかしなことになってくる数値になります。

補足

ありがとうございます。
最小二乗法において解析精度をあげるためには実験方法，装置や解析法を全体的に再検討する必要がありそうですね。

再検討してみたいと思います。

回答No.1

tono-todo

私の知識が乏しいのか、質問の意味が分かりません。

「同定」とは何？

「一般化最小二乗法」とは何？最小二乗法と何が違う？

「式誤差の影響」式には誤差が無いと思うが、あるのは測定値だが？

最小二乗法なら簡単です。
Ｆが従属変数で、３つの独立変数（変位、速度、加速度）の一次結合で表される、
従属変数も独立変数もすべて測定可能で測定されている。
この一次結合式を求めれば良いのではないですか。
教科書にいくらでもあります。

補足

説明不足ですみません。
F，dx^2/dy^2，dx/dyとxが既知で，最小二乗法でM,B,Kを算出することを同定としています。

理解が乏しいのですが・・・・
私の解釈では最小二乗法において，残差の二乗和が小さくなるようにパラメータ，及び推定値を算出しますが，残差が白色でないことから算出する推定値に影響（参考文献「信号解析とシステム同定」ではかたよりと記載されていました。）を与えてしまうらしいです。
そこで，一般化最小二乗法は残差を白色化し，パラメータがより精度良く求める手法のようです。

通常用いられている最小二乗法はすでに試していますが，より精度よくパラメータを算出する方法を探しています。