読者です 読者をやめる 読者になる 読者になる

粗大メモ置き場

個人用,たまーに来訪者を意識する雑記メモ

RLS(逐次最小二乗法) program for MATLAB

メモ メモ-MATLAB コード コード-MATLAB

2016/7/4 追記

1入力1出力の系での推定は次の記事プログラムを参考にしてください。
ossyaritoori.hatenablog.com

背景:MATLABで使いやすい逐次最小二乗法のコードがない

僕が無知なだけの可能性は高いですが,
MATLABで逐次最小二乗法をしたいなってなった時にいい感じのコードがないんですよね。

組みました

はい。今回は関数ファイルとクラスファイルの2通りで実装しました。
MATLABオブジェクト指向の書き方をしたのは初めてで友人にいろいろ教わりながらやりました。

逐次最小二乗法は結構書き方にいろいろ流派があるように思いますが,
今回は簡単のため以下のURLのPDFに準拠しています。
http://yoh.ise.ibaraki.ac.jp/edu/graduateschool/Project1.pdf

ある程度原理を理解すれば関数自体を組むのは簡単ですが
グローバル,あるいはローカルで逐次変化する変数を保存する必要があり,
コードがバラけそうなのでクラスファイルで定義しています。

ざっくりした説明?

ここでは次のような線形モデルを仮定します。
f:id:ossyaritoori:20160628235418p:plain
ytが出力,ztが過去の入力と出力で,θがその係数です。
vは雑音を示しています。

逐次最小二乗法の目的は,ざっくり言うと
「ytとztを与えてθを推定する」ことです。

具体的な式

推定式の導出はやや面倒で結果を示すと以下のような式で更新をしていきます。
f:id:ossyaritoori:20160628234854p:plain
ただし,ここでは
f:id:ossyaritoori:20160628234924p:plain
としています。


実質,現在の出力:Yn と
過去の参照入出力 Zn' = [Yn-1,Yn-2,...,Un-1,...,Un-m] (2M×1のベクトル)
がわかれば現在の推定値 θnを更新していけるということになります。

その他、忘却係数について

初期値は今回結構雑でもいいのですが忘却係数の扱いを少し特殊にしました。
忘却係数は0~1の間の数で、1に近いほど過去の結果を尊重した慎重な検出になり、
1より小さくしていくとパラメータの時間変化に強くなりますがノイズに弱くなります。

本プログラムでは最初は忘却係数を小さめにすぐに結果をアップデートするようにし、
時間がたつにつれ1に近づけることで定常地付近でのノイズへの耐性をつけています。
f:id:ossyaritoori:20160629000350p:plain


ふぅ。文章打つほうが大変ですね。

追記>
なお、こんな感じに時間で変数の推定が収束しているのがわかります。
ノイズが入ってもきちんと推定できます。

f:id:ossyaritoori:20160629005058p:plain

ソースコード (関数版)

ソースコードは今回Gistというサービスを使ってみました。
はてなブログでソースコードを貼り付ける方法を調べてみた - ウェブと食べ物と趣味のこと

毎度毎度長いのを貼ってしまうのでこれからは長い奴はこういう形式にしようかと


逐次的に保存すべき変数:共分散行列Pnとθn
入力:Yn,Zn

最初の行を削除すれば忘却係数Rhoは消してもいいかと思います。

ソースコード (クラスファイル版)

クラスで定義したソースコードがこちら。
オブジェクトが逐次変化する変数を保存してくれるので実行時のコードが格段に短くなります。

上が定義したクラスファイルで
下がテスト用のファイルになってます。

珍しくいろいろコメント入れましたが,やっぱりあんまり読みやすくはないですね。
コメントや修正などありましたら受け付けます。(聞くだけw)