2016/7/4 追記
1入力1出力の系での推定は次の記事プログラムを参考にしてください。
ossyaritoori.hatenablog.com
組みました
はい。今回は関数ファイルとクラスファイルの2通りで実装しました。
MATLABでオブジェクト指向の書き方をしたのは初めてで友人にいろいろ教わりながらやりました。
逐次最小二乗法は結構書き方にいろいろ流派があるように思いますが,
今回は簡単のため以下のURLのPDFに準拠しています。
http://yoh.ise.ibaraki.ac.jp/edu/graduateschool/Project1.pdf
ある程度原理を理解すれば関数自体を組むのは簡単ですが
グローバル,あるいはローカルで逐次変化する変数を保存する必要があり,
コードがバラけそうなのでクラスファイルで定義しています。
ざっくりした説明?
ここでは次のような線形モデルを仮定します。
ytが出力,ztが過去の入力と出力で,θがその係数です。
vは雑音を示しています。
逐次最小二乗法の目的は,ざっくり言うと
「ytとztを与えてθを推定する」ことです。
具体的な式
推定式の導出はやや面倒で結果を示すと以下のような式で更新をしていきます。
ただし,ここでは
としています。
実質,現在の出力:Yn と
過去の参照入出力 Zn' = [Yn-1,Yn-2,...,Un-1,...,Un-m] (2M×1のベクトル)
がわかれば現在の推定値 θnを更新していけるということになります。
その他、忘却係数について
初期値は今回結構雑でもいいのですが忘却係数の扱いを少し特殊にしました。
忘却係数は0~1の間の数で、1に近いほど過去の結果を尊重した慎重な検出になり、
1より小さくしていくとパラメータの時間変化に強くなりますがノイズに弱くなります。
本プログラムでは最初は忘却係数を小さめにすぐに結果をアップデートするようにし、
時間がたつにつれ1に近づけることで定常地付近でのノイズへの耐性をつけています。
ふぅ。文章打つほうが大変ですね。
追記>
なお、こんな感じに時間で変数の推定が収束しているのがわかります。
ノイズが入ってもきちんと推定できます。
ソースコード (関数版)
ソースコードは今回Gistというサービスを使ってみました。
はてなブログでソースコードを貼り付ける方法を調べてみた - ウェブと食べ物と趣味のこと
毎度毎度長いのを貼ってしまうのでこれからは長い奴はこういう形式にしようかと
逐次的に保存すべき変数:共分散行列Pnとθn
入力:Yn,Zn
最初の行を削除すれば忘却係数Rhoは消してもいいかと思います。