修士論文 要旨

リスク鋭敏型分布を用いた状態空間モデルのオンラインパラメータ推定
数理情報学専攻 48166233
指導教員

翁 啓翔
駒木 文保 教授

1 はじめにリスク鋭敏型フィルタを用いた状態空間モデルにおけるパラメータのオンライン推定方法について考える．
リスク鋭敏型フィルタの概要について紹介し，どのようにパラメータ推定に応用するかについて述べる．以下
では xt , yt は次のモデルに従うものとする：

x1 ∼ fθ (· | ∅) and xt | (xt−1 = x) ∼ fθ (· | x) ,yt | (xt = x) ∼ gθ (· | x) .
図 1. カルマンフィルタとリスク鋭敏型フィルタの推

2 既存研究定誤差の比較

2.2

2.1 リスク鋭敏型フィルタ

対数周辺尤度最大化によるオンライン推定

線形ガウスモデルにおける状態変数を推定する代表

次に y1 , y2 , · · · が逐次的に観測される場合に θ をオン的な手法であるカルマンフィルタは真の状態変数との 2ラインで推定することを考える．Doucet & Tadic [3]，乗誤差の期待値を最小化するものとして書き表される．

Poyiadjis et al.[5] では粒子フィルタを用いて周辺対数

この 2 乗誤差を一般化することで得られるリスク鋭敏尤度 log pθ̂ (y1:t ) を近似的にオンラインで最大化するア型フィルタ（Risk-Sensitive Filter）がある．
α > 0 をリスク鋭敏パラメータと呼び，パラメータ α(α)のもとで，観測 y1:s を用いた xt を x̂t|sとおくことにする．Vk を正定値対称行列として，リスク鋭敏型フィルタは以下のように定義される：

ルゴリズムを提案した．各時刻での推定値 θ̂t に対して(( )αJ(x1:t , x̂1:t−1, ξ) : = expt−1∑αα(α) 2xk − x̂k|k+ ∥xt − ξ∥2VtVk2 k=12pθ̂ (y1:t ) を計算するには O (t) の時間がかかるため，便宜的に次の擬似尤度関数 pθ̂1:t (x1:t , y1:t ) を考える :)pθ̂1:t (x1:t , y1:t ) :=(t∏k=1fθ̂k (xk | xk−1 ) gθ̂k (yk | xk ) .
)新たに xt+1 , yt+1 , θ̂t+1 が与えられるとき，と前の時(α)α)
x̂(t|t: = argmin E [J(x1:t , x̂1:t−1, ξ) | y1:t ]ξ

このようにして求まるリスク鋭敏型フィルタは仮定したシステムの分布が真の分布からずれている時，フィルタによる状態変数の推定誤差は通常のフィルタより刻までに計算された値と t + 1 時刻目に新たに観測された量の積を使って計算されるので，O (1) で更新されることがわかる．次にこの擬似尤度関数のパラメータに関する勾配を考える．x1:t , y1:t の擬似対数尤度関数の勾配はも経験的に小さくなるという性質があることが Boel etal. [2] によって指摘されている．このことを数値的に表∇θ̂1:t log pθ̂1:t (x1:t , y1:t ) =したのが図 1 である．
また，xt の推定値のみならず，以下のようにその分布も定義することが出来る [1]．ただし，γk(α) (xk , ξ) :=)(2exp α2 ∥xk − ξ∥Vk とおく．

+∇θ̂k log gθ̂k (yk | xk )}t {∑∇θ̂k log fθ̂k (xk | xk−1 )k=1と定義する．ここで，フィシャーの等式を用いると，擬似対数周辺尤度の勾配を St とおけばSt = ∇θ̂1:t log pθ̂1:t (y1:t )p(θα) (xt | y1:t−1 ) ∝}∫ {∑∫t ()(α)=∇ log fθ̂k (xk | xk−1 ) + ∇ log gθ̂k (yk | xk )fθ (xt | xt−1 ) pθ (xt−1 | y1:t−1 ) dxt−1 ,k=1()α)× pθ̂1:t (x1:t | y1:t ) dx1:tp(θα) (xt | y1:t ) ∝ γt(α) xt , x̂(t|tgθ (yt | xt ) p(θα) (xt | y1:t−1 )リスク鋭敏型分布を用いた状態空間モデルのオンラインパラメータ推定 (翁 啓翔)となる．この関数を擬似スコア関数と呼ぶこととする．
通常の最急降下法では ∇ log pθ̂ (y1:t ) を更新方向としてパラメータを更新するが，オンライン推定では∇ log pθ̂ (yt | y1:t−1 ) = ∇ log pθ̂ (y1:t ) − ∇ log pθ̂ (y1:t−1 )が平均的に ∇ log pθ̂ (y1:t ) /t になることから，これを確率的勾配とみなして更新方向とする．以上から次のアルゴリズムが得られる:

θ̂t+1 = θ̂t + ηt ∇θ̂1:t log pθ̂1:t (yt | y1:t−1 )= θ̂t + ηt {St − St−1 } .

図 2. 線形ガウスモデル, 初期値 ϕ0 = 0.1, 真の値ϕ∗ = 0.85 のときの収束の様子

3 本研究で提案するアルゴリズム
既存手法では (x1:t , y1:t ) の対数尤度をフィルタ分布によって積分したものを更新方向としてパラメータを勾配法によって推定していた．積分に用いるフィルタ分布とは更新中のパラメータ θ̂t を用いて与えられるものである．仮に真のパラメータ θ ∗ をプラグインした時のフィルタ分布が pθ∗ (x1:t | y1:t ) が与えられたとして，この分布を用いて上記の更新方向を計算すると既存手法に比べ，早く真の値に収束することが観測された．
本研究ではリスク鋭敏型フィルタから誘導される誘導されるリスク鋭敏型予測分布 p(α) (x1:t | y1:t ) が従来
θ̂t用いられた用いられた分布 pθ̂t (x1:t | y1:t ) よりも経験的に真のフィルタ分布に近いことを用いて，リスク鋭敏型予測分布を用いて期待値計算を行うアルゴリズムを提案した：

図 3. 線形ガウスモデル，初期値 ϕ0 = 0.1，真の値

ϕ∗ = 0.9 のとき，50 回系列を発生させた場合の各系列における誤差．縦軸：提案手法の誤差,横軸：既存手法の誤差θ̂t+1 = θ̂t + ηt d(θ̂α) (y1:t ) ,1:t(α)

d(θ̂α) (y1:t ) := Sθ̂ (y1:t ) − Sθ̂(α) (y1:t−1 ) ,1:t1:t−1∫ 1:tSθ̂(α) (y1:t ) := ∇θ̂1:t log pθ̂1:t (x1:t , y1:t )1:t× pθ̂ (x1:t | y1:t ) dx1:t(α)

1:tとパラメータを更新する手法を提案した．また，数値的には Poyiadjis et al. [5] によるパラメータの手法とOrguner et al. [4] のリスク鋭敏型の粒子フィルタを用いた計算を組み合わせることで提案手法を実装した．
数値実験では線形ガウスモデルと確率的ボラティリティモデルにおいてパラメータを推定する実験を行い，初期値によるものの両モデルにおいて既存手法より早く真のパラメータに収束する例があることが確認された（図 2, 3）．また，時系列モデルではないものの，状態変数が独立同一分布に従う潜在変数モデルの単純な場合において提案アルゴリズムが従来手法よりも早く真のパラメータに収束する解釈性のある理由を与えた．

参考文献[1] S. Bhaumik, S. Sadhu, and T. Ghosal, “Alternativeformulation of risk-sensitie particle filter,” in Proceedings of 2006 Annual IEEE India Conference, 2006,
pp. 1–4.
[2] R. Boel, M. James, and I. Petersen, “Robustness andrisk-sensitive filtering,” IEEE Transactions on Automatic Control, vol. 47, pp. 451–461, 2002.
[3] A. Doucet and V. Tadic, “Parameter estimation ingeneral state-space models using particle methods,”
Annals of the Institute of Statistical Mathematics,
vol. 55, pp. 409–422, 2003.
[4] U. Orguner and F. Gustafsson, “Risk-sensitive particle filters for mitigating sample impoverishment,”
IEEE Transactions on Signal Processing, vol. 56, pp.
5001–5012, 2008.
[5] G. Poyiadjis, A. Doucet, and S. Singh, “Particle approximations of the score and observed information
matrix in state space models with application to parameter estimation,” Biometrika, vol. 98, pp. 65–80,
2011.