論文・研究レポート

生態学的予測モデルに対するカルマンフィルタによる近似法の研究
上田徹，中村剛

川…llllll …‖‖‖‖‖＝‖‖‖＝‖‖＝‖＝‖‖州………ll llIl‖‖‖＝＝＝＝‖‖‖州 I…＝＝‖‖川LLLl＝＝川LLLllL…＝＝IuHLLHL LLLLlllLLl川 LLL（（（IulIMll…11u MlutLILuIlu LLLLlllLL l… uM‖＝tIu M…＝……………Il州Il＝＝It Lttlt…l川＝は状態Ⅹ（乃）を逐次的に推定するアルゴリズムである．

1．まえがきトレンド成分，周期成分に関してはカルマンフィルタ［1，5］は予測誤差を調整しながら∇27、（刀）予測を進めるため短期間の予測では大きな誤差を避け…（T（乃）一丁（乃−1））−（r（乃−1）一丁（乃−2））ることができる．また，ロジスティ ックモデルや＝〟（乃）Bassモデルのような生態学的予測モデル［2］はモデル

ぶ−1の係数に定性的な意味があり，長期の構造を吉見明する

ど＝0

（3）

∑S（乃−オ）＝〃（乃）（4）が成り立つものとする．ただし，α（乃），〃（乃）は互いにのに適しているが，短期的には大きな誤差が生じる可独立かつ平均0で，それぞれ分散r2，♂2（未知）の正台㌢性がある．そこで，これらの方法の長所だけを組み規性ホワイトノイズと仮定し，周期成分の周期は5と合わせることができないかを検討する．すなわち，カする．式（3）のトレンド成分表現を滑らかなトレンドモ

ルマンフィルタを適用できるモデルは基本的には線形デルと呼ぶ．このとき，モデルであるが，カルマンフィルタを非線形モデルのⅩ（乃）＝（T（乃），T（乃−1），5（乃），…，5（乃一ざ＋2））f近似として適用できないかを検討する．そして，2種（5）の競合関係を考慮したモデル［3，4，6］にカルマンフと置くと式（1），（2）のF（乃），C（5），〟（乃），u（乃）およぴ

ィルタを適用できるかも検討する．

U（乃）の分散共分散行列￠（乃）は2 −1 0 0 … 0 02．カルマンフィルタの概要

1 0 0 0 … 0 0

時点乃でのトレンドや季節性などに関する状態を表す変数Ⅹ（乃）を推定・予測したいものとする．時点
F＝

乃から時点（乃十1）への状態変化は，雑音U（乃）を考慮して構造（システム）方程式
Ⅹ（わ＋1）＝F（乃）Ⅹ（乃）＋G（乃）u（乃）0 0 −1−1‥・

一1−1

0 0 1 0 …

0 00 0 0 1 …

0 00 0 0 0 … 1 0（1）で与えられるものとする［乃＝1，2，…，〃］．時点乃で1 0の観測値y（循）は，トレンド成分r（乃）と周期成分S（乃）と雑音紬（乃）【分散山2］の和，すなわう観測方程C＝式

0 0〟（乃）＝（1，0，1，0，…，0）亡0 1u（乃）＝（祝（乃），む（乃））f（6）0 0y（乃）＝r（乃）＋S（乃）＋抑（乃）＝〟（乃）f∬（乃）＋抑（乃）輔▲ご2］（2）で表されるものとする（式（5）のⅩ（乃），式（6）の〟（耽）参照）．このように構造方程式と観測方程式で表現されるモデルを状態空間表現といい，カルマンフィルタ0 0である．推定すべきパラメタβ＝（〟2，r2，♂2）のうち〟2は1と置くことができ，その他は最尤法を用いて推定できる．

うえだ とおる，なかむら つよし（現在，NEC）
成渓大学工学研究科〒180−8633武蔵野市吉祥寺北町3−3−1
受付02．4．8 採択02．11．7
Tg4（42）時点乃2までの観測値を用いたz（乃1）の推定値をz（乃1l乃2）とすると，（順方向では）© 日本オペレーションズ・リサーチ学会. 無断複写・複製・転載を禁ず.

オペレーションズ・リサーチ

【時間更新】r（乃）一丁（乃−1）＝α2r（ヰー禦）y（乃l乃−1）＝〃（乃） x（乃l乃−1）（1カ式（10）と同様に斤（乃l雅一1）…1ねγ【y（乃l乃−1）】。腑＝1−＝〝（乃）や（乃l乃−1）〃（乃）＋ひ2
x（乃l乃−1）＝Fx（乃−1l乃一1）

ヱ也
3‖監椚と置き換えてP（乃l和一1）…Vαγ【Ⅹ（乃l乃−1）］r（乃−1）＝（1−α2C桝l）r（乃）＝ダP（乃−1l乃−1）F ＋C￠C

（7）

【観測値による更新】（14）とする．
式（ll），（14）に誤差項を加えたものを，式（3）および逆向Ⅹ（乃l乃）＝Ⅹ（乃l乃−1）＋瓜（y（乃ト〃（乃） x（乃l光一1））きカルマンフィルタのトレンド項表現の代わりに用い

P（乃l乃）＝P（乃l雅一1）−瓜〃（乃） P（乃l乃−1））

る．すなわち，式（6）における行列ダ（乃）の左上の2×2行列の部分を通じて予測値y（刀l雅一1）を逐次，得ることができる．

勅）＝［∴1］

ただし，lねγ（ズ）はズの分散を意味する．ここで，卜を，順方向のカルマンフィルタでは乃に依存する﹈
0

乃C0

向きに走らせることにより推定できる．逆向きカルマ

α1

十l

初期値Ⅹ（Olo）はカノ㌧マンフィルタを最終時点から逆1ンフィルタでもF，C，〃，uは変わらないが，状態ベ﹈

01

乃の卜に置き換え，逆方向のカルマンフィルタではクトルは

C＋01

Ⅹ（乃）＝（r（乃），r（乃＋1），S（刀），・‥，S（乃＋5−2））としなければならない．に置き換える．

3．各種モデルのカルマンフィルタによる

3．2 Bassモデルl
新製品拡散モテリレとして知られているBassモデル

近似［2，3］は，研を総（最終）需要とするとここでは生態学的予測モデルを用いてトレンド項

dr（′）励＝（椚一丁（′））（仇＋（β．／∽）r（り）

r（乃）を近似することを考える．と表現される．ロジスティックモデルと同様に左辺の3．1ロジスティックモデル［2］微分を差分で表現し，右辺の時点調整を行うと需要予測によく使われるモデルであるロジスティッT（乃＋1）一丁（乃）クモテリレは次式で表される．

（19＝（椚一丁（乃一1））（仇＋（βl／∽）r（乃））肇粁仇（1一望），研：総（最終）需要＝となる．ここで上式の微分を差分で表現し，時点乃での増加率は時

c乃＿1＝∽一丁（乃−1）（1￠と置き換えると点（乃−1）で決まるとすると，

T（乃＋1）＝（1十（β1／椚）c〃＿1）r（乃）＋α1C乃＿1

r（州トr（乃）＝のr（ヰーヱ怒吐）（9）となる．とする近似が考えられる．そして次に逆向きカルマンフィルタのシステムもc乃−．＝1一T（乃）一丁（乃−1）（10）＝（∽一丁（乃＋1））（勉＋（戯／研）r（乃））と置きかえると（用と変更し，r（乃＋1）＝（仇C〃＿1＋1）r（乃）（11）r（乃−1）＝（1−（戯／椚）c乃吊）r（乃）一助C侶lと表すことができる．時点乃でc乃＿1を確定値として扱うとT（乃）を定数倍したものがr（乃＋1）というこ（1幼とする．
式（1れ（用を式（11），（14）の代わりに用いる．すなわち，とになる．

﹈0︶

〃 C＋0

伽1

凡（乃）の代わりに，順方向のカルマンフィルタでは＋l

逆向きカルマンフィルタ用のシステムを変更する．式ト

次に，カルマンフィルタで初期値を求めるために，（9）と同様に2002年12月号（1乃© 日本オペレーションズ・リサーチ学会. 無断複写・複製・転載を禁ず.

（43）丁95

を用い，逆方向のカルマンフィルタでは﹈0

れ（乃十1）＝れ（乃）（1＋ん乃＿1）

01−（＆／椚）c乃吊と置くと，l

右（乃＋1）＝右（乃）（1＋ん，乃＿1）［となる．そこで，状態ベクトルをを用いる．

Ⅹ（乃）＝［Ⅹ．（乃），Ⅹ2（乃）］また，式（lれ（1鋸こおける右辺第2項のために，

Ⅹf（乃）＝【r（乃），れ（乃−1），Sf（乃），…，5f（乃−∫十2）］

Ⅹ（乃＋1）＝F（乃）Ⅹ（乃）＋C（乃）＋C（乃）U（乃）C（乃）f＝（広，0，0，…，0）のように節3．3までの状態ベクトルを二つ繋げ，であり，式（17）ではF（乃）もダ（乃）ヂこ乃）凡て乃，］（九＝仇C乃＿1
式（lq）では1＋ん乃＿10 0 0 … 0 0（九＝一助C桝1である．
3．3 BassモデルIl

式（川の右辺第2項はr（乃−1）の1次式なのでその占（乃）＝まま用いることとし，

10 0 0 ‥・0 000 −1−1…

−1−100 1 0 … 0 00 0 1 … 0 000 0 0 … 1 00

r（刀＋1）＝（1＋（β1／椚）c乃＿1）r（乃）＋仇（研一r（乃−1））餉と置く．とし，他の行列，ベクトルも行，列の長さを同様に倍逆向きカルマンフィルタもシステムを変更し，長にすればよい．なお，0は占（乃）と同じ大きさで要素がすべて0の行列である．

r（刀−1）＝（1＋（＆／椚）c桝．）r（乃）十勉（研一r（乃＋1））帥

4．モデル評価と置く．
式帥，飢）を式（lれ（用の代わりに用いる．すなわち，4．1評価に用いるデータと尺度

j㌔（乃）の代わりに，順方向のカルマンフィルタでは

本節では総務庁統計局作成の「家計調査年報」の「1世帯当たり1ヶ月の品目別購入数量（全世帯）」の酒類項目における「ビール」と「清酒」のデータを用［1＋干β1㌣）c乃￣1￣ご1］を用い，逆方向のカルマンフィルタではいる．1983−1997年のデータをパラメタ推定用に用い，1998年のデータを予測誤差評価用に用いる（乃＝1￣（β2）c乃￣1

1：1983年1月，乃＝Ⅳ…180：1997年12月）

冨］［を用いる．

二乗誤差和が最小となるようにパラノタ推定を行っ

3．4 競合を考慮したモデル

た．例えばロジスティックモデルでは，競合の起きる2種の観測値をれ（乃），右（乃）とする∈（乃）＝r（乃＋1トr（乃トのT（乃）（1−r（乃）／椚）と，競争的2種モデルはとして，1983−1997年のデータに対してdれ値＝れ（α−∂n−C右）JV−1

且0＝∑ど（乃）2
乃＝1

d右ノ肋＝右（d−eれ−′右）が最小となるようにの，椚を決定した．となる．ここでも微分を差分で表現し，右辺の（）内他のモデルについても同様である．

は1時点前の値で決まるとするとなお，r（乃）は最終需要研固定時のカルマンフィ

れ（乃＋1）−れ（乃）＝れ（乃）（α−∂れ（雅一1）−C了ち（乃−1））㈹

易化のため，中心化12ヶ月移動平均

右（乃＋1）一右（乃）＝右（乃）（d一〆れ（乃−1ト′右（乃−1））となる．ここで，んか1＝α−れ（乃−1）−C右（乃−1）

ルタのトレンド項を用いることもできるが，手順の簡錮T（症7）＝去如（乃＋ブ）旬（乃＋出））を用いた．
カルマンフィルタ近似によるモデルの適合度を評価

ん刀−1＝d一〆n（乃−1）−′箭（乃−1）
丁96（44）© 日本オペレーションズ・リサーチ学会. 無断複写・複製・転載を禁ず.

オペレーションズ・リサーチ

のパラノタが増加，すなわち，（久一♪2＝−3）と考えするのにると，（β∫＝0）となる斤む＝Eり／Euの値は表2のよう〃−1 且＝∑（y（乃＋1l乃）−〟（乃＋1））2になる．表1のEげを用いて計算された斤びの値を表〝＝l〃−23に示す．表3の値が表2の値よりも大きければAIC

g2＝∑（y（乃＋2l乃）一計（乃＋2））2
乃＝1の観点からモデル才の方がモデル1よりも優れている
12 与′＝∑（y（〃＋乃l〃）−y（Ⅳ十乃））2ことになる．

月＝1を用いる．ただし，

表1，3より，節3．1−3．3で提案したモデルは滑らy（乃＋椚l乃）＝〃 Ⅹ（乃＋∽l乃）＝〃 Fmx（乃l乃）であり，Fの要素には逐次，推定値を用いる．

かなトレンドモデルに優るとも劣らない予測力を持っているといえる．

4．2 評価結果（清酒の場合）（1）ビールの場合

清酒の場合には，表4から分かるように，Erl，Ef2

ロジスティックモデルでは，最小二乗誤差和且）を実現するパラメタ値はに関してはどのモデルも滑らかなトレンドモデルより劣っているが，最も重視したいEぴに関してはモデル

椚＝6700，仇＝0．00501，α2＝0．0044

3，4（BassモデルI，II）ともにモデル1（滑らかであったが，そのとき，Eo＝173384であった．最終なトレンドモデル）より改善されている．また，この需要∽として切りのよい∽＝10000を用いたときに例は過去のデータによく合うモデルが予測でも優れてはいるとは限らないことを示している．の＝0．00302，α2＝0．0029，且）＝173556

であり，E。の差は小さいのでこれらの値もモデル比

表2 仇＝0となる斤むの値

較時のパラメタ値として用いる．

節2のモデルをモデル1，節3．1，3．2，3．3のモデルをそれぞれモデル2，3，4と呼ぶ．モデル才の誤差

ノ

月む 1．0341

2

f1．03431．0869＆をgむとし，各モデルの誤差al，Ef2，Eぴを表1に示す．これから，モデル4（BassモデルⅠⅠ）のE2
表3 斤びの値を除き，モデル1（滑らかなトレンドモデル）よりも誤差が小さい．

2

ノ

モデル1とモデルブとのAICの差は

f尺2ノ（〝桓6700）

1．3Z951．3015

1．8Z94となる．ただし，凡は誤差旦′を求める際のデータ数

尺2ノ（m＝10000）

1．3Z74

1．31771．812Z

であり，♪fはモデルオのパラメタ数である．予測に尺3ノ（m＝6700）

1．32Z81．29291．0604尺3ノ（m＝10000）

1．3Z45

1．Z9501．3174尺。ノ（爪＝6700）

1．3ZZO

0．9776

1．76Z6

私ノ（∩咋＝10000）

1．3Z451．00891．7953仇＝凡log（且ノ／Eび）＋2（♪．−か）使うパラメタは研とのであるが，α2も含めて3個
表1各モデルの誤差（ビール）ら1節Z（f＝1）

ら2 ら

4918 4882 386
表4 各モデルの誤差（清酒）節3．16700 3699 3751 Zll（／＝Z）

10000 3705 3705 213

節3．Z

6700 3718 3776 364

節2（ノ＝1）（／＝3）

10000 3713 3770 Z93

節3．3
（／＝4）

2002年12月号

昂1

昂2

居げ

204

205

26．6

節3．1（J＝2） 211

221

26．8

6700 37ZO 4994 219

節3．2（J＝3） 206

216

23．9

10000 3713 4839 Z15

節3．3（J＝4） 218

255

23．1

© 日本オペレーションズ・リサーチ学会. 無断複写・複製・転載を禁ず.

（45）丁9丁表5 節3．4のモデルの場合独立に扱ったモデルで十分と思われる．今回は「ビール」と「清酒」のデータでモデルを比較したが，他の

gl

卑ビール

3676

313．3

清酒

206．1

20．79

データでも適合度や予測力を比較する必要がある．また，競合型Bassモデル［3，4］についても同様のカル
マンフィルタ近似が可能であるが，今回は比較対象にしなかった．同様のカルマンフィルタ近似は様々な非線形モデルに適用可能だと考えている．

なお／本研究は科学技術融合振興財団助成金による（ビールと清酒の競合考慮時）

ビールと清酒が競合しているものとして，節3．4の支援を受けたことを付記する．

方法を用い，ビールと清酒のトレンド成分をそれぞれ参考文献

れ（乃）∴㍍（乃）とした結果，表5が得られた．
ビール単独の場合と比較すると，節2の（滑らかなトレンド）モデルよりはよいが，節3．1の（ロジスティック）モデルとは与′に閲しで劣っている．清酒単独の場合と比較すると，節3．1の（ロジスティック）モデルよりはよいが，節2の（滑らかなトレンド）モデルとはElに関して劣っている．

［1］上田徹：「予測手法（1）：時系列予測法」オペレーショ
ンズ・リサーチ（1994．6）．

［2］上田徹：「予測手法（2）：生態学モデル」オペレーションズ・リサーチ（1994．7）．

［3］Tohru Ueda：

Demand Forecasting Models forMarketswithCompetition，ZTC−13，JensenandV．B．

Iversen（Editors），EIsevier Science Publishers，Pp．
26ト265（1991）．

5． まとめ［4］Tohru Ueda：
パラノタ固定時の二乗誤差の例E。＝173556と比べるとElの値ははるかに小さくカルマンフィルタによる精度向上が明白となった．また，節3で提案したモ

A Srudy of a Competitive Bass

ModelWhichTakesintoAccountCompetitionamong
Firms，Joumalqfthe（妙entionsReseaYChSocie＆〆ノ卸α乃，Vol．33，No．4，pp．319−334（1990）．
［5］阿部威郎，上田徹：「状態空間モデルを用いた電話収入

デルは既存の滑らかなトレンドモデルに優るとも劣らない結果が得られたのではないかと思う．特に，E2，旦′に関してトレンドモデルよりも誤差が小さく，長期の予測への対応力がトレンドモデルよりも大きかったといえる．しかし，最後に提案−した競合を考慮にいれたモデルは，モデルの複雑性が増す割には予測力，予測」，電子通信学会論文誌，Vol・J68−A，No・5，Pp・437二443（1985）．

［6］R．ハーバーマン著，稲垣宣生訳：「個体群成長の数学モデル」，pp．130−138，現代数学社（1981）．
［7］「家計調査年報」，総務庁統計局．

適合度が圧倒的に優れているわけではなく，それぞれ798（46）© 日本オペレーションズ・リサーチ学会. 無断複写・複製・転載を禁ず.

オペレーションズ・リサーチ