北九州工業高専添田満 北九州工業高専古谷忠義

九州工業大学研究報告（工学）No．63 1991年9月       g

欠落のあるデータを用いた自己回帰モデルの推定

      （平成3年5月28日 原稿受付）

北九州工業高専添田満

北九州工業高専古谷忠義 設計生産工学科松岡清利

On Parameter Estimation of Auto−regressive Models       without Complete Data

by Mitsuru SOEDA

  Tadayoshi FURUYA   Kiyotoshi MATSUOKA

Abstract

       It is often required to estimate parameters of nonphysiとal models， where some observations arρ      absent． Data absence at some sampling instants happens in such situations as enviro㎜ental data

．       In this paper， a method is proposed for parameter estimation of auto−regressive processes with．

     out some points of data． The sampling points， whose data are lost， are assumed to be random with      aprobability as a Bemoulli sequence． This method consists of two parts；one is the part of pre−

     treatment of missing data， the other is the part of parameter estimation with pertreated data． In      the pretreatment， missed data are extrapolated with simple statistical values by the four different      ways． And the parameters of AR processes are estimated by the conventional least squares      method．

       It is shown that this method is practical and effective by the examples of the prediction of con．

     centration levels of pollutants of air pollution．

       まま用いたのではうまくモデリング・予測を行うことが  1．はじめに

       出来ない。

 データ処理において，観測データに非物理モデルの一   このような観測データの欠落は周期的に発生するもの つである自己回帰モデルを用いてモデリングを行い，そ  とランダムに発生するものが考えられる。周期的に欠落 の状態の予測・制御を行う方法がある；1）（2）これらの方  のあるデータに対してはスペクトル密度や分散の推定に 法は，基本的には観測データがある一定の時間間隔毎に   関する方法が報告されている④（5）。確率的に欠落があ 得られる条件の下に行われる（3）。      る場合，その発生がベルヌーイ過程に従うものとしてス  しかし実際の場合，データが明らかに異常で除外せざ  ペクトル密度の推定法が報告されている（6）。しかしこ

るを得なかったり，大きな雑音が発生したり，観測器の  れらは周波数領域におけるスペクトル解析が中心で環境

故障やその他の原因で観測時点における観測データが欠  制御など予測法を中心とするような処理には向いていな

落していることがよく起こる。こうした場合，一定時間  い。ランダムに欠落があるデータに対する時間領域での

間隔毎に観測データが得られないため従来の手法をその  処理としては，欠落をサンプリング間隔の変化としてと

らえる同法定（7），欠落データとシステムパラメータを

剛推定する方法…ぷ分をあてはめフ，ルタを繊、システム

^{@ 出力時系列｛σ ｝（ ＝1・2・ ）}

するなどの統計的方法（9）（10）などが提案されている。し かし実際の場合には欠落の有無により，これらの推定法 と欠落の無い時の推定法を切り換えて併用しなければな

らず，アルゴリズムが複雑になるなどしてオンラインの         観測器の故障 時系列予測が中心となる問題にはむいていない。一方，，        明かなデータ異常 欠落の無いデータに対しては多くの優れた推定法が提案

されている。

 そこで本論文では・確率的に欠落を生じる時系列デー      観測値の欠落 タに対して，その欠落部分を簡単な量で補間する前処理

を行い，前処理後のデータに欠落の無いときの推定法を

用いて自己回帰モデルをあてはめる方法について報告す      欠落観測値 ｛XJ（ゴ＝1，2，…）

る・補間する量としては簡単な4つの統計量を与えてい        。、データ嫡

る。そしてこの方法を大気汚染レベルの実データに摘要 し，モデリングの精度を汚染レベル値の予測精度の面か ら比較検討した。

2．観測値の欠落と前処理       図一1 観測値の欠落  一定の時間間隔でシステムから観測されるデータを接

続してできた時系列を｛σ9｝（ ＝1，2，…∧りとする。

 これらの観測値がランダムに欠落を生じ，データが入   平均値を補間し 手できるか否かはベルヌーイ過程に従うものとすると図       1M

、のように実際の観測甑｝（ 一、，，，…めはつぎ 炉7滅濁

のように表現される。       とする。ここで凡は入手されたデータ個数である。

  X ＝4ゴ・ひ      方法3

ここでは｛4f｝は｛ひ｝とは独立なベルヌーイ過程で0と1   欠落を生じている時点 には，移動平均値を補間し の値をとり，その確率は

  P。｛己一。｝一、−P         尤一詰塩（・〉・）

で与㌻鳳｝：：−P）が欠箒を裁   x＝÷紅（・≧ ≧1）

 このような欠落を生じたデータに対してモデリングを   とする。ここで5は移動項数である。

行うために，前処理として欠落箇所に以下に示す4種類   方法4

の量で補間を行う。      たとえば環境データには1日24時間毎にある程度周期        性があるであろうという観点から，、このような観測デー 方法1       タについては，あらかじめ入手されている過去のデータ  観測データの欠落を生じている時点τには，その直前   を周期毎に分割し，1周期における各時刻（1時から24

に得られた観測値で補間し      時）ごとに平均値を求めておく。そして欠落を生じた時

  X，＝X∫．、      点τにはそれに対応する1周期における時刻の平均値を

とする。      補間する。すなわち欠落時点τがK周期めの時刻∫に対

方法2      応するとき，時刻∫の全周期にわたるデータ平均値Mj

 欠落を生じている時点rには，入手出来た全データの   を求め

欠落のあるデータを用いた自己回帰モデルの推定       11

  X｝＝M／       1のベルヌーイ系列を順次発生させ，0の時それに対応 とする。       する時点のデータは欠落しているものとして欠落データ        を作成した。

 以上の方法により補間されてできた新たな時系列デー   欠落率（1−P）を変化させ種々の欠落を有したデー タ｛X、｝（i＝1，2，…めを用いてモデリングを行う。   タを作成し，これらのデータに前述の各前処理を施し，

      自己回帰モデルをあてはめ，得られたモデルから汚染レ

3・モデリング       ベルの予測を行った．モデリングに使用したデータは予

 このようにして得られたデータは一般に非定常性を有   測する日の前日までの1カ月間，744個のデータである。

すると考えられる。そこでデータの非定常成分を除去す   データに欠落が無いときの従来のARモデルあては るため，時系列の1次差分を計算することにより，新た   めによるレベル値の予測結果を図2に示す。図3は欠落 にデータ系列｛K｝を作る。      率30％の欠落があるデータに対して前処理を施さず，そ   酩＝X−X，．、       のまま従来のARモデルあてはめにより予測を行った

この時系列の平均値      結果である。このようにデータの欠落が多くなってくる r一 ﾒ     と曝蒜よ㌶㌫三芸二霊鑑

を求め      前処理を行った後にARモデルをあてはめ予測を行っ

．  酩＝yi一γ   ＝1，2，…，」V）

 より，あらためて時系列｛酩｝を作る。

  この時系列に自己回帰モデル

   酩＿￡。，．恥＋、8       −t・u・value

     ノ＝1      ・・・…  predicted value

    ここでθ、は正規性白色系列

をあてはめる。       δ いま｛幻の共分散

  ＆−N』κ署酩・呂袖

を求めると，次数mにおけるARパラメータ推定値∂，

はレビンソンのアルゴリズムにより計算される。また最 適次数初＝吻。は，AIC基準値を求め，これが最小にな

σう

るように決定する！11）       Time［h］

 こうして得られたモデルを用いるとYiの予測値は        図一2 データに欠落のない時の予測結果

  へ          γ，＝Σ〜7」γi一ゴ

により計算される。

 従ってシステムの出力時系列の予測値は         N

  ＿      ＿  ＿      O

  x《＝x，＿1十γ，十yP       uり により求められる。

4．数値実験例

 昭和49年5月および6月に徳島県那賀川町で1時間ご とに毎日観測されたSO2の大気汚染レベルデータを用 いて実験を行った。これらの実データは欠落の無いデー

タであり，実験ではパソコン内に一様乱数を用いて0と     図一3 欠落を有する時の予測結果（前処理なし）

た結果を図4，図5に示す。前処理を施すことにより欠    次に，各方法によるモデリング精度の比較検討を行う 落の無い場合に近いモデリング・予測が行われているこ   ため，30日間にわたって予測を逐次行い，その予測誤差

とがわかる。       の評価を次式により行った。

       true value

δ       δ oり    1       ．       の

     〜         it   n      l         ∫      i，      ・     ξ              8     言      ：      ，      」，

      ．λγ     r  、．◆ 、 ∫ ．・

       ．  ・  、     、 ξ   ・  c  ・・

δ

24     48     72     96     120       24     48     72     96     120

   方法1による予測       方法2による予測

      パ

       i    i  了        ｝・       ∫

．グ  ． γ・ ．2  『

       了  、  、3   ぎ

24     48     72      96     120      24     48     72     96     120

   方法3による予測       方法4による予測         図一4 欠落率10％における予想結果

      一 true value       ・・・…  predicted value

δ      δ

の       i      ，       め

 ． ．タ ．・一・．・ ｛ ・

δ

 、，・タ』γ・・、ハ日・．

駕  ．              、

24     48     72     96     120       24     48     72     96     120

   方法1による予測       方法2による予測       δ

      の

．．、γ・陣h．パ烈

．，、タぽ5．！ ？・、

24     48     72     96     120       24     48     72     96     120

   方法3による予測      方法4による予測

        図一5 欠落率30％における予想結果

欠落のあるデータを用いた自己回帰モデルの推定       13

    ムド       

    Σ（X −X・）2      5．おわりに   力＝同N

      評2         観測値に確率的な欠落を生じるデータに対して，欠落

また環境のデータ処理などにおいては，予測を行って警   部分を簡単な量で補間することによりモデリングを行う 報を発する場合など，予測値が実測値より大きく予測さ   方法を提案した。補間量としては4つの統計量を考え，

れる傾向を持つ方が望ましい。そこでこのことを含めた  実際に欠落のある環境データに適用しモデリング・予測 評価を行うため次式の評価値を計算した。測定値より大   を行った。そしてそれらの有効性を予測精度の面から比

きめに予測された場合，このJ2は小さくなる性質を持  較検討した。その結果，データに欠落がある場合，従来 っ。       のような複雑な処理をしなくても，これらの簡単な前処      N  一      理を行うことでモデリング・予測が可能なことがわかっ     Σ（Xi−X ）2

  五＝宗  ＿      た・

    Σ（X〜十Xゴ2）

    i＝1

 種々の欠落率に対する各方法の1時間先の予測精度      参 考 文 献

Jl，J2を表1，表2に示す。      （1）Y・Sawaragi， et aL：The Predictions of Air Pollu一  いずれの方法も欠落率が5％程度まではほぼ同程度の     tion Levels by Nonphysical Models Based on Kal一

精度でモデリングされるが，さらに欠落率がふえてくる     man Filtering Method， Trans． ASME， Joumal of と方法1が最も精度がよい。また他の方法と比べて予測     Dynamic Systems， Measurement and Control，

値も大きめに求められ，環境データの処理などに有効で     series G， vol．98，375／386，1976

あることがわかる。方法2と方法4については欠落率が   （2）大恵：空間時系列モデルの大気汚染レベル値予測へ 30％程度まではその効果はあまり変わらないが，さらに    の適用，システム制御情報学会論文誌，vol．4，

欠落率がふえてくると方法4の方が精度がよいことがわ・    No．2，61／69，1991

かる。      （3）G．E． P． Box and G． M． Jenkins：Time Series       Analysis， Holden．Day，1970

      （4）R．H． Jones：Spectral Analysis with Regularly Mis−

      sed Observations， Ann． math． Statist．，33，455／461，

      1962

      表1 欠落率と予測誤差評価J1        （5）EParzen：On Spect．al Analysis with Missing       Observations and Amplitude Modulation， ASNK−

      HYA， ser． A，25，383／392，1965

      （6）P．A． Scheinnok：Spectral Analysis with randomly       missed observations．the binomial case， A皿． Math       Statist．，36，971／977，1965

       （酬  （7）RBM・giffi・・ndD・N・PM…h・・P・…・…

      estimation for auto−regressive systems with mis−

      sing observations， Int． J． SYSTEMS SCI．， Vol．11，

      表2 欠落率と予測誤差評価J2

      No．9，1021／1034，1980

      （8）稲葉，塩谷：不完全データによるARモデルの同       定問題，第26回システムと制御発表講演会，69／70，

      1982

      （9）酒井，徳丸：欠落データのある自己回帰過程のパラ       メーター推定，第10回確率システムシンポジウム，

       （E．。2）   77／80・1983

欠  落  率（％）

方法

0

20 30 40 50

1

4．46 4．46 4．48

4．51

5．06 5．41 5．76

4．47

4．63

4．57 4．83 5．07

5．60

欠 落 率（％）

方法

0

20 30 40 50

1 2．31 2．32 2．32 2．32 2．32

2．33

_{2．31 2．34 2．39 2．54 2．69}

2．31 2．32 2．34 2．38 2．41 2．47 2．53

2．31 2．35 2．40 2．55 2．70

3．01

⑩ 立花：情報理論とその応用研究会第2回研究討論会   資料，35／42，1979

（1D 木村：ディジタル信号処理と制御，昭晃堂，1982

（吻 添田，片平，葉賀：日本機械学会講演概要集，Nα