技術レポート: 時系列解析事例集

技術レポート: 時系列解析事例集
知識工学部 長沼 茂太
平成 21 年 11 月 11 日
1
概要
昨今の情報科学技術の発展により大量のデータを扱うことが容易になってきたことに伴い、そこ
から有意な情報を抽出する解析技術の需要は高い。本レポートでは、そのような解析技術の中でも
時系列解析に着目し、現在開発中の時系列解析エンジンを用いて実際の時系列データに対する解析
例を示すことで、時系列解析およびエンジンの有用性を示すことを目的としたものである。
なぜ時系列解析か
時系列データとは、ある期間にわたって一定の観測対象の値を時刻情報と共に
採取しつづけたものである。本稿でもいくつか例を示すが、世に蓄積されているデータ・問題とな
るデータで、時系列データと呼べるものは非常に多い。およそ ’ ログ’ と名の付くデータには時刻
が記載されているであろう (web アクセスログ、出退勤記録 etc ...)。
通常の解析では物理現象を定量化した観測値のみが対象で、個々のデータ間に順序は考えない
が、時系列解析ではそこに時間という観念を定量化した時刻情報を付与し過去・現在という流れの
中で観測値を捉え記述することにより現象理解およびその予測により有益な情報を得ることができ
る。このことは次のような簡単な例からでもわかる。以下の表は、ある地点の 4 時間ごとの気温お
よび湿度を 2 日間に亘って示したものである。
時刻
1
5
9
13
17
21
1
5
9
13
17
21
気温
17.8
77
15.8
79
22.4
51
25.3
44
23.1
60
21.4
70
17.7
86
16.8
85
18.4
72
21.4
61
20.9
62
19.2
62
湿度
気温と湿度が、一方が変化する方向と逆の方向に変化するという負の相関をもっていることは、
その観測値からだけでも明らかである。しかし、次の 1 日でこれらはどう変化するのか、という問
いに対してはそれだけではわからない。そこで時刻を考慮することにより、これらの値に昼夜の周
期があるということがわかり、予測も出来得るだけの現象に対する知識を獲得できる。本稿では、
より詳しく時系列解析の要諦を示すために、汎用的な解析例やデータ固有の特徴・問題に則した応
用例を示す。
なぜエンジンを開発したか
このことについては 6 章で詳述するが、時系列解析をできる機能を
もった統計解析エンジンは R, VMS, S-plus と様々あるが、操作が煩雑で、少し突っ込んだ解析を
しようと思えば、ユーザーが自らスクリプトを書かなければならない。よって、インタラクティブ
な環境で簡便に、しかも柔軟に機能を追加できるエンジンが在れば、それは意義あるものと考え、
エンジン開発に至った。
1
2 章では例に用いる時系列データについて、3 章から 4 章まででそれらの時系列データを用いた
実際の解析例を、6 章では時系列解析エンジンについて、7 章では総括・今後の展望を述べ、まと
めとする。
本レポートによって、時系列解析とエンジンに興味を抱いてもらえれば幸いである。
目次
1
概要
1
2
時系列データ
3
3
UKgas
3.1 記述 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5
5
3.1.1 acf . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.1.2 トレンドの推定 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
モデリング・予測 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5
5
7
AR モデル . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
HoltWinters モデル . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7
8
RTT
4.1 移動平均 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
11
11
4.2
4.3
13
15
3.2
3.2.1
3.2.2
4
5
OBS
5.1 追従しなくなる故障
5.2
6
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
17
17
追従の仕方が異常である故障 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
19
時系列解析エンジンについて
6.1
6.2
7
ccf . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
periodgram . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
24
開発概要 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
意義 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
まとめ
24
24
27
2
2
時系列データ
本稿で解析対象として用いるデータは UKgas, RTT, OBS の三種類である。これらを図示する
ことによって以下に示す。図示することは時系列解析の第一歩として、そのおおまかな特徴を捉え
る上で重要であり、今後どのような解析を行なうかの指針を立てることもできる。
• UKgas
UKgas(図 1) とは 1960 年から 1986 年までの各季節ごとのイギリスにおけるガスの使用量
のデータである。その値には、大局的には徐々に増加していく傾向と、毎年季節ごとに同じ
ような変動を繰り返しており、その変動幅が増大していくのが見て取れる。
図 1: UKgas
• RTT
RTT(図 2) とはネットワーク上の各ルータに対して、一定時間間隔 ( 5 分) で ping を打ち、
反応が返ってくるまでの時間を採取したデータである。図はある一つのルータのデータであ
る。一見したところ、取り得る値に下限値がありそうだが、変動が激しく、増減に明らかな
周期性は見られない。
図 2: RTT
3
• OBS
OBS とは自動車の各部品に取り付けられたセンサー値 (パラメータ) を、故障発生時に 15 秒
間 ( 0.2 秒毎 ) にわたって採取したデータである。さまざまな故障について、OBS は約 14000
件収集されているが、本稿で着目するのは、パフォーマンス系の故障である。パフォーマン
ス系の故障とは、正常時においてはある二つのパラメータ値が、一方のパラメータに反応し
て追随するものが、故障時に全く追随しなくなったり、反応が鈍くなったりするもののこと
である。図 3, 4 は OBS データにおけるパラメータ AF SENSOR と HO2S S2 の正常時と
異常時の例である。正常時は HO2S S2 の値は AF SENSOR の値に同程度の変化幅で時間
差なく追随するが、異常時にはその反応の度合いが鈍くなっていることがわかる。
図 3: OBS:正常時の HO2S S2 と AF SENSOR
図 4: OBS:異常時の HO2S S2 と AF SENSOR
4
3
UKgas
まずは典型的な時系列解析の流れを UKgas(図 1) を用いて説明する。概略は以下である。
1. 記述
時系列を図示したり、自己共分散・自己相関・ピリオドグラムなどの記述統計量を用いて時
系列の特徴を簡潔に表現する。
2. モデリング
与えられた時系列に対し、その変動を表現する時系列モデルを構成する。記述によってその
時系列の特徴は様々に表現されるので、対象・目的に応じた適切な時系列モデルを選択し、
モデルに含まれるパラメータを推定する。
3. 予測
モデリングによって得られたモデルと、現在までに得られた情報から今後の変動を予測する。
3.1
記述
2 章で、時系列データ UKgas の値には季節成分・トレンド成分があるように見られると述べた。
そのような観察は、時系列解析においては次に示すような方法で記述される。
3.1.1
acf
時系列データの特性を表す基本的な統計量として、平均・自己共分散・自己相関などが良く知ら
れているが、ここでは自己相関関数 (以下 acf ) を用いて、UKgas の周期性を表現する。
データ間の関連を表す一般的な統計量として相関係数があるが、acf は相関係数に時間差（ラグ）
の概念を導入したものである。すなわち元のデータと、時間を一定幅ずらしたデータとの相関係
数を求めることになる。もしデータが周期 t で同じような変動を繰り返すならば、ラグ t の整数
倍ごとに相関係数は高くなるはずである。例えば sin 関数のような 2π ごとに同じ波形を繰り返す
データであれば、2π ずらす (ラグ 2π) ごとに元のデータと全く重なるので、自己相関係数はラグ
2π ごとに高くなる。
実際に UKgas に対して、自己相関関数を求めた結果が図 5 である。横軸がラグを表し、縦軸が
相関係数値である。図 5 を見てみると、ラグ 4 ごとに相関係数値が高い値を示している。UKgas
の時間単位は季節なので、4 季節すなわち 1 年ごとの周期があると考えられ、観察を裏付ける結果
となっている。
3.1.2
トレンドの推定
時系列データには長期間にわたって増加や減少が続く傾向が見られるデータがある。UKgas(図
1) も全体的に増加傾向にある。このような長期的傾向変動を トレンド と呼ぶ。時系列データの趨
勢の把握にはトレンドを推定し、記述することが有効である。そのために時系列モデルとしてトレ
ンドモデルを利用した方法があり、時系列データのトレンドを簡潔かつ具体的に表現することがで
きる。
5
図 5: UKgas の自己相関関数
図 6, 7 は UKgas にトレンドモデルを適用しトレンドを記述した結果である。青い線が UKgas
で、緑の線がトレンドである。トレンドモデルにはいくつかパラメータがあり、それを調整するこ
とによって、図 7 のように短期間の動きに反応するようなトレンドを表現したりすることもでき
る。どのようなトレンドが妥当かは観察者が取捨選択することになる。UKgas に関しては、季節
変動が激しく、トレンド自体はゆるやかに単調に増加しているので、図 6 の方がより妥当と考えら
れる。
図 6: UKgas のトレンド 1
図 7: UKgas のトレンド 2
6
3.2
モデリング・予測
時系列モデルを用いることにより、時系列の特徴を簡潔に表現し、シミュレート・予測すること
ができる。しかしモデルおよびモデルにおけるパラメータは解析対象データに則したものを選ぶ
べきであり、その判断のためにいくつか指標がある。以下では、時系列解析における基本的な AR
モデルと、より UKgas に則した時系列モデル HoltWinters モデルによるモデリングを示し、比較
検討する。
3.2.1
AR モデル
AR モデルとは、別名自己回帰 (AutoRegression) モデルと呼ばれ、時系列の各値を過去の値に
よる 1 次多項式で表すものである。すなわち係数パラメータを {a1 , . . . , ak }, 時系列を yn とする
と、次の式で定義される。
yn =
k
∑
ai yn−i
i=1
次数 k (どこまで回帰するか) の値に対し係数パラメータ {ai } を推定するアルゴリズムとして、
yule-walker 法, burg 法があり、推定されたモデルを評価する指標として AIC がある。AIC とは
簡単に言えば ’ モデルの尤もらしく無さ’ と ’ モデルのムダ’ を足し合わせたもので、その値が低
いほどモデルと実際の時系列は適合度が高いと考えられる。表 1, 2 は、それぞれ次数 20 までの
yule-walker 法, burg 法で求めたモデルの AIC である。
Order
AIC
Order
AIC
Order
AIC
Order
AIC
0
1501.27
11
1097.15
0
1501.27
11
1298.38
1
2
3
1460.96
1456.80
1173.92
12
13
14
1095.52
1097.52
1099.38
1
2
3
1462.69
1459.58
1313.48
12
13
14
1299.80
1301.70
1303.42
4
5
1118.77
1118.38
15
16
1101.31
1099.15
4
5
1295.85
1296.60
15
16
1305.39
1307.29
6
7
8
1115.61
1112.01
1096.43
17
18
19
1101.00
1098.82
1099.72
6
7
8
1290.36
1292.24
1294.23
17
18
19
1309.19
1311.13
1313.11
9
10
1095.01
1096.90
20
1099.88
9
10
1294.51
1296.41
20
1313.91
表 1: burg 法による AR モデル
表 2: yule-walker 法による AR モデル
次数 0 から AIC は減少していき、yule-walker 法は次数 6, burg 法では次数 9 で最小を示し、
その後は増加している。さらに AIC の値をアルゴリズム間で比較すると、全般的に burg 法の方
が値は小さく、ukgas を推定する上では yule-walker よりも適したアルゴリズムといえる。
AR モデルによる予測 構成した AR モデルによって、UKgas が今後どのような値をとっていく
か予測することもできる。図 8 は次数 9, burg 法 による AR モデルで、30 年先までの予測を行
なった結果である。青線が UKgas で、緑線が AR モデルによる値、縦の青い点線以降が予測値で
7
ある。予測以前までは良く UKgas の値に モデルの値が沿っており、平均絶対誤差も 6.45 と低い
値であるが、では予測値はとなると、トレンドは減少傾向にあり、周期的な変動幅も減少し、最終
的に値は UKgas の平均値 337.63 に収束してしまっている。実際のデータを観察する限りではこ
の逆すなわちトレンド・周期的変動共に増加傾向を示しており、一つの値に収束することは到底考
えられず、この予測値が妥当であるとは言えない。
実は AR モデルは定常的1 な時系列に対して有効であり、UKgas のようにトレンドがあるよう
な非定常なデータの予測には適さない。
図 8: UKgas と AR モデル
よって次にトレンドや周期性をもった時系列データに対して有効な HoltWinters モデルについ
て述べる。
3.2.2
HoltWinters モデル
HoltWinters モデルは別名指数平滑化法と言い、過去のデータが現在のデータにおよぼす影響
が、時間的に遠いほど指数関数的に減少するので、そのように呼ばれる。
HoltWinters モデルでは、データ値は三つの成分、基準値・トレンド値・周期値から成り立って
いると考え、各成分にパラメータ α, β, γ を設定することによって、過去のデータからの各成分へ
の影響度合いを決定している。
各パラメータは 0 から 1 までの値をとり、1 に近いほど直前の観測値を、0 に近いほど長期の変
動を重視してモデルの値が算出される。図 9, 10, 11 にパラメータ α によるモデルの差を示した。
青線が観測値で、緑線がモデル値であるが、α が 1 に近づくにつれ直前値を重視した値が、0 に
近づくにつれ長期の変動を重視した値がモデル値となっていることがわかる。問題点としては、1
に近すぎると特異値が予測値に大きな影響を与えてしまい、0 に近づけすぎると観測値の変動に対
して予測値の追随が遅れてしまうということがある。
α, β, γ の決定 適した 3 つのパラメータ α, β, γ を決定するために時系列解析エンジンでは、誤
差評価関数によって観測値とモデル値の誤差を算出し、自動的に最も小さい誤差関数値を出力する
1 非定常データとは時間によって、統計的特徴が変化するデータのこと。よって定常的とは、時間と共に平均または平均
の周りの変動が変化しないことを言う
8
図 9: α = 0.1
図 10: α = 0.5
図 11: α = 0.9
組み合わせを選んでいる。
観測値 {yi , . . . , yn }, モデル値 {m1 , . . . , mn } とすると、表 3 は、誤差評価関数として
n
∑
M SE =
(yi − mi )
i=1
M AP E =
n
∑
n
|(yi − mi )/{(yi + mi )/2}|
i=1
n
∑
RAE = ∑
n
|yi − mi |
i=1
|yi − mean(y)|
i=1
を用いてパラメータを選定した結果である。seasonal というのは、データ値を構成するときに季
節成分をどう寄与させるかについてで、加法的 (additive) または乗法的 (multiplicative) がある。
例えば、ある商品の販売量が、毎季節ごとに一定量ずつ増えていく場合には additive、一定割合ず
つ増えていく場合には multiplicative のモデルを選ぶのがよい。
誤差関数
seasonal
誤差関数値
α
β
γ
MSE
MSE
MAPE
additive
multiplicative
additive
1276.80
1087.69
0.09
0.04
0.02
0.11
0.97
1.00
0.86
0.99
0.80
1.00
MAPE
RAE
multiplicative
additive
0.08
0.13
0.02
0.09
0.89
1.00
0.92
1.00
RAE
multiplicative
0.11
0.05
0.31
0.95
表 3: いろいろな誤差評価関数による比較
表 3 を見ると、ほとんどの場合で α は 0 から 0.1、β および γ は 0.8 から 1 となっており、
誤差関数および seasonal によって違いはあまり見られない。誤差関数値は UKgas の取り得る値
9
が 84.8 から 1163.9 であることとそれぞれの定義から判断すると十分小さい値となっており、そ
れぞれ観測値に適合している。seasonal に関しては multiplicative がどの場合でも誤差関数値が
additive より小さい値となっており、multiplicative が適していると判断できる。図 12 は誤差評
価関数 を MSE, seasonal を multiplicative として 10 年先まで予測をしたものである。AR モデ
ルとの違いが予測値に顕著である。トレンドおよび周期成分の変動が増加傾向であることがよく表
現されており、AR モデルと比較しても大幅にこちらの予測の方が妥当である。
図 12: UKgas と HoltWinters モデル
さらに誤差を小さくする方法として、パラメータの精度を細かくする方法がある。表 3 では精
度を 0.01 すなわち小数点以下第二位までとして算出していたが、精度を 0.002 までにして、誤差
評価関数 を MSE, seasonal を multiplicative として求めた結果、α = 0.022, β = 1, γ = 0.798 と
なり、誤差関数値は 1087.01 であった。精度を上げると計算すべきパラメータの組み合わせも増え
てしまうため、0.01 だと数分だが、0.002 とすると数時間かかってしまう2 。改善される誤差関数
値 ( 0.5 程度 ) からすれば、UKgas の場合 0.01 で十分であるが、パラメータの数値的最適化の手
法は時系列解析エンジンにいずれ加えなければならないだろう。
2 計算量は、三つの各パラメータで
5 倍となるので合計で 125 倍である
10
4
RTT
ネットワーク上のあるルータに対して ping を打ちそれが戻ってくるまでの時間を RTT と呼ぶ。
2500 個のルータに対し、300 秒間隔で 7746 回 RTT を採取した。それらのルータに時系列解析を
応用することで、なんらかの特徴を見出すことが本章の目的である。主に解析対象とするのはその
1, 2, 207 番目のルータである。それぞれ RTT-1, RTT-2, RTT-207 と呼ぶことにする。
図 2 は実は RTT-1 のデータである。先にも述べたが、波形を見ても変動が激しいため、下限値
の存在程度のことしか洞察できない。一般的に非定常なデータをある程度定常的にするための処理
として、対数変換・差分・移動平均などがある。RTT-1 のような変動の激しいデータは、移動平
均をとり滑らかにすることで、全体的な増減の傾向が明確になることがある。
4.1
移動平均
移動平均とは、ある時間におけるデータの値をその前後一定範囲の値の平均として決めるもので
ある。観測値 {y1 , . . . , yn }, 一定範囲幅 k, 移動平均値 {m1 , . . . , mn } とすると
k
∑
mi =
yi+j
j=−k
2k
, i = 1, . . . , n
である。移動平均をとれば、もともとの値が周りと比して特異な値であったとしても、移動平均
をとることでその近辺での平均値に値は引っ張られる。よって全体の結果としては、あらかじめ決
めた k に応じた長さの期間での動きが反映されている。
今は、RTT-1 について何もわからないので、短・中・長期間さまざまな k でその動きを見てみ
たい。図 13 は RTT-1 に対して、k = 5, 20, 100, 200, 500, 1000 としてそれぞれの移動平均をグ
ラフ化したものである。
k = 5, 20, 100 のグラフでは、k が増えるに従って、時間にして約 300 の間隔での増減の繰り返
しが滑らかに見られるようになり、その後 k = 200, 500, 1000 では、平坦になっていき、1000 で
はほどんど一定値で 13 から 13.5 の値をとっている。RTT-1 の平均が 13.63 であることからこれ
以上長い期間での動きは RTT-1 にはないと考えられ、さらに、約 300 の周期が RTT-1 の大勢を
占めていると推測できる。
3.1.1 節で述べたように、時系列データに周期性があれば、それは acf で裏付けることができる。
周期は 300 と予想されるので、三周期分ラグ 900 ぐらいまでの自己相関を求めた結果が図 14 であ
る。顕著にラグ 288 ごとに相関係数値が高い値を示しており、周期 288、 実際の時間に換算して
ちょうど 1 日であることがわかる。
11
図 13: RTT-1 におけるいろいろな移動平均
図 14: RTT-1 の acf
12
4.2
ccf
いま RTT-1 の周期は判明したので、ほかのルータとの関係も見てみたい。時系列データにおけ
る各変量間の関係を見る手法として、相互相関関数 ccf がある。ccf は acf と似ており、ラグごと
の相関係数を求めるものであるが、ccf は二つの変量間に定まる値で、一つの変量ともう片方の変
量の時間をずらしたものとの相関係数である。
図 15 は RTT-1 と RTT-2 についてラグ 600 まで ccf を求めたものである。一般的に acf は自
身との相関係数なのでラグの正負で対称となる。よっていままでラグは正の部分しか見てこなかっ
たが、ccf はそうはならないので図 15 ではラグ -600 から 600 までを求めている。
図 15: RTT-1 と RTT-2 の ccf
しかしながら図 15 はほとんど左右対称で、ラグ 288 ごとに相関係数値が高い値を示しており、
RTT-1 と同様の周期が推測される。実際に図 16 で RTT-2 の acf を見てみると、RTT-1 と非常
に似た形で、ラグ 288 ごとに相関係数値が高い値を示しており、周期 288 であることがわかる。
図 16: RTT-2 の acf
他の番号のルータに関しても同様の結果が得られたが、いくつか相関が崩れているものもあっ
13
た。図 17 は RTT-1 と RTT-207 の ccf であるが、図 15 とは明らかに違い、規則性はなく相関係
数の値も小さな値である。RTT-207 の acf ( 図 18 ) を見てみると相関係数値は ラグ > 0 で 0.04
以下と非常に小さい値をとっている。
図 17: RTT-1 と RTT-207 の ccf
図 18: RTT-207 の acf
RTT-207 のような acf のグラフを見て思い当たるのが、白色雑音 (white noise) である。白色雑
音とは簡単に言えば規則性がないノイズのことであり、白色雑音であれば、acf のグラフは ラグ 0
で 1, ラグ > 0 で 0 となる。RTT-207 は白色雑音に近いと考えられるが、次節ではそのことにつ
いてもう少し詳しく見てみる。
14
4.3
periodgram
観測値は色々な周波数を持つ波 (以下 周期成分) を合成することで表出する値である、と見たと
きにどの周期成分がどれだけ強いか、を解析する方法をスペクトル分析と言う。例えば観測値に顕
著に見られる周期性があったりとすると、スペクトル分析をすれば、その周波数をもつ周期成分の
強度が著しく強くなる。特に規則性がない白色雑音は、スペクトル分析をすると、どの周期成分の
強度も同じとなる3 。
時系列解析エンジンではピリオドグラムを求めることによって周期成分の強度を一定間隔の周波
数について推定することができる。データ量が多いほど、求められる周波数の精度も細かくなる
が、ピリオドグラムはスペクトルの一致推定量4 ではないため、ここではデータを分割して平均を
とることによって真の値に収束させるようにしている。
図 19 は RTT-1 のピリオドグラムである。他の周波数と比べて特に 7/2048, 7/1024, 21/2048
の周期成分の強度が強く、これらはそれぞれ周期 約 292, 146, 97 である。縦軸は周期成分強度の
対数値をとっているので、一番強度が強い 10.42 の値をとっている周期 292 の周期成分が RTT-1
では支配的と考えてよいだろう。これは移動平均や acf による解析から判明した周期 288 という
結果に非常に近い値であり、1 日周期であるということを強く裏付ける。
図 19: RTT-1 のピリオドグラム
一方、図 20 は、RTT-207 のピリオドグラムである。RTT-1 のピリオドグラムのときのように
ピークは見られず、強度は最大でも 6.7 であり RTT-1 に比して小さい。最小・最大幅も 4.53 と
これも RTT-1 の 7.9 という値よりも狭く、白色雑音の特徴5 により近くなっている。
表 4 は、 1000 番までのルータについてピリオドグラムを算出し、最大強度の周期成分の周期と
強度について集計したものである。RTT-1 のように周期 292, 強度 10 以上というのが 901 個あ
り、大半のルータは RTT-1 のような 1 日周期であると考えられる。
表 5 は、周期 292 の周期成分が最大強度をとらなかったルータ (RTT-207 含む) に対して、強
度の最小最大幅を集計したものである。ほとんどが RTT-207 と同様に狭い値を示している。これ
らのルータの RTT が白色雑音だったとすれば、ping の遅れに規則性がない、すなわちそのルー
タへの負荷には規則性がないことになる。
3 このような特徴をもった波形の光は白く見えることから、白色雑音と呼ばれる。
4 一致推定量とは、ある推定量で、標本が増えるほど真の値に収束するもののことである。
5 周波数によって強度が変化しない
15
図 20: RTT-207 のピリオドグラム
強度
周期
∼ 15
∼ 10
∼9
∼7
∼6
∼5
2048
1024
682
0
0
0
0
0
0
0
0
0
3
1
0
1
0
2
3
0
0
7
1
2
512
292
256
0
1
0
0
900
0
0
17
0
0
1
0
0
4
2
1
1
0
1
924
2
85 ∼
0
1
0
2
27
33
63
1
901
17
7
36
38
1000
表 4: 1000 番までのルータの最大強度周期成分集計表
強度の幅
∼7
∼6
∼5
∼4
∼3
個数
1
2
8
46
19
表 5: 強度の最小最大幅
RTT-1 のように 1 日周期というような規則性があれば、ルータという物の役割上、恐らく人間
の活動と大きく関わっている・使われているルータと思われる。逆に RTT-207 のようなルータは
まったく人間と関わっていない・使われていないルータなのではないかと推測される。
16
5
OBS
2 章でも述べたが、OBS とは自動車の各部品のパラメータ値を故障時に 15 秒間 ( 0.2 秒毎 ) 採
取した時系列データである。手元には約 14000 件の OBS データがあるが、特にパフォーマンス
系の故障を判別するのが困難であり、それを時系列解析を用いて見出そうというのが本章の目的で
ある。
パフォーマンス系の故障とは、すなわちパラメータの追従関係の崩れであるが、これは大きく分
けて二つに考えられる。すなわち
• 追従しなくなる故障
• 追従しているが、仕方が異常である故障
である。それぞれについてその正常・異常な状態が識別できるか試みた。
5.1
追従しなくなる故障
パラメータには二種類あって、ある部品に対して指令を出しているもの (指令値) と、実際にそ
の部品がどのような値を示しているか (センサー値) を示すものがある。正常であれば図 21 のよ
うに指令値パラメータ (CMD) をセンサー値パラメータ (SENSOR) が一定間隔 k (図では 18) で
追従するようになっている。CMD と SENSOR について相互相関係数 ccf を求めれば、図 22 の
ようにラグ -k (図では -18) で最大値となる。これは SENSOR 値を左 k だけずらせば、CMD と
SENSOR は非常に相関が高いということを示しており、事実と合致している。
図 21: 指令値 (CMD) とセンサー値 (SENSOR)
OBS におけるパラメータ CMP CTRL CMD と CMD CTRL がちょうどそのような指令値と
センサー値の関係になっている。この二つのパラメータを例にとって考えて見る。採取した OBS
はさまざまな故障について集めたものなので、全体で見れば CMD CTRL CMD と CMP CTRL
は正常な関係が多数を占めている。そこで、 CMP CTRL CMD と CMP CTRL をパラメータに
もつ OBS (4009 個) について、これらのパラメータの ccf を求め、相関係数の最大値 (絶対値) を
集計したところ、表 6 のようになった。
標準偏差が 0 で、相関係数が算出できなかったもの 2484 個あるが、これについては後述する。
それらを除外して考えれば相関係数 0.836 以上と高い値を示している obs が 1649 個（算出できた
17
図 22: CMD と SENSOR の ccf
最大相関係数値
OBS 数
0.836 - 1.000
1649
0.673 - 0.836
0.509 - 0.673
0.345 - 0.509
53
29
45
0.182 - 0.345
-0.146 - 0.018
23
1
-0.309 - -0.146
-0.473 - -0.309
-0.637 - -0.473
3
3
9
算出不可 (標準偏差 0)
2484
表 6: CMP CTRL CMD と CMP CTRL の最大相関係数集計表
中では 90 % 以上）ともっとも多く、指令値に対するセンサー値の追従をとらえていると考えられ
る。よって CMP CTRL CMD を CMP CTRL が追従している・いないの閾値を 0.836 としてみ
て、実際の故障例に当てはめてみる。
OBS には故障の識別子が付いている。幸いなことに採取した OBS の中で、 CMP CTRL CMD
に対して CMP CTRL が追従しなくなってしまった故障についてはあらかじめ判明しており、識
別子 P0011 (5 個) と P1009 (119 個) がそれである。表 7, 8 はそれらの識別子をもつ OBS に対
して、同様に相関係数の最大値を集計したものである。
双方合わせて、90 % 以上の場合で、先述の閾値 0.836 よりも低い値となっており、追従して
いるしていないに関しては、相互相関係数の最大値を求める方法が有力であると考えられる。
算出不可 (標準偏差 0) について
どちらかのパラメータ値が 15 秒間全く変化が無いときに、そ
のパラメータの標準偏差が 0 になり、相関係数が算出できない。このような状況は以下の二つに
分けられる。
• 片方のパラメータ値には変化がある
18
最大相関係数値
最大相関係数値
OBS 数
0.222 - 0.297
1
0.146 - 0.222
-0.457 - -0.382
2
1
算出不可 (標準偏差 0)
1
表 7: P0011 の最大相関係数集計表
OBS 数
0.812 - 0.969
8
0.655 - 0.812
0.497 - 0.655
26
15
0.340 - 0.497
0.182 - 0.340
-0.447 - -0.290
35
11
1
-0.605 - -0.447
2
算出不可 (標準偏差 0)
21
表 8: P1009 の最大相関係数集計表
• 両方変化がない
時系列解析では各パラメータに対して parcor 係数というものを求めることができる。これはど
の程度パラメータ値が変動しているかを表しており、これを用いれば、前者と後者を区別すること
が可能である。
前者のパターンで、片方が動いている度合いが高ければ、追従していないすなわち異常であると
考えてよいだろう。後者のパターンでは追従しているかいないかは判断できない。このような場合
はどのような手段を用いても OBS から判断はできない。
5.2
追従の仕方が異常である故障
2 章で紹介した図 3, 4 は、追従しているが仕方が正常と異なる類の故障に含まれる。ここでは
この HO2S S2 と AF SENSOR を例にとって考える。
前節においては追従の有無を検出することが目的であったが、この場合はどれくらいの間隔（ラ
グ）で追従しているかも重要になる。よって前節で ccf により最大相関係数値を求めたが、どれく
らいのラグのときに最大値をとっているかということも集計したのが表 9 である。
相関係数値
OBS 数
( ラグ . OBS 数 ) OBS 数の多い順
0.800 - 1.000
0.600 - 0.800
1457
950
(0 . 959) (-1 . 252) (-2 . 103) (-3 . 56) (-4 . 18) (3 . 10) . . .
(0 . 316) (-1 . 195) (-2 . 109) (-3 . 84) (-4 . 45) (-5 . 24) . . .
0.400 - 0.600
0.200 - 0.400
894
344
(0 . 165) (-1 . 118) (-2 . 102) (-3 . 60) (-4 . 35) (-5 . 34) . . .
(0 . 37) (-1 . 25) (-3 . 22) (-2 . 18) (-5 . 15) (-8 . 10) (-4 . 9). . .
0.000 - 0.200
-0.200 - 0.000
-0.400 - -0.200
15
18
492
(-4 . 2) (-61 . 2) (18 . 2) (-10 . 2) (-25 . 1) (25 . 1) (13 . 1) . . .
(1 . 4) (26 . 3) (-1 . 2) (-29 . 1) (7 . 1) (-6 . 1) (-15 . 1) . . .
(1 . 64) (2 . 62) (4 . 27) (0 . 18) (3 . 18) (5 . 18) (8 . 14) . . .
-0.600 - -0.400
-0.800 - -0.600
962
630
(1 . 152) (2 . 135) (3 . 74) (6 . 43) (4 . 42) (0 . 39) (5 . 33) . . .
(1 . 138) (2 . 102) (3 . 64) (0 . 58) (4 . 32) (5 . 29) (7 . 17) . . .
-1.000 - -0.800
285
(0 . 134) (1 . 45) (2 . 36) (3 . 16) (4 . 12) (-1 . 12) (5 . 6) . . .
算出不可 (標準偏差 0)
2531
表 9: HO2S S2 と AF SENSOR の最大相関係数集計表
19
OBS 全体では正常な状態が多数であることから考えると、比較的高い相関係数値を示している
部分 0.8 - 1.0 や 0.6 - 8.0 では、半分以上がラグ 0 で最大値を示している。従って、正常な状態で
は HO2S S2 はほとんど遅れなしで AF SENSOR を追従していると考えられる。
AF SENSOR と HO2S S2 の関連が崩れることに起因する OBS 識別子も判明しており P2A00
や P1172 である。これらの OBS では相関係数値が軒並み低くなり、ラグが 0 でないものが増え
るのではないかと予想される。表 10, 11 は P2A00, P1172 について同様に最大相関係数値とその
ラグの分布を調べたものである。
相関係数値
OBS 数
( ラグ . OBS 数 ) OBS 数の多い順
0.809 - 0.997
0.622 - 0.809
371
179
(0 . 262) (-1 . 51) (-2 . 44) (-3 . 7) (1 . 3) (3 . 1) (4 . 1) . . .
(-1 . 59) (0 . 49) (-2 . 36) (-3 . 20) (-4 . 5) (12 . 2) (9 . 2). . .
0.435 - 0.622
0.247 - 0.435
48
3
(-1 . 19) (-2 . 9) (-3 . 7) (-4 . 5) (0 . 5) (17 . 2) (31 . 1). . .
(2 . 3)
-0.502 - -0.315
-0.690 - -0.502
-0.877 - -0.690
4
19
7
(7 . 1) (1 . 1) (54 . 1) (33 . 1)
(9 . 3) (43 . 2) (20 . 2) (6 . 2) (17 . 2) (31 . 1) (24 . 1) . . .
(12 . 1) (11 . 1) (6 . 1) (19 . 1) (-56 . 1) (41 . 1) (16 . 1). . .
算出不可 (標準偏差 0)
29
表 10: P2A00 の最大相関係数集計表
相関係数値
OBS 数
( ラグ . OBS 数 ) OBS 数の多い順
0.823 - 0.993
0.653 - 0.823
182
33
(0 . 78) (-1 . 47) (-2 . 29) (-3 . 15) (-4 . 8) (5 . 2) (1 . 1) . . .
(-3 . 10) (-1 . 8) (0 . 5) (-2 . 4) (-5 . 2) (-4 . 2) (7 . 1) (14 . 1)
0.482 - 0.653
0.312 - 0.482
-0.540 - -0.370
8
1
1
(-2 . 2) (-1 . 1) (0 . 1) (20 . 1) (23 . 1) (19 . 1) (-4 . 1)
(21 . 1)
(4 . 1)
-0.711 - -0.540
2
(4 . 1) (16 . 1)
算出不可 (標準偏差 0)
21
表 11: P1172 の最大相関係数集計表
以上のように、OBS 全体で集計した場合と同様で、相関係数値は総じて高く、しかもラグ 0 が
そのうちで最多となってしまっている。そこで実際にいくつか P2A00, P1172 と平常時の OBS の
波形を見比べてみたところ、図 23,24 に示すような違いが見られた。そこでは AF SENSOR に対
して HO2S S2 の反応の遅れよりも、反応度合いの違いの方が顕著である。平常時 (図 23) のもの
は、AF SENSOR の変化幅と同程度に HO2S S2 も変化しているが、P2A00(図 24) のものは、AF
SENSOR の変化幅に対して HO2S S2 の変化幅は鈍化しているように見受けられる。反応の遅れ
はそこまで見られない。
原因として、相互相関関数では値が標準化されるので、このような反応度合いの違いを捉えるこ
とは、できない。
20
図 23: 平常のものと思われる波形と相互相関関数
21
図 24: P2A00 の波形と相互相関関数
22
相互共分散関数
反応度合いの違いを捉える方法の一つとして、値を標準化しない相互共分散関
数を使う方法がある。例えば図 25 上部に示すように、指令値 (CMD) に対して同程度に遅れて反
応している SENSOR-1 と SENSOR-2 があるが、反応度合いが SENSOR-1 に対して SENSOR-2
は半分である。これら CMD と SENSOR-1、CMD と SENSOR-2 に対して相互相関関数を用い
ても、同じ結果しか得られないのだが、相互共分散関数を使えば、図 25 下部に示すように、その
反応度合いに応じて値が変わるので、違いを捉えることができる。しかしこの方法の問題点は、相
互相関関数のように値の大小によってその相関の有無を判断することはできない。
図 25: CMD, SENSOR-1, SENSOR-2 とその相互共分散関数
よって、追従関係の有無およびその遅れを見出すために、最大相関係数値・ラグを利用し、それ
に反応度合いの違いを捉えるために、最大共分散値を併用して、三つの値の正常な状態を算出し、
判断するのがよいのではないだろうか。
23
時系列解析エンジンについて
6
本章ではここまでの時系列解析を実際に行なうために用いてきた時系列解析エンジンについて、
その概要・意義について述べる。
6.1
開発概要
本エンジンは 2009 年 4 月から開発を始め、開発言語は全て Common Lisp を用いている。その
ため基本的に動作プラットフォームに制約されることなく、アプリケーションを使用することがで
きる。たとえば WindowsXP32bit 版 (図 26) や Ubuntu64bit 版 (Linux ベースの OS、図 27) でも
特に手を加えることなく動作させることができる。
処理系はフランツ社製6 の Allegro Common Lisp (以下 ACL) を用いており、その ACL から提
供されている GUI 開発環境によって、ユーザーとのインタラクティブな環境を実現している。こ
れも当然 Common Lisp 準拠であるから前述の如く、プラットフォームを選ばない。
図 26: WindowsXP 上で動作させた画面例
6.2
図 27: LinuxOS 上で動作させた画面例
意義
本エンジンの意義は、開発における柔軟性およびその GUI 環境にある。統計解析エンジンとし
て世界的に広く知られている R7 と比較してみよう。統計解析機能の多様性という意味では、開発
規模・期間からいって本エンジンはまだ R と比較できるようなものではないが、その処理系・開
発環境によって、柔軟にユーザーフレンドリーな機能を追加することができる。一例として、3.1.2
節で挙げたトレンドの推定について本エンジンと R を比較してみる。
R には道具としてトレンドを推定するようなものは揃っているが、一手でトレンド推定を行な
えるようなものは無く、トレンド推定を行なってグラフを描画するまでのコードを書かねばならな
い。例えば、東京都の気温を 16ヶ月にわたって記録したデータを対象に、一次のトレンド推定を
6 http://www.franz.com/
7 http://www.r-project.org/
24
試したところ、以下のようなコードを書くことになった。この結果出力されるグラフが図 28 であ
り、赤線がトレンドを示している。
> temp <- read.csv("tokyo-temperature.csv", header=TRUE)
> temp1<-window(temp,start=c(1,1),end=c(8,1),frequency=30)
> temp2<-window(temp,start=c(8,2),end=c(12,30),frequency=30)
> temp3<-window(temp,start=c(13,1),frequency=30)
> time1<-seq(1:length(temp1))
> time2<-seq(1:length(temp2))
> time3<-seq(1:length(temp3))
> kaiki1<-lm(temp1~time1)
> kaiki2<-lm(temp2~time2)
> kaiki3<-lm(temp3~time3)
> trnd1<-ts(kaiki1$fitted,start=c(1,1),frequency=30)
> trnd2<-ts(kaiki2$fitted,start=c(8,2),frequency=30)
> trnd3<-ts(kaiki3$fitted,start=c(13,1),frequency=30)
> ts.plot(temp,trnd1,trnd2,trnd3,type="l",lty=c(1:4),col=("black","red","red","red"))
一方、開発したエンジンでは、データを読み込み、図 29 に示すような画面で、モデル次数 k お
よび モデル分散 T 2 を入力するだけである。これをコードで書いたとしても、
: (setq temp (time-series-data
(read-data-from-file "sample/tokyo-temperature.sexp")))
: (trend-prediction temp :k 2 :T^2 0.000001)
: (draw (merge-ts temp *))
と簡便であり、図 30 のような、滑らかなトレンド予測ができる。滑らかな予測を R で実現でき
ない訳ではないが、さらに煩雑なコードを書かなければならないであろう。
さらに他の統計解析エンジン、例えば弊社の VMS または S-Plus の時系列解析機能についても
同様で、インタラクティブ環境化されている機能としてトレンド推定はなく、R と同じようにユー
ザーがスクリプトを書かなければならない。VMS, S-Plus にはインタラクティブ環境化されてい
る機能は少なく8 、トレンド推定以外の時系列解析でも、ユーザーがスクリプトを書かなければな
らない状況は免れ得ないであろう。さらに言えばこれらはプラットフォームとして Windows 上で
しか動作しない。
以上のことから、本エンジンのように有用と思われる機能を素早くインタラクティブな環境で実
現し、ユーザーにわかりやすく、動作プラットフォームに拘ることなく使用できることは意義があ
ることであり、強みであると考えられる。
8 自己共分散・相関、ARIMA
モデル、ラグプロット、スペクトル分析の４つ
25
図 28: R によるトレンド推定結果例
図 29: トレンド推定入力画面
図 30: 本エンジンによるトレンド推定結果例
26
7
まとめ
本レポートでは、時系列解析エンジンを用いて、3 つの時系列データを時系列解析の俎上に乗
せ、それによってどのようなことができるかを主に述べてきた。
まず初めにイギリスのガス使用量の記録 UKgas を例に用いて、時系列解析の凡例を以下の手順
で示した。
1. 記述 (acf, ccf, トレンドの推定, 移動平均 など)
2. モデリング (AR モデル, HoltWinters モデルなど)
3. 予測
そして、実際の時系列データにおける時系列解析の応用を、時系列データ RTT と OBS につい
て試み、以下のような結果となった。
• RTT
ルータの ping の遅れを表す時系列データ。移動平均をとることでその周期性が洞察され、acf
によって大半のルータは 1 日周期を持つことが判明した。ピリオドグラムを求めることで周
期成分のない白色雑音のようなルータがいくつかあることが判明し、ほとんど使用されてい
ないルータなのではないかと推測される。
• OBS
自動車の各部品の車載時系列データ。パフォーマンス系の故障について ccf を用いることで
識別を試みた。結果、追従の有無に関する故障については高い確率で識別できることが判明
したが、追従はあるが、通常と異なる故障に関しては識別できなかった。原因として、この
系統の故障は、追従の遅れよりも、追従の反応度合いがより顕著に変化し、それは ccf では
評価できないことにあった。対策としては自己共分散関数を用いることが考えられる。
現状の時系列解析エンジンでは、時系列データの特徴を記述する基本的な道具に関してはほとん
ど備わったと考えており、今後はモデリングの充実を図る。今回示したように、時系列データには
さまざまなパターンがあるので、多様なモデルを構成できたり、その検定・推定などができるとい
うことは時系列解析エンジンの性能と直結する。具体的な方向性としては、非ガウス型・非定常型
時系列モデルを考えている。実際の時系列データはガウス分布でないもの・定常的でないものが多
くあると考えられるからである。
参考文献
[1] 北川 源四郎 『時系列解析入門』
(岩波書店、2005)
[2] 金 明哲 『R によるデータサイエンス』
(森北出版、2007)
[3] 田中 孝文 『R による時系列分析入門』
(シーエーピー出版、2008)
27
[4] 松原 望、縄田 和満、中井 検裕 『基礎統計学 I 統計学入門』
(東京大学出版会、1991)
[5] Brockwell, Davis “Introduction to Time Series and Forecasting, Second Edition”
(Springer, 2001)
[6] James D. Hamilton “Time Series Analysis”
(Princeton University Press, 1992)
[7] 『WIKIPEDIA』
http://ja.wikipedia.org/wiki/
28