放電加工機の時系列解析（その３）
1 はじめに
1.1 時系列の取組について放電加工機の時系列解析に取り組みます。
時系列解析は、以下の４ステップから構成されます。
① ディスクリプション (description)
時系列の特徴を把握します。定常性の確認、成分分解等を行います。
② モデリング（modeling）
時系列モデルを構成し、パラメータを推定します。
③ 予測(prdiction)

現在までの情報から今後の変動を予測します。
④ 信号抽出必要な信号や情報を取り出す。
（異常検知とか）
今回の報告書は、上記の モデリング（modeling）と 予測 (prdiction)についての取組内容です。モデルについて、代表的な手法として、自己回帰モデル と 状態空間モデルがありますが、今回は、状態空間モデルについての報告です。

1.2 どのようなデータを扱ったか
放電加工機は、加工中に障害が発生すると自己診断して 3 段階のレベルに分け、レベルに応じたメッセージを表示します。内容は以下です。
・エラーメッセージ：続行不可能な障害が発生したときに表示し、動作を中断する。
・ハルトメッセージ：再開可能な障害が発生したときに表示し、一時停止する。
・コメントメッセージ：続行可能な障害で発生し、注意を促す。
また、メッセージの記録は、USB により CSV データとして取り出すことが可能です。
レベルに対応してデータ処理できるように、3 段階のレベルについて障害のレベルが大きいと、点数は大きくなるように点数を設けました。
これらは、メッセージが発生した時の時間毎のデータで、点過程データと呼ばれるもので、
1

時系列データではありません。故に、一日分の平均値を取って時系列データとして取り出すようしました。本データを以下に示します。2022 年 4 月 から 2023 年３月までのデー
タです。
上段が発生時毎のデータ（点過程データ）で、1 年で 37000 回発生してます。
1 点の発生は、一定時間の間隔ではありません。
下段は点過程データを、一日毎に平均して、平均したデータです。点の間隔は 1 日毎ですので、時系列データと言えます。

2

1.3 プログラムについて今回 以下三つのパッケージで計算して、予測しています。
① R の 時系列解析パッケージ TSSS で計算します。
１変量 AR モデルのあてはめ → arfit
状態空間モデルによる時系列の予測と補間 →https://jasp.ism.ac.jp/ism/TSSS/

tsmoothを参照ください。

② Python の statsmodels.tsa.statespace で計算します。
モデルのあてはめ→

UnobservedComponents()

但し ローカルレベルモデル（次章で説明します）に対応したものです。

③ Python の Pykalman と scipy.optimize の minimize 関数、で計算します。
KalmanFilter

計算 →パラメータの推定 →Pykalmanscipy は数値解析ソフトウェア・ライブラリーでoptimize 最適化を求めるツールで,minimize は最小化を行います。

本報告書にはプログラムは記載していません。

～参考文献、WEB～「R による時系列モデリング入門」 北川源四郎 著 （株）岩波書店
時系列解析 -自己回帰モデル・状態空間モデル・異常検知- 島田直希 著 共立出版（株）第 3 回 時系列分析～状態空間モデル編～
https://note.com/yiida/n/n5e5db1bea552

3

2 状態空間モデルについて先に報告した自己回帰モデルによる解析は、観測値を直接モデル化しました。
状態空間モデルは、状態と観測値に分解してモデル化する手法で、単一のモデルというより、様々なモデル解析を統合したフレームワークといえます。
2.1 状態空間モデルとは状態値と観測値の関連を図で表示しました。X が状態値で、Y が観測値となり、状態値は1 時点前の状態に影響され、観測値は、その時点の状態値に影響されます。

さらに、数式では以下となります。
を l 変量の時系列として、以下を状態空間モデルと呼びます。

𝑦𝑡

𝑥𝑡= 𝐹𝑡 𝑥𝑡−1+ 𝐺𝑡 𝑣𝑡(状態モデル)𝑦𝑡= 𝐻𝑡 𝑥𝑡+ 𝑤𝑡（観測モデル）𝑥𝑡 :観測できない k 次元のベクトルで、状態と呼びます。

𝑣𝑡 :状態ノイズで、平均ベクトル０，分散共分散行列𝑄𝑡 に従う m 次元の正規白色雑音

𝑦𝑡 :観測可能な時系列データで、l 次元ベクトル

𝑤𝑡 :観測ノイズで、平均ベクトル０，分散共分散行列𝑅𝑡 に従う l 次元の正規白色雑音

𝐹𝑡

𝐺𝑡

𝐻𝑡 : それぞれ 係数行列(k x k),係数行列(k x m),観測行列(l x k) の行列です。

状態空間モデルは、非定常、非線形でもモデル化が可能とされています。
次に、状態空間モデルで Local Level モデルというものがあり、以下の数式となります。
𝑥𝑡= 𝑥𝑡−1+ 𝑣𝑡(状態モデル)𝑦𝑡= 𝑥𝑡+ 𝑤𝑡

（観測モデル）係数を１として、計算を簡略化しています。
（statsmodels.tsa.statespace に使用しているモデルです。）
4

2.2 状態の逐次推定とパラメータ推定推定について述べますが、線形ガウス型モデルについての内容です。
状態の逐次推定は、以下の図で説明されます。

状態𝑥𝑡 を推定することを状態推定と呼び、観測時点 j と・1 期先予測 : j < tp(𝑥𝑡 | 𝑌1:𝑗 )・フィルタ : j = tt の関係から、三つに分けます。

・平滑化 : j > tとは、時点１からｊまでの観測データ集合Ｙから、t 時のｘの確率分布を求めることを意味し、カルマンフィルタというアルゴリズムを用います。
また、p(𝑥𝑡 | 𝑌1:𝑗 )は正規分布に従うので、平均𝑥𝑡|𝑗 と 分散共分散行列とになります。
𝑥𝑡|𝑗= E[𝑥𝑡 | 𝑌1:𝑗 ]𝑉𝑡|𝑗= E[(𝑥𝑡 − 𝑥𝑡|𝑗 ) (𝑥𝑡 − 𝑥𝑡|𝑗 )𝑇 | 𝑌1:𝑗 ]となり、カルマンフィルタにより以下の計算式となります。
・一期先予測・フィルタ

𝑥𝑡|𝑡−1 =𝐹𝑡 𝑥𝑡−1|𝑡−1

𝑉𝑡|𝑡−1= 𝐹𝑡 𝑉𝑡−1|𝑡−1 𝐹𝑡𝑇+𝐺𝑡 𝑄𝑡 𝐺𝑡𝑇

𝐾𝑡 = 𝑉𝑡|𝑡−1 𝐻𝑡𝑇 (𝐻𝑡 𝑉𝑡|𝑡−1 𝐻𝑡𝑇 + 𝑅𝑡 )𝑇𝑥𝑡|𝑡= 𝑥𝑡|𝑡−1 + 𝐾𝑡 (𝑦𝑡 - 𝐻𝑡 𝑥𝑡|𝑡−1 )𝑉𝑡|𝑡=( 𝐼𝑘- 𝐾𝑡 𝐻𝑡 ) 𝑉𝑡|𝑡−1

I は単位行列・平滑化−1𝑇
𝐴𝑡 = 𝑉𝑡|𝑡 𝐹𝑡+1
𝑉𝑛+1|𝑛𝑥𝑛|𝑁 = 𝑥𝑛|𝑁 + 𝐴𝑛 (𝑥𝑛+1|𝑁5- 𝑥𝑛+1|𝑛 )を求めるこ𝑉𝑛|𝑁 = 𝑉𝑛|𝑛・長期予測+ 𝐴𝑛 (𝑉𝑛+1|𝑁- 𝑉𝑛+1|𝑛 )𝐴𝑇𝑛i = 1,…,j について、
𝑥𝑛+𝑖|𝑛 =𝐹𝑛+1 𝑥𝑛+𝑖−1|𝑛𝑇
𝑉𝑛+𝑖|𝑛 = 𝐹𝑛+1 𝑉𝑛+𝑖−1|𝑛 𝐹𝑡+1

𝑇+ 𝐺𝑛+1 𝑄𝑛+𝑖 𝐺𝑡+1

これらの計算式で、求めていくわけですが、F,G,H,Q を決めないと計算できません。
（仮に、想定した数列で計算を進めることもあります）
F,G,H,Q は、パラメータの推定により、 を計算していきます。
パラメータの推定について、先回の技術報告書に、数式として紹介しました。
今回も繰り返しになります。
パラメータは θ = (𝑎1 , ………… 𝑎𝑚 , 𝜎 2
対数尤度は、)𝑇とします。

l(θ) = ∑𝑀𝑛=1 𝑙𝑜𝑔𝑝( 𝑦𝑛 |𝑦1 ,…, 𝑦𝑛−1 )となります。

さらに、l(θ) を変換すると、 （参考文献を参照ください）𝑚2l(θ) = -((N-M)/2)log2π𝜎 2 - 1/2𝜎 2 ∑𝑁
𝑛=𝑀+1( 𝑦𝑛 - ∑𝑖=1 𝑎𝑖 𝑦𝑛−𝑖 )となり、l(θ) を 最大とする、パラメーターを算出します。
まずは、l(θ) を最大にする分散 𝜎 2 を求めると、（参考文献を参照ください）1

𝑚2𝜎 2 = 𝑁−𝑀 ∑𝑁
𝑛=𝑀+1( 𝑦𝑛 - ∑𝑖=1 𝑎𝑖 𝑦𝑛−𝑖 )自己回帰係数 𝑎1 , … 𝑎𝑚となりの対数尤度はl(𝑎1 , … 𝑎𝑚 ) = -((N-M)/2)log2π𝜎 2 - (N-M)/2 となります。
分散 𝜎 2を最小化する 自己回帰係数を求めることになります。

今回は、①Ｒの TSSS 解析パッケージにおいては、ユール-ウォーカー法、で計算しています。
②Python の statsmodels.tsa.statespace においては、そもそも、係数は１です。
③Python のscipy.optimize の minimize 関数 においては、最適化手法として、SLSQP法を、計算しています。

6

2.3 予測の結果
2.3.1 R 時系列解析パッケージ TSSS での予測観測モデルを AR モデルに従うとして、状態空間モデルを予測してます。
𝑦𝑛 = ∑𝑚𝑖=1 𝑎𝑖 𝑦𝑛−𝑖+ 𝑣𝑛（観測モデル） となり、状態モデルを、𝑥𝑛 = (𝑦𝑛 , 𝑦𝑛−1 , … 𝑦𝑛−𝑚+1 )𝑥𝑛= F𝑥𝑛−1 +𝑦𝑛= H𝑥𝑛

G𝑣𝑛𝑇と定義すると、という関係になりで観測モデルが得られます。

本モデルは、状態モデルで、観測値で、完全に決定され、観測ノイズは０となるモデルと定義されています。この条件で、TSSS を用いて、計算しプロットしたグラフが以下です。
出力されたグラフです。線を太くしたかったのですが、線太さの変更が R では難しく、デフォルトで出力してます。黒線が、観測値で、赤線が予測です（3 月分で x 軸で 335 から365 の間）。また、青線の幅が、標準予測誤差区間です。

m=30m=15m=23357ｍは次数です。

m=2, m=15, m=30 と変化させてグラフを作成しました。

以下は、次数毎の対数尤度と AIC の値です。
この表からは、m=2 が最適となりました。

対数尤度

AIC

m=2-142.294-278.689（対数尤度は、大きい数値、AIC はm=15-147.128-262.255

小さい数値の方が、最適と考えます）

m=30-151.663-241.326

予測したグラフからは、335 から 337 までは予測できてますが、それ以降は、予測精度は、落ちてしますが、m が大きい程、振幅の減衰は少なく、より良く予測できてるように見えます。

2.3.2

Python statsmodels.tsa.statespace での予測

拡大

青線が観測値、赤線が推定値、黒線が予測値 です。
予測値は、観測値の振幅に追従することはなく、一定値で推移しています。
Local Level モデルで、係数を１にしているので、予測値精度は落ちているようです。

8

2.3.3 Python Pykalman, scipy.optimize,での予測観測モデルを、トレンド成分（時間とともに単調に増加/減少する変動）と季節成分（同じ周期で規則的に繰り返される変動）と不規則成分との合計と考え、不規則成分は、AR モデルに準じると考えます。
数式で示すと= 𝑡𝑛 + 𝑠𝑛 + 𝑟𝑛 + 𝑤𝑛

𝑦𝑛𝑤𝑛𝑡𝑛(𝑘)= ∑𝑘𝑖=1 𝑐𝑖 𝑡𝑛−𝑖 + 𝒗𝒏𝟏

𝑠𝑛= - ∑𝑖=1 𝑠𝑛−𝑖𝑟𝑛=𝑣𝑛1 ~ N ( 0, 𝜏12 )トレンド成分+ 𝒗𝒏𝟐

𝑣𝑛2 ~ N ( 0, 𝜏22 )

季節成分∑𝑞𝑖=1 𝛷𝑖 𝑟𝑛−1 + 𝒗𝒏𝟑

𝑣𝑛3 ~ N ( 0, 𝜏32 )AR 成分

𝑝−1となり、状態モデル 𝑥𝑛
𝑥𝑛~ N ( 0, 𝛿 2 )は次のようなります。

= (𝑡𝑛 , 𝑡𝑛−1 ,… 𝑡𝑛−𝑘+1 , 𝑠𝑛 , 𝑠𝑛−1 ,… 𝑠𝑛−𝑝+2 , 𝑟𝑛 , … 𝑟𝑛−𝑞+1 )𝑇

2-1 項にある状態空間モデルの式で、進めます。
𝑥𝑛= 𝐹𝑛 𝑥𝑛−1+ 𝐺𝑛 𝑣𝑛(状態モデル)𝑣𝑛~ N( 0, 𝑄𝑛 )𝑦𝑛= 𝐻𝑛 𝑥𝑛+ 𝑤𝑛（観測モデル）𝑤𝑛~N( 0, 𝑅𝑛 )ここで、トレンド成分の係数行列を

𝐹1 , 𝐺1 , 𝐻1 とし,季節成分の係数行列を

𝐹2 , 𝐺2 ,𝐻2 とし, AR 成分の係数行列を 𝐹3 , 𝐺3 , 𝐻3 とすると、
𝐹𝑛= 𝐹10 00 𝐹200 0 𝐹3𝐺𝑛 = 𝐺10 00 𝐺2𝐺3𝐻 𝑇 = 𝐻1𝑇0𝐻2𝑇0 0𝐻3𝑇で表現されるので、データの特徴から 行列を定義します。
（パラメータの推定は未だです）
トレンド成分は、傾きは無いとして、次数を１としました。係数は１です。
季節成分は、月単位の変化を想定して、次数を 30 としました。係数は１です。
AR 成分は、次数を２としました。係数は 0.5 です。

9この条件で、Pykalman で計算して、グラフを作成しました。
青線が、観測値で、橙色が推定値で、333 から 364 は予測値になります。
（4 月 1 日を０から起算しているので、1 日引いた数字になってます）

青線と橙線は、ほぼ重なり
訓練データへのあてはめは出来てますが、予測は外れてしまいました。

300

333

10364

さらに scipy.optimize の minimize 関数を用いて、係数の推定を行い、グラフを作成しました。

300

333

訓練データへのあてはめは出来てますが、予測は、分布の平均値であり、振幅は、追従できてません。

11364

グラフを成分別に分けたグラフが以下です。

下から、季節成分（緑線）、AR 成分（赤線）、トレンド成分（橙線） となります。

3. まとめと今後の取組
今回は、①R の TSSS、 ②Python の statsmodels.tsa.statespace、③Python の Pykalman と scipy.optimize による計算で、いずれも、係数を推定し、カルマンフィルタにより推移を推定しますが、いずれも、提供されるパッケージにある、ブラックボックスと言える関数による計算です。
①と③は、訓練データのあてはめは出来たのですが、予測がうまくできませんでした。
②は、ローカルレベルモデルのためか、①、③ と比較して、精度は劣るようです。
今後、モデルの定義、係数の見直し等に、取り組む予定です。

12