概要
次のような時系列データがあります。
上図のような時系列データに対して、Pan-Matrix Profileは次のようになります。
下部のヒートマップがPan-Matrix Profileで、横軸が時間、縦軸が部分時系列の長さになります。★で示す頂点（？）を参照することで、時系列に内在する頻出パターンの出現位置(横軸)とその長さ(縦軸)が一目で分かります。
Pan-Matrix Profileとは
前回、時系列データ解析の革新的技術としてMatrixProfileについて紹介しました。Pan-Matrix Profile(PMP)は、簡潔に言えばMatrixProfileの部分時系列長$m$をある範囲で全て調べ、行列としたものです。
先ほどの画像の例を見てみましょう。
Pan-Matrix Profileの行列の各要素が何を意味するのか説明します。例として、(504,80)の要素を考えます。先ほど述べた通り、この行列は横軸が時間、縦軸が部分時系列の長さを表します。まず、時系列の504点目～504+80点目の部分時系列(上側の図の青色の部分時系列の左の方)について、その部分時系列と最も似ている、同じ長さの部分時系列(上図、青色、右の方)との距離を計算します。ここで、距離の尺度はピアソンの相関係数を用います。今、その距離の値が0.95なので、(504,80)の要素には0.95が入ります。
さて、この行列を力任せ探索で作ろうとすれば実に$O(n^4)$(ただし$n$は時系列長)とかいう計算時間がかかる訳ですが、文献[1]では様々な高速化テクニックにより$O(n^2r)$(ただし$r$は部分時系列長の候補数)まで削減されます。
#Pan-MatrixProfile用ライブラリ"matrixprofile"
matrixprofile用のライブラリはたくさんあって迷ってしまいますが、今回はmatrixprofileという名のPython用ライブラリを使います。そのまんまですね。
pip install matrixprofile
実装例①：人工データへの適用
まずは、ライブラリに標準でついてくるいくつかの時系列データを合成された人工データで実装してみましょう。
from matplotlib import pyplot as plt
import numpy as np
import matrixprofile as mp
# データの読み込み
dataset = mp.datasets.load('motifs-discords-small')
X = dataset['data']
# データの可視化
plt.figure(figsize=(18.0, 6.0))
plt.plot(np.arange(len(X)), X, color="k")
plt.xlim(0, len(X))
plt.title('Synthetic Time Series')
Pan-Matrix Profileを作成し、可視化するコードは次のようになります。
#Pan-Matrix Profileおよびその他の解析
profile, figures = mp.analyze(X)
# PMP表示
figures[0]
mp.analyze(X)で簡単にPMPが作成出来ます。profile内にはPMPの行列の他、PMP作成時に必要な各部分時系列に対する最近傍部分時系列のインデックスの行列(PMPI)、MotifやDiscord等の解析情報やその他もろもろが入っています。
そして、figuresには解析を可視化した画像が入っており、mp.analyze(X)の実行時点で表示が実行されますが、figures[0]でPMPの部分のみを表示することが出来ます。
実装例②：ミトコンドリアのDNA配列への適用
続いて実データへ適用してみましょう。元論文[1]で紹介されていたミトコンドリアのDNA配列へ適用します。（こちらは厳密には時系列ではないですが…。）
データはこちらからmat形式のものが入手できます。
#データの読み込みその２
import scipy.io as sio
dataset = sio.loadmat('termite_DNA_circular_shift')
X = dataset['t2'].reshape((-1,))
#データの可視化
plt.figure(figsize=(18.0, 6.0))
plt.plot(np.arange(len(X)), X, color="k")
plt.xlim(0, len(X))
plt.title('Mitochondrial DNA Sequence')
Pan-Matrix Profileを作成してみましょう。今回は数分かかります。
#Pan-Matrix Profileおよびその他の解析その２
profile, figures = mp.analyze(X)
#PMP表示その２
figures[0]
PMPが無駄に大きいような気がしますね。実は、mp.analyzeには閾値パラメータthresholdがあり、そちらを調整することで部分時系列長上限値の探索範囲を調整出来ます。パラメータフリーの文言はどこ行ったって思いますが、まぁ部分時系列長$m$ほど敏感ではないらしいです。
実装例③：眼電図(EOG)への適用
こちらも元論文[1]にて用いられていたEOGのデータです。
#データの読み込み３
dataset = sio.loadmat('eog_multiple_scale_example')
X = dataset['testdata'].reshape((-1,))
#データの可視化
plt.figure(figsize=(18.0, 6.0))
plt.plot(np.arange(len(X)), X, color="k")
plt.xlim(0, len(X))
plt.title('EOG')
デフォルトのままで可視化すると以下のようになります。
#Pan-Matrix Profileおよびその他の解析その３
profile, figures = mp.analyze(X)
#PMP表示その３
figures[0]
さて、今回は非常に縦幅の小さなPMPとなりました。拡大して、適当にMotifを取り出してみます。
さて、一番長いMotifでこんな形です。元論文で紹介されていたこのデータのMotifはもっと長かったはずです。一体何が問題なのでしょうか。
まず、matrixprofileライブラリの可視化部分のコードを覗いてみると、PMPにおける距離が1以上のものは全て1で埋められています。つまり、PMPの画像において、黄色い部分は全て1です。こんなことしちゃってよいのでしょか。
直観的に、部分時系列間の距離は、部分時系列長が長くなれば長くなるほど大きくなる傾向にあるように考えられます。つまり、この方法では長いMotifを取り出すことが出来ません。
ちなみに、元論文によると、（理解出来ているか自信がないのですが）Top-k個のMotifを取り出す方法はPMPの値が小さいものを取り出しているっぽいです。matrixprofileライブラリのMotif抽出コードもそのようになっています。これでは、短いMotifほど過大に評価されると思うのですが・・・。
↑2022年9月追記
論文をよく読むとPMPを計算するときは時系列間の距離の尺度として相関係数を用いるのですが、この実装ではz-ユークリッド距離になってしまいます。実装が良くなかったです。言い訳は長くなるので省略します。
まとめ
という訳で、今回は時系列データ解析の最新手法Pan-Matrix Profileの紹介でした。 どうも疑問の残る終わり方となってしまいましたが、そのうちこの疑問解決に取り組んでみたいと思います。 ←解決しました。
参考文献とか
- [1]Madrid et. al. Matrix Profile XX: Finding and Visualizing Time Series Motifs of All Lengths using the Matrix Profile. IEEE Big Knowledge 2019.
- UCRのMatrixProfileサイト：https://www.cs.ucr.edu/~eamonn/MatrixProfile.html
- PMPのライブラリ：https://matrixprofile.docs.matrixprofile.org/installation/Linux_Installation.html
- データセット：https://sites.google.com/view/pan-matrix-profile/datasets?authuser=0