時系列分析入門【第５章 その１】２つの時系列データの相互相関・移動相関・動的時間伸縮をPythonで実践する
この記事は、テキスト「RとStanではじめる 心理学のための時系列分析入門」の第５章「時系列データ同士の関係の評価」のRスクリプトをお借りして、Python で「実験的」に実装する様子を描いた統計ドキュメンタリーです。
取り扱いテーマは、２つの時系列データの類似度に関する次の指標です。
・相互相関／移動相関
・動的時間伸縮法（DTW）
時系列データ間の距離感って、どんなことなのでしょう？
さあ、時系列分析を楽しんでいきましょう！
テキストの紹介、引用表記、シリーズまえがきは、このリンクの記事をご参照ください。
テキストで使用するデータは、R・Stan等のサンプルスクリプトで指定されています。
サンプルスクリプトは著者GitHubサイトからダウンロードして取得できます。
第５章 時系列データ同士の関係の評価
この記事は「5.2 相互相関」「5.3 動的時間伸縮法」を実践いたします。
インポート
この章で用いるライブラリをインポートします。
### インポート
# 基本
import numpy as np
import pandas as pd
# Rデータセット読み取り
import statsmodels.api as sm
import rdata
# WEBアクセス／googleトレンド
import requests
from pytrends.request import TrendReq
# Zipファイル操作
import zipfile
# データ前処理
from sklearn.preprocessing import StandardScaler
# 日付処理
from datetime import datetime, timedelta
from dateutil import tz, parser
# DTW処理
# pip install dtw_python
from dtw import dtw
# 信号処理
from scipy.signal import spectrogram, coherence, periodogram
import scipy.signal as signal
# 方向データの統計処理
# pip install pingouin
import pingouin as pg
# 階層的クラスタリング、デンドロビウム描画
from sklearn.cluster import AgglomerativeClustering
from scipy.cluster.hierarchy import dendrogram
# VAR・AR
from statsmodels.tsa.vector_ar.var_model import VAR
from statsmodels.tsa.ar_model import AutoReg
# グラフ描画
import matplotlib.pyplot as plt
import seaborn as sns
plt.rcParams['font.family'] = 'Meiryo'
# ワーニング表示の抑制
import warnings
warnings.simplefilter('ignore')【補足】パッケージのインストール
利用環境に応じてパッケージをインストールしてください。
以下はインストールの参考コードです。
### インストール
## dtw-python
# anaconda
conda install -c conda-forge dtw-python
# pip
pip install dtw-python
## pingouin
# anaconda
conda install -c conda-forge pingouin
# pip
pip install pingouin
相互相関と移動相関
① データの読み込み
RのrMEAパッケージに含まれる「臨床対話中のクライアントとセラピストの身体動作」に関するデータセットをWebサイトより読み込みます。
10分間の臨床対話のビデオより、クライアントとセラピストの動きをMotion Energy Analisysで評価した値が記録されています。
### txtデータの取得
# データをダウンロードしてpandasのデータフレームに格納
url = r'https://raw.githubusercontent.com/kleinbub/rMEA/master/inst/extdata/' \
'normal/203_01.txt'
df_cor = pd.read_csv(url, sep=' ')【実行結果】
１秒あたり25フレーム、全15000フレームのデータです。
patはクライアントの動き値、thはセラピストの動き値です。
データを描画します。
テキストの図5.1左に相当します。
pandasのplotを利用します。
### データの描画 ※図5.1（左）に相当
df_cor.plot(figsize=(10, 5), xlabel='Time', ylabel='動きの大きさ', legend=False,
color=['black', 'red'], lw=0.5, ls='dashed')
plt.legend(['クライアント', 'セラピスト'], loc='upper right');【実行結果】
② 相互相関の算出と描画
描画ライブラリである matplotlib.pyplot に搭載された、相互相関を算出して描画する xcorr() を利用します。
テキストの図5.1右に相当します。
### 相互相関 matplotlibのxcorrを利用 ※図5.1（右）に相当
# データの標準化
df_cor_std = (df_cor - df_cor.mean()) / df_cor.std()
# 相互相関の描画（同時に相互相関関数をxcorr_resへ格納）
xcorr_res = plt.xcorr(x='pat', y='th', data=df_cor_std, normed=True,
maxlags=100)
# 修飾
plt.xlabel('Lag')
plt.ylabel('相互相関関数')
plt.show()【実行結果】
テキストによると、ラグ０付近でクライアントとセラピストの動作に正の相関があるので両者が体を動かすタイミングは類似している、ということのようです。
また、ラグ－５０：－２秒あたりでは負の値になっており、セラピストの動きがクライアントの動きに先行しているようです。
変数 xcorr_res に代入した相互相関の値を確認してみましょう。
### 相互相関関数の表示
xcorr_res[1]【実行結果】
表示は途中省略されています。
計算値を保持できることがわかりました。
③ 移動相関
一定期間を区切る「窓」を移動させて相関をはかる「移動相関」を算出します。
窓長1000（40秒）、ラグ50（2秒）で計算します。
移動相関を算出できるPythonパッケージが見当たらなかったので、matplotlib.pyplot の xcorr() を用いた計算関数を作成しました。
### 移動相互相関の計算関数の定義
def calc_moving_xcorr(df):
# 初期値設定
lag = 50 # ラグの長さ
window = 1000 # 窓の長さ
fps = 25 # 1秒間の映像コマ数：25フレーム/秒
xcorr_list = [] # 一時リスト
# データの標準化
df_std = (df - df.mean()) / df.std()
# 窓を移動させて指定ラグ長の相互相関関数を算出
for i, win in enumerate(df_std.rolling(window=window)):
# 期間が窓未満の場合は処理しない
if i < window-1:
pass
else:
# matplotlibのxcorrで窓ごとの相互相関関数を計算して一時リストに格納
xcorr_res = plt.xcorr(x=win['pat'], y=win['th'], maxlags=lag)
xcorr_list.append(xcorr_res[1])
# xcorrの描画を表示しない指定
plt.close()
## データフレーム化
moving_corr = pd.DataFrame(xcorr_list).T
# 列名を分単位に補正
moving_corr.columns = moving_corr.columns / (60 * fps)
# インデックスを秒単位に補正
moving_corr.index = (moving_corr.index - lag) / fps
# 行を逆順にソート
moving_corr = moving_corr.iloc[::-1, :]
# 戻り値：移動ラグ相関のデータフレーム
return moving_corr【実行結果】なし
移動相関の計算を実行しましょう。
グラフ描画ライブラリ用いた処理ということもあり、処理時間がかなりかかります。
自環境では約7分かかりました。
### 移動相関関数の計算の実行 7分程度
result = calc_moving_xcorr(df_cor)
# データフレームの表示
display(result)【実行結果】
行：ラグ（秒）、列：時間（分）、値：移動相関です。
移動相関をヒートマップで描画します。
テキストの図5.2に相当します。
なお、（原因は不明ですが）テキストの描画結果とかなり異なっています。
移動相関の計算に誤りがあるかもしれません・・・。
### 身体動作の移動ラグ相関の描画 ※図5.2に相当 ★テキストの図と異なっている
# 描画領域
fig, ax = plt.subplots(figsize=(12, 5))
# ヒートマップの描画
sns.heatmap(result, cmap='jet', ax=ax, xticklabels=1500, yticklabels=50)
# 修飾
ax.set(xlabel='Time(min)', ylabel='Lag(sec)');【実行結果】
テキストによると、Lag（ラグ）が正のエリアで高い相関（赤色に近い）がある場合、セラピストの動作がクライアントに先行していることを示しているそうです。
移動相関の計算結果を格納したデータフレームをcsvファイルに保存しましょう。
### resultの保存
file = r'ch05_xcorr.csv'
result.to_csv(file)
１次元データの動的時間伸縮
動的時間伸縮法（DTW）は、時間間隔の異なる時系列データの類似度を評価できる指標です。
テキストの図5.3やDTW関連のWebサイトで、通常の距離（ユークリッド距離など）とDTWによってアラインメントされた距離のグラフを見比べて、動的時間伸縮のイメージを掴んでおきましょう。
① データの読み込みと前処理
Rの標準パッケージに含まれる「BJsales」データセットについて、①Webサイトより読み込み、②statsmodelsを利用した読み込みにより、２つのデータを取得します。
ts1は120期の先行指標、ts2は150期の販売実績です。
期間の長さが相違します。
### データセットの読み込みと前処理
# BJsales.leadデータセットの読み込み：WEBサイトから直接取得
url = 'https://vincentarelbundock.github.io/Rdatasets/csv/datasets/'\
'BJsales.lead.csv'
ts1 = pd.read_csv(url)
ts1 = ts1.drop(columns='rownames').iloc[:120]
# valueの標準化
ss = StandardScaler()
ts1['value'] = ss.fit_transform(ts1['value'].values.reshape(-1, 1))
# BJsalesデータセットの読み込み：statsmodelsを利用してRデータセットを取得
# 第１引数：Item, 第２引数：Package
ts2 = sm.datasets.get_rdataset('BJsales', 'datasets').data
# valueの標準化
ss = StandardScaler()
ts2['value'] = ss.fit_transform(ts2['value'].values.reshape(-1, 1))
# データセットの表示
print('ts1.shape:', ts1.shape)
display(ts1.head())
print('ts2.shape:', ts2.shape)
display(ts2.head())【実行結果（一部）】
値（value）は標準化されています。
② DTWの計算：基本形
Rのdtwパッケージの利用感に近そうな、Pythonのdtw-pythonパッケージを利用します。
dtw() の引数に２つのデータセットを指定して実行します。
### アラインメントを調べる
alignment = dtw(ts1['value'].values, ts2['value'].values, keep_internals=True)【実行結果】なし
DTWにより調整（アラインメント）された２つの指標の関係をプロットしましょう。
dtwで得た変数 alignment に対してメソッド plot を適用します。
引数 type でグラフの種類を指定します。
図5.4に相当します。
### 時点間の対応関係の描画 ※図5.4に相当
alignment.plot(type='threeway')
plt.suptitle('Timesereis alignment');【実行結果】
左と下に２つのデータセットが配置されています。
大きなグラフは、左・下の２つのグラフの時点を結ぶ軌跡になっています。
DTW距離を表示します。
dtwで得た変数 alignment に対してメソッド distance を適用します。
### DTW距離の表示
print(f'DTW(global): {alignment.distance:5.2f}')【実行結果】
テキストのDTW距離 26.25 に近似しています（一致しない理由は不明です。。。）。
２つの時系列データに関するDTWの対応関係を描画しましょう。
plot() の引数 type='twoway' を指定します。
### 点ごとの比較の描画
alignment.plot(type='twoway')
plt.suptitle('Alignment rawdata');【実行結果】
黒線が先行指標 ts1、青線が販売実績 ts2 です。
対応する時点を点線で結んでいます。
テキスト図5.5に相当するグラフを描画します。
描画関数を定義します。
### 整列後のQueryとReferenceの描画関数の定義
def alignment_plot2(alignment, title):
plt.figure(figsize=(8, 3))
plt.plot(range(alignment.index1.size), ts1.iloc[alignment.index1]['value'],
lw=0.9, label='Query')
plt.plot(range(alignment.index1.size), ts2.iloc[alignment.index2]['value'],
color='red', ls='--', label='Reference')
plt.title(title)
plt.xlabel('Time')
plt.legend()
plt.show()描画します。
### 整列後のBJsalesとBJsales.leadの描画 ※図5.5に相当
alignment_plot2(alignment, 'Global Alignment')【実行結果】
青い線が先行指標 ts1、赤い点線が販売実績 ts2 です。
DTWの対応関係を用いて、２つの時系列データを「長さの等しい時系列データ」に変換しています。
dtwで得た変数 alignment に対してメソッド index1、index2 を適用して変換後のデータを取得しています。
③ DTW：探索範囲の窓長に制限を加える
窓長を２つの時系列データの長さの差（30）に設定する方法を検討します。
テキストの DTW距離 30.43 のケースです。
ひとまずコード化しましたが、適切な内容かどうかは分からない状態です。。。
### アラインメントの探索
alignment_window = dtw(ts1['value'].values, ts2['value'].values,
window_type='sakoechiba',
window_args={'window_size': abs(len(ts1)-len(ts2))},
keep_internals=True)【実行結果】なし
DTW距離を算出します。
### DTW距離の表示
print(f'DTW(windowed, global): {alignment_window.distance:5.2f}')【実行結果】
テキストの 30.43 と比べて、若干相違します。。。
２つの時系列データの対応グラフを描画しましょう。
### 点ごとの比較の描画
alignment_window.plot(type='twoway')
plt.suptitle('Windowed Alignment');【実行結果】
整列後のデータをプロットしましょう。
### 整列後のBJsalesとBJsales.leadの描画
alignment_plot2(alignment_window, 'Windowed Alignment')【実行結果】
④ DTW：ローカルな整列を行う
パターンが似ている部分だけを抜き出して揃える方法を検討します。
テキストの DTW距離 20.08 のケースです。
時間の長さが120になります。
ひとまずコード化しましたが、適切な内容かどうかは分からない状態です。。。
引数 open_end と open_begin をTrue にしました。
### アラインメントの探索
alignment_openend = dtw(ts1['value'].values, ts2['value'].values,
step_pattern='asymmetric',
open_end=True, open_begin=True,
keep_internals=True)【実行結果】なし
DTW距離を算出しましょう。
### DTW距離の表示
print(f'DTW(local): {alignment_openend.distance:5.2f}')【実行結果】
テキストの 20.08 と比べて、若干相違します。。。
２つの時系列データの対応グラフを描画しましょう。
### 点ごとの比較の描画
alignment_openend.plot(type='twoway')
plt.suptitle('Open.beginとOpen.end');【実行結果】
整列後のデータをプロットしましょう。
テキストの図5.6に相当します。
### 整列後のBJsalesとBJsales.leadの描画 ※図5.6に相当
alignment_plot2(alignment_openend, 'Local Alignment')【実行結果】
２次元データの動的時間伸縮
２次元平面の移動の軌跡データを動的時間伸縮法（DTW）で整列（アラインメント）します。
① データの読み込みと前処理
RのMARSSパッケージに含まれる「loggerhead」データセットを利用します。
アメリカ東海岸の複数のアカウミガメの移動を記録した時系列の位置データ（緯度・経度）です。
いったんローカルＰＣにダウンロードして、pandasのデータフレームに読み込みます。
### rdaファイルのダウンロード
# 設定：ダウンロードファイルの格納先パス＋ファイル名
filename = r'./loggerheadNoisy.rda'
# ダウンロードの実行
url = r'https://github.com/cran/MARSS/blob/master/data/' \
'loggerheadNoisy.rda?raw=true'
res = requests.get(url)
with open(filename, mode='wb') as f:
f.write(res.content)【実行結果】なし
## loggerheadNoisy Rdataの読み込み https://pypi.org/project/rdata/
# ファイルの読み込み・解析
parsed = rdata.parser.parse_file(filename)
# データの変換：辞書型（中のデータ形式はxarray）
converted = rdata.conversion.convert(parsed)
# データフレーム化
logger_df = pd.DataFrame(converted['loggerheadNoisy'])
display(logger_df)【実行結果】
turtle列は８頭のウミガメの名称です。
DTW計算では、「BigMama」と「MaryLee」の２頭の位置データを利用します。
その他の項目は、日付：月・日・年、位置：経度lon・緯度latです、
テキストの図5.7に相当する８頭のウミガメの移動軌跡を描画します。
seabornのlineplotを利用します。
### 描画 ※図5.7に相当
# 描画領域
fig, ax = plt.subplots(figsize=(6, 6))
# seabornの折れ線グラフでウミガメごとに軌跡（経度lon, 緯度lat）を描画
sns.lineplot(data=logger_df, x='lon', y='lat', hue='turtle',
lw=0.8, ax=ax)
# 修飾
ax.set(xlabel='経度: lon', ylabel='緯度: lat')
ax.legend(bbox_to_anchor=(1, 1), title='Turtle')
ax.grid(lw=0.3)
plt.show()【実行結果】
南西（左下）・北東（右上）の範囲にとどまるウミガメ、大きく移動するウミガメの様子が分かります。
② DTWの計算
引き続き dtw-python パッケージを利用します。
まず、BigMama と MaryLee のデータを抽出します。
### データセットの準備
ts1 = (logger_df[logger_df['turtle']=='BigMama'][['lon', 'lat']]
.reset_index(drop=True))
ts2 = (logger_df[logger_df['turtle']=='MaryLee'][['lon', 'lat']]
.reset_index(drop=True))【実行結果】なし
２頭の軌跡データに含まれる期間を見てみましょう。
### ２頭の軌跡データの期間
print('BigMama:', len(ts1))
print('MaryLee:', len(ts2))【実行結果】
４期のずれがあることがわかりました。
続いて、dtw()を利用して、２頭のDTWを計算します。
### 時系列点の対応関係を調べる
# 多変量時系列の距離はlocal-distance matrixを用いて計算される
alignment = dtw(ts1, ts2, keep_internals=True)【実行結果】なし
DTW距離を表示します。
### DTW距離の表示
print(f'DTW Distance: {alignment.distance:5.2f}')【実行結果】
テキストのRスクリプトで計算したDTW距離は 80.556 。
一致しました。
③ DTWの対応関係の描画
２頭のウミガメの軌跡データの対応関係を描画します。
テキストの図5.8上に相当します。
### ワープ曲線の描画 ※図5.8上に相当
alignment.plot(type='alignment')
plt.title('Warping curve plot')
plt.xlabel('BigMama')
plt.ylabel('MaryLee');【実行結果】
４期ずれていた２頭の軌跡データが整列：アラインメントされました。
④ DTWの対応関係に基づく２頭の各時点の描画
２頭のウミガメの整列済みデータを用いて、時点（数）と位置（ｘ軸，ｙ軸）を描画します。
テキストの図5.8下に相当します。
### 描画 ※図5.8下に相当
## データの前処理
# BigMama
b = ts1.iloc[alignment.index1]
b['row_num'] = range(1, len(alignment.index1)+1)
# MaryLee
m = ts2.iloc[alignment.index2]
m['row_num'] = range(1, len(alignment.index2)+1)
## 描画
plt.figure(figsize=(5, 5))
# BigMamaの時点row_numの描画
for (x, y, num) in zip(b['lon'], b['lat'], b['row_num']):
plt.plot(x, y)
plt.text(x, y, str(num), color='tab:red')
# MaryLeeの時点row_numの描画
for (x, y, num) in zip(m['lon'], m['lat'], m['row_num']):
plt.plot(x, y)
plt.text(x, y, str(num), color='tab:blue')
# 凡例
plt.text(-73, 38.5, '凡例 Turtle:', color='black')
plt.text(-73, 38, ' 1 :BigMama', color='tab:red')
plt.text(-73, 37.5, ' 1 :MaryLee', color='tab:blue')
# 修飾
plt.grid(lw=0.2)
plt.xlim(-82.5, -73.5)
plt.ylim(29.5, 40)
plt.xlabel('経度: lon')
plt.ylabel('緯度: lat')
plt.show()【実行結果】
南西（左下）から北東（右上）へ移動する時点（数）が似通っていることが分かります。
以上で終了です。
お疲れ様でした。
第５章 その２ へ続く
■次回の取り組みテーマ
周波数領域の類似度に関連する「パワースペクトル」「コヒーレンス」「ウェーブレットクロススペクトル解析」「位相差」「レイリー検定」を実践します。
次の記事
前の記事
目次
ブログの紹介
note で４つのシリーズ記事を書いています。
ぜひ覗いていってくださいね！
１．のんびり統計
統計検定２級の問題集を手がかりにして、確率・統計をざっくり掘り下げるブログです。
雑談感覚で大丈夫です。ぜひ覗いていってくださいね。
統計検定２級公式問題集CBT対応版に対応しています。
２．実験！たのしいベイズモデリングをPyMC Ver.5で
書籍「たのしいベイズモデリング」の心理学研究に用いられたベイズモデルを PyMC Ver.5で描いて分析します。
この書籍をはじめ、多くのベイズモデルは Ｒ言語＋Stanで書かれています。
PyMCの可能性を探り出し、手軽にベイズモデリングを実践できるように努めます。
身近なテーマ、イメージしやすいテーマですので、ぜひぜひPyMCで動かして、一緒に楽しみましょう！
３．Python機械学習プログラミング実践記
書籍「Python機械学習プログラミング PyTorch & scikit-learn編」を学んだときのさまざまな思いを記事にしました。
この書籍は、scikit-learn と PyTorch の教科書です。
よかったらぜひ、お試しくださいませ。
４．データサイエンスっぽいことを綴る
統計、データ分析、AI、機械学習、Python のコラムを不定期に綴っています。
「統計」「Python」「数学とPython」「R」のシリーズが生まれています。
最後までお読みいただきまして、ありがとうございました。
この記事が気に入ったらサポートをしてみませんか？