空間斉次自己回帰モデルの乱数生成とそれに基づく実験慶應義塾大学、早稲田大学力丸佑紀（株）データサイエンスコンソーシアム、慶應義塾大学、早稲田大学柴田里程一般に、弱定常空間過程の対数尤度は厳密に求めることは困難である。歴史的には Whittle(1954) や Guyon(1986), Robinson(2006) などがさまざまな形の近似尤度を提案しているが、いずれも複素関数の積分計算が含まれており実用的ではない。しかし、空間斉次自己回帰モデルに限れば、簡単な計算で求められる近似尤度が存在することを発見し、一昨年の当学会で報告した。理論的にはこの近似尤度を最大化することで（フィッシャー情報量が正則である限り）漸近有効性をもつパラメータ推定量が得られるが、数値的な検証をおこなうためには、空間斉次自己回帰モデルに従う乱数を実際に生成する必要があった。しかし、われわれの調べた限り、その生成法はいまだ確立されていないようである。そこで、本報告では、われわれの考えた空間斉次自己回帰モデルに従う乱数の生成法を紹介するとともに、その乱数を用いたパラメータ推定実験の実験結果を報告する。以下、簡単のため、2 次元空間上の空間斉次自己回帰モデル P (T1 , T2 )Xv = εv , εv ∼ N (0, σ 2 I) (1) に限り、n1 × n2 の大きさの長方格子上で観測値 {Xv } が得られたとする。ただし、弱定常性を保証するため P (z1 , z2 ) は |z1 | = |z2 | = 1 上で根をもたない 2 変数の（負のべき乗も含む）多項式であるとする。時系列の自己回帰モデルと異なり {εv } が独立ではあるもののイノベーションではないため、(1) の回帰式を直接利用した乱数生成法は収束する保証がなく、実際うまく機能しない。また、共分散行列もパラメータの複雑な関数となるため、共分散関数を利用した乱数生成も実用的ではない。そこで、私たちは、伝達関数が P (z1 , z2 ) = P1 (z1 )P2 (z2 ) のように分解できる場合に限れば、次のような乱数生成法が存在することを発見しその有効性を検証した。ただし、P1 (z1 ), P2 (z2 ) は負のべき乗も含む多項式である。具体的には、P1 (z1 ) を P1 (z1 ) = c p1 −1 ∏ ∏ (1 − αj z1 )(1 − αk z1−1 ) （ただし、c は定数） j=1 k=−p2 のように因数分解し、その各項を |αj | < 1 ならば (1 − αj z1 )−1 = ∞ ∑ (αj z1 )ℓ ℓ=0 |αj | > 1 ならば (1 − αj z1 )−1 = −(αj z1 )−1 ∞ ∑ (αj z1 )−ℓ ℓ=0 のように展開する。(1 − αk z1−1 ) については、z1 を z1−1 で置き換えるだけである。同様に P2 (z2 ) についても逆演算をべき級数に展開し、順次、(1) の両辺に施すことにより、左辺の演算子を簡略化し、最終的に Xv = P (T1 , T2 )−1 εv を得る。実際には級数は適当な次数で打ち切り近似することで、 {εv } から {Xv } を作り出す。この方法で生成した乱数のペリオドグラムとスペクトル密度を比較したところ、非常に近い形状となり、もとのモデルに従う乱数を生成できていることが確認できた。また、この乱数を用いた数値実験により n1 = n2 = 30 程度でも近似最尤推定量が十分な精度をもつことが確かめられ、その実用性が検証できた。