放電加工機の時系列解析（その４）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 プログラムについて粒子フィルターを実行してくれるライブラリーを見つけることはできなかったです。
参照したプログラムは、参考文献にありますが、Class を 定義して、 量子フィルターを実行するプログラムを構築されています。
以下の関数を使用しています。
scipy.stats import multivariate_normal, cauchy, invgamma
数学・科学・工学のための数値解析ソフトウェア モジュール
統計機能関連の確率変数等のサブパッケージです。
多変量正規分布に従う乱数
コーシー分布
逆ガンマ関数from numpy.random import uniform
乱数を発生する関数

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

～参考文献、WEB～
時系列解析 -自己回帰モデル・状態空間モデル・異常検知- 島田直希 著 共立出版（株）
Python によるパーティクルフィルタの実装と状態空間モデルへの適用
https://qiita.com/kenmatsu4/items/c5232b1499dfd00e877d32. 非線形非ガウス型モデルと粒子近似
今回は、状態の確率分布が非ガウス性があり、状態と観測値との間に線形の関係が無いとした場合の近似に取り組みます。
2.1 非ガウス型とは （ガウス分布と正規分布は同じです）
下の分布図で、左が、ガウス分布（平均=0,標準偏差=1）であり、右が非ガウス分布です。
右図は 平均=-2,標準偏差=1 と平均=2,標準偏差=0.8 のガウス分布を、重ねただけで、少し変ですが、2 峰であり非ガウス分布となります。
ガウス分布

非ガウス分布

下図のような時系列データを考えるとき、t-1 時点のの確率分布をガウス分布、非ガウス分布のいずれで、解析するかということです。
時系列データ

X
x-1

t-1tT

ガウス分布であれば、平均と分散であらわすので、解析も扱いやすです。分布が非ガウス分布の場合（実際には時系列データはガウス分布ではないことが多いと思います）、なんらの近似が必要で、３つの方法が考えられます。
・混合ガウス近似 いくつかのガウス分布の和として近似します。
・階段関数近似階段関数（ヒストグラム）として近似します。

・粒子近似

分布を粒子の数として近似します。(粒子フィルターと呼ばれます)

今回は、粒子近似に取り組みます。
4

2.2 粒子フィルタ（partticle filter）とは2.2.1 概要確率分布を点の集まりと考えます。確率が高いところは粒子の密度を高くします。

粒子フィルターは、任意の分布を、に示す粒子の密度に近似して解析します。
2.2.2 量子フィルターの実際
まずは、量子フィルターを、数式で表現します。
𝑥𝑡 = 𝐹𝑡 ( 𝑥𝑡−1, 𝑣𝑡)

𝑣𝑡:システムノイズ𝑦𝑡 = 𝐻𝑡 ( 𝑥𝑡−1, 𝑤𝑡)

𝑤𝑡：観測ノイズ

量子フィルタの仕組みを下図で説明します。次ページより 1 フレーム毎、説明します。
量子数は説明のため、6 個にしてます。
（実際は 1000 個以上が推奨されています）

5

STEP1：
t 時点を、量子フィルターで近似しますので、t-1 時点が初期値となります。初期値は、その 1 個前がないので、乱数を生成して作成します。
y は観測値で、x は状態値になります。

上述した数式 𝑥𝑡 = 𝐹𝑡 ( 𝑥𝑡−1, 𝑣𝑡)

により、t 時点の 𝑥𝑡 値を計算します。
システムノイズ 𝑣𝑡は ガウス分布とします。

STEP3:各粒子の 𝑥𝑡が得られるので、時刻 t の観測値𝑦𝑡 から 各粒子の尤度を計算します。
尤度を重みとして、計算します。
円の大きさでその重みを表しています。

6

STEP4:重みを使って加重平均したのが、緑の丸で、フィルタリングした 𝑥𝑡の値です。

STEP5:
尤度の低い粒子は消滅させて、尤度の高い粒子は複数個に分裂し、新たな粒子の確率分布として、扱います。

STEP1 から STEP5 を繰り返して、推定、長期予測を計算します。
プログラムには、前述しましたが、定型のライブラーは無いようです。
参考文献を参照ください。

二つの Class が定義されており、ガウス分布を持つノイズ と 確率密度 と システムノイズの更新 を行う Class,上図のリサンプリングまでの計算 と 推定と 長期予測 を行う Class です。

7

モデルのあてこみとして、以下の入力を行います。
（）内が今回の設定です。
・観測値の次元数(1)

・トレンド成分の次元数(1)・季節成分の次元数(3)・AR 成分の次元数(2)・粒子数(1000)・システムノイズの分散共分散の各要素の値 (1500)
・観測ノイズの分散共分散の各要素の値 (10)・システム関数 ：transition_func = lambda x, F:F.dot(x) と定義

・観測関数

observation_func = lambda x, H:H.dot(x) と定義：

・各ノイズオブジェクトのインスタンスは Class から計算します。
2.3 予測の結果
2.3.1 観測値と推定値のグラフと成分分解モデル

青色が観測値で、橙色が推定値です。
線形ガウス型モデルのカルマンフィルターでの推定値は、観測値と一致していましたが、粒子フィルタは、傾向は同じですが、一致まではいきませんでした。
次ページに、トレンド成分、季節成分、AR 成分に分けたグラフを提示します。
青色：観測値、橙色：季節成分、 緑色：季節成分、 赤色：AR 成分です。

8

2.3.2 予測のグラフ （333 以降を予測しました。
）333 から少しの間は、予測してるように見えますが、観測値からの推定が、できないためか、中心値に収束しています。

93 まとめ
いままで 2 回計算したグラフと比較します。いずれも、先の予測ができてないです。
今後も、色々と取り組む予定です。

3.1 線形ガウス分布の状態空間モデルで計算
推定はカルマンフィルタ

300

333

3.3 線形ガウス分布の自己回帰モデル
で計算 推定は最尤法

10364