最小二乗法の逐次解法 †最小二乗法で代表されるシステムのパラメータ推定手法について、紹介する。特に、時系列データーのように、入出力データーが逐次的に得られる場合のパラメーター推定の方法を示す。 線形離散システムのパラメーター推定 †線形システムの次数n を既知と仮定すると,その入出力関係は次式のような差分方程式で 表される. xt = Σaixt-i +Σbiut-i i=1,n eq1 yt = xt + vt eq2 但し、utは入力、ytは出力 Vtは平均値ゼロの離散値確率雑音。入力と無相関 E[utvt] = 0
ベクトル形式で表現 †2式を1式に代入して整理する。 yt = Σaiyt-i +Σbiut-i + rt 但し rt =Σaivt (a0 = 1) eq3 上式は次のようなベクトル形式で表現できる。 yt = z't θ + rt z't = [-yt-1,....,-yt-n, ut-1,...., ut-n] θ' = [a1,...., an, b1,...., bn] 誤差最小化問題と推定値 †パラメータを推定するため,次のような二乗誤差評価を最小にすることを考える。 データ数はNとしよう。 J(θ) =Σρ^(N-t)rt^2 =Σρ^(N-t)(yt - z't θ)^2 ただし,0 < ρ <= 1 は忘却係数である. 忘却の効果を小さくすれば,ゆっくり変動する時変パラメータの推定にも適用できる。 J を最小にする必要条件はJ(θ)の一階偏微分ベクトル = 0 である.すなわち, ∂J/∂θ= -2Σρ^(N-t)ztyt + 2Σρ^(N-t)ztz't θ = 0 これより、θの最小二乗推定値は θ*=[Σztz't]^(-1)・Σρ^(N-t)ztyt 上式が成り立つためには,上記の逆行列が存在する必要があるが、一般に入力信号が十分多くの周波数成分を含んでいれば成立する。 逐次推定の考え方と方法 †データ数がNの場合のθの最小二乗推定値をθ*(N)と置く。 パラメータ推定のための評価規範を、下記のように置いてみよう。 J(θ;N) =(1/N)Σρ^(N-t)rt^2 =(1/N)Σρ^(N-t)(yt - z't θ)^2 これは、次のように書きかえられる。 J(θ;N) =C(N)-2θ'b(N)+θ'R(N)θ ただし b(N)=(1/N)Σρ^(N-t)ytzt R(N)=(1/N)Σρ^(N-t)ztz't C(N)=(1/N)Σρ^(N-t)(yt)^2 ここでJ(θ;N) をθ に関して微分して0とおく。すると、正規方程式(normal equation) と呼ばれるθに関する連立一次方程式が得られる。この解が推定値である。 R(n)θ*(N)=b(N) :正規方程式
前節で説明した最小二乗法をもとに、再帰的最小二乗法(RLS アルゴリズム)を導出してみよう。 正規方程式は次のように、分解表現できる。 θ*(N)=R(n)^(-1)・b(N) =[Σρ^(N-t)ztz't]^(-1)・(Σρ^(N-t)ytzt) 5式 まず、行列P(N) を P(N)=[Σρ^(N-t)ztz't]^(-1) 6式 とおき、これを共分散行列と呼ぶ。すると、 P(N)^(-1) =Σρ^(N-t)ztz't =P(N-1)^(-1)+zNz'N 7式 が得られる。同様にして b(N)=Σρ^(N-t)ytz't + yNz'N 8式 6から8式を5式に代入して θ*(N)=P(N)・(Σρ^(N-t)ytzt) =P(N)・[P(N-1)^(-1)θ*(N-1)+yNzN] =P(N)・{[P(N)^(-1)-znz'n]θ*(N-1)+yNzN] =θ*(N-1)-P(N)zn[z'Nθ*(N-1)-yN] =θ*(N-1)+P(N)zn[yN-z'Nθ*(N-1)] 9式 9式と7式が、再帰的最小二乗法(RLS アルゴリズム)である。 しかし、P(N)の逆行列 をオンラインで計算することは困難である。そこで、逆行列補題(matrix inversion lemma) を用いて、(1.17)式をオンライン計算が可能な形式に変形しよう。
7式に逆行列補題を適用すると、次式が得られる P(N)=[P(N-1)^(-1)+zNz'N]^(-1) =P(N-1)-P(N-1)(zNz'N)P(N-1)/{1+z'NP(N-1)zN} 10式 9式も右辺第二項に含まれるP(N)zNを上式を代入することで、さらに簡単に修正できる。 P(N)zN={P(N-1)-P(N-1)(zNz'N)P(N-1)/[1+z'NP(N-1)zN]}zN = P(N-1)zN/[1+z'NP(N-1)zN] 11式 上のP(N)zNを9式に代入して整理する。 θ*(N)=θ*(N-1)+{P(N-1)zN/[1+z'NP(N-1)zN]}・[yN-z'Nθ*(N-1)] 12式
推定値の漸近的性質 †データ数が多くなるにつれて,確率極限 N-->無限大の時、plim θ*= θ が成立するとき,θ* は一致推定量であるという。
plim θ*= {plim (1/N)Σztz't}^(-1) plim (1/N)Σzt[z't θ + rt] = θ +{plim (1/N)Σztz't}^(-1)plim (1/N)Σzt rt が得られるので、一致推定量であるためには,ベクトルzt と式誤差rt は無相関でなければならない。 plim (1/N)Σzt rt = 0 が成り立つ必要がある。
補助変数法によるパラメータの一致推定方法 †ベクトルzt と同サイズの補助変数ベクトルmt を導入し,次のように補助変数法による推定 値 θ** を得る。 θ**=[Σmtz't]^(-1)・Σρ^(N-t)mtyt 前節と同様に式変換することで、補助変数ベクトルmt が次のような条件を満足すれば,パラメータ推定値は一致推定量であることが言える。
この2つの条件を満たす補助変数mtの選び方は、いろいろある。
最小二乗法および補助変数法は,データを蓄えておき,データ行列の逆行列を 計算して推定値を得るバッチ処理である.このような方法は繰り返し計算を必要としないこと からオフライン法とも呼ばれている.
逐次推定アルゴリズム †N 個のデータから得られる推定値をθN とする.そこで,θN を過去の推定値及び時刻N でのデータより求めることを考える.新しくデータが得られる度に直前の推定値を修正していく方式(逐次計算式)が実現されれば,オンライン推定や実時間推定への適用が可能である. 最小二乗法と補助変数法に関する逐次推定アルゴリズムは以下のようになる. ただし,PN は2n x 2n の行列で,LN は2n x 1 のベクトルである. 上式において,逐次最小二乗法は次のように与えられる. また,逐次補助変数法は次のように与えられる.
練習問題:2次の時変係数の線形系 †2 次の時変システムを考える。逐次最小二乗法で上記のシステムの時変パラメータを推定せよ.時刻を横軸として,各推定パラメータの挙動を,パラメータの真値とあわせて図示せよ.忘却係数及び観測雑音によるパラメータ推定値への影響を考察せよ. 参考になる文献 †
|