#freeze
*同時推定の考え方 [#ac3166df]
線形離散時間モデルのパラメータと状態を同時に推定する方法を解説する。
尤度方程式を解くとパラメータに関する線型方程式と状態に関する線形方程式が得られる。この方法は２つの方程式を交互に繰り返して推定する。
-本方法は、逐次解法であるので、オンラインの適応制御に有用である。
-1入力1出力系で示しているが、多入力多出力の高次系にも、拡張できる。
ここで提案した方法は漸近的に有効な推定値を与える。
--漸近有効性とは、ある母数θ に対し、漸近正規推定量のうち分散nが最小のものを、漸近有効推定量と呼び、この良好な性質をもつ推定方法である。
*線形離散時間モデル [#e1262ae6]
入力utと出力ytの観測データに基づいて、離散時間モデルを逐次推定する方法を紹介する。
線形動的システムのモデルとして、次の最も簡単な1次モデルをまず考えよう。
xt+1=fxt+gut+vt 状態推移式
yt=xt+wt 観測式
-xtが状態であり、vtはシステムノイズ、wtは観測ノイズである。
-未知パラメータは(f,g)で定数である。
-ノイズは、それぞれ独立な正規確率変数であり、互いに無相関であると仮定する。
E(vt)=0,E(vtvj)=σv^2・δtj
E(wt)=0,E(wtwj)=σw^2・δtj
E(vt,wj)=0,E(x0,
-初期状態x0は、正規分布N(μx0,σ0)に従い、ノイズと無相関とする。
E(x0)=μx0 E((x0-μx0)^2)=σ0
*尤度関数 [#b3f54393]
パラメータが所与のもとで、n期までの観測データの確率密度関数を次のように表記し、この確率密度関数を最大とする推定値を求める。パラメータベクトルをΘ=(f,g)'とする。
Ln(Θ)=P(YDn,XDn|Θ)
パラメータベクトルをΘ=(f,g)'とする。
YDnはn期までの出力データ：YDn={yt}t=1,~n
XDnはn期までの入力データ：xDn={xt}t=1,~n
ここで、n期までのデータに基づく推定値を次のように表記する。
f*n,g*n:n期までのデータに基づくｆとgの推定値。
Θn*=(f*n,g*n):n期までのデータに基づくパラメータの推定値
x*k|n :n期までのデータに基づく事後的な状態推定値。
--状態推定値x*k|nは、k<nのとき内挿値と呼び、k=nのときフィルタリング値と呼ぶ。
*簡単な試算 [#ue4c41a9]
t期の状態変数が、初期状態とt期までの入力とシステムノイズから、正規確率密度関数であらわされることを試算してみよう。
x1=fx0+gu0+v0
x2=fx1+gu1+v1=f^2x0+fgu0+gu1+fv0+v1
x3=fx2+gu2+v2=f^3x0+f^2gu0+fgu1+gu2+f^2v0+fv1+v2
=f^3x0+(f^2gu0+fgu1+gu2)+(f^2v0+fv1+v2)
.....
より
xt=f^tx0+ Σf^(t-1-i)ui+ Σf^(t-i)vi t=0,1,2,...
Σ は i=0からt-1までの合計
で表される。
これは、確率変数x0,viの加法和であり、互いに無相関で、それぞれ独立な正規分布であった。
uiは既知の入力系列である。そこで、t期の状態変数が正規確率密度関数で表されることがわかる。
また、観測地も観測方程式より、正規確率密度関数で表されることがわかる。
*尤度関数の計算 [#hd411cb3]
逐次推定を行うので、尤度関数を書き直して、正規分布の積の形に表せることを示す。
ベイズの定理より、尤度関数の漸化式がえられる。
Ln(Θ)=P(YDn,XDn|Θ)
=p(yn,xn|Θ,XDn-1,YDn-1)P(XDn-1,YDn-1|Θ)
=P(yn|xn)P(xn|Θ,xn-1)Ln-1(Θ)
尤度関数のn=0での同時確率は、次式の通り。
L0(Θ)=p(yo,x0|Θ)=P(y0|x0)p(x0)
尤度関数の第1項の計算は、観測式yt=xt+wtより、xnが所与のときのynの確率密度関数は
P(yn|xn)=N(xn,σw^2) n=0,1,2,...
で表される。
尤度関数の第2項の計算は,
P(xn|Θ,xn-1)=N(fxn-1+gun-1,σv^2)
これらを用いて、尤度関数の漸化式から、最終的な尤度関数は、正規分布の積で表され
Ln(Θ)=P(yn|xn)P(xn|Θ,xn-1)Ln-1(Θ)
=N(xn,σw^2)N(fxn-1+gun-1,σv^2)Ln-1(Θ)
これより、''求める尤度関数は以下の正規分布''である。
Ln(Θ)=(2πσw)^(-n-1)・σv＾(-n)・σ0^(-1)
EXP{-Σ[(yi-xi)^2/(2σw^2)]-Σ[(xi-fxi-gui)^2/(2σv^2)]-(xo-μx0)^2/(2σ0^2)｝
ただし、最初のΣはi=0~nの合計で、次のΣはi=1~nの合計
*パラメータと状態の同時推定 [#g081990a]
上記の尤度関数を最大とする状態xiとパラメータΘを求めればよい。
そこで、尤度関数の対数をとり、これをパラメータと状態に関して一回微分をとりゼロと置く。
-初期状態で微分
(y0-x0)/σw^2-(x0-μx0)/σ0^2 + f(x1-fx0-gu0)/σv^2=0
-xk(k=1~n-1)で微分
(yk-xk)/σw^2-(xk-fxk-1-guk-1)/σv^2 + f(xk+1-fxk-guk)/σv^2=0
-xnで微分
(yn-xn)/σw^2-(xn-fxn-1-gun-1)/σv^2=0
-パラメータf で微分
Σ {xi-1(xi-fxi-1-gui-1)/σv^2}=0 Σはi=1~nの合計
-パラメータg で微分
Σ {ui-1(xi-fxi-1-gui-1)/σv^2}=0 Σはi=1~nの合計
上記の尤度方程式を連立して解くことで、最尤推定できる。
しかしながら、全体としては非線形であり収束計算が必要となる。
-状態で微分した上側の3式については、x0,x1,...xnに関して線形である。
-パラメータf,gで微分した下側の2式については、f,gに関して線形である。
そこで、次の２つの問題に分けて、相互依存的な線形問題として解くことを考える。
*問題１：状態推定問題：パラメータ所与 [#n3dcb076]
パラメータ所与のもとでの状態の推定（平滑：スムージング）問題である。
-状態で微分した上側の3式をxtの線形関数として整理してみよう。
-第1の式は
(y0-x0)/σw^2-(x0-μx0)/σ0^2 + f(x1-fx0-gu0)/σv^2=0
より
box0+b2x1=q0
ただし
b0=1+(σw^2/σ0^2)+f^2(σw^2/σv^2)
b2=-f(σw^2/σv^2)
q0=y0+μx0(σw^2/σ0^2)+fgu0(σw^2/σv^2)
となる。
-第2の式は
b2xk-1+b1xk+b2xk+1=qk
ただし
b2:同じ
b1=1+(1+f^2)(σw^2/σv^2)
qk=yk+(uk-1-fuk-1)g(σw^2/σv^2)
-第3の式は
b2xn-1+b3xn=qn
b2:同じ
b3=1+(1+f^2)(σw^2/σv^2)
qn=yn+un-1g(σw^2/σv^2)
以上より、''パラメータからb0,b1,b2を与え、q0,qk(k=1~n-1),qnを観測データから計算しておけば
逐次内挿値Xk,k=1~nを求めることができる。''
*問題２：パラメータ推定問題：状態推定値所与 [#j16afe57]
最尤方程式の4式と5式を使って、パラメータを求めよう。
この場合内挿値xkが問題1を解いて与えられているものとする。
n個の状態内挿値の内積を定義する。
x'x=Σxi^2 Σはi=0~n-1の合計
x'u=Σxiui Σはi=0~n-1の合計
u'u=Σui^2 Σはi=0~n-1の合計
これから、下記の2つの式の形に、5式,6式が整理できる。
Σ(xi^2)f+Σ(xiui)g=Σ(xi+1xi)
Σ(xiui)f+Σ(ui^2)g=Σ(ui+1xi)
f,gは、上記の2個の線形連立方程式を解けばよい。
*逐次推定 [#w6795385]
すべての過去の状態の内挿値が、各期で必要になることが、難点であろう。
そこで、次の繰り返し計算を提案する。
-1.パラメータΘ=(f,g)の初期推定値を適当に与え、問題1で状態の推定を行う。
-2.求められた状態推定値から、パラメータの再推定を行う。
-3．上記を繰り返して、次の収束判定を行う。
||Θ*m-Θ*m-1||<ε εは十分小さい値。（例えば、0.001)
良好な収束性を示したとの報告が、研究論文にある。
パラメータが、十分に収束した時点では、通常のカルマンフィルタ を用いた推定値で、問題1の最適内挿値の代用をしても良いであろう。
*参考 [#seeac7f4]
-y.BAR-SHAROM:"Optimal Simulations State Estimation and Parameter Identification in Linear Discrete-Time Systems",IEEE AC-17,no.3,JUNE(1972)
![[PukiWiki]](image/pukiwiki.png "[PukiWiki]")
