【Kaggle】tsfreshを使って多次元時系列データを特徴量エンジニアリングしてみた


# 0. はじめに

突然ですが、時系列データって扱いが難しいですよね。
しかも、変数が増えるとますます心が折れそうになると思います。
でも、「時系列データから特徴量さえ抽出してしまえば、あとは何とでもなる！」
って人も多いのではないかと思います。

今回はそのような方に向けて、**tsfresh** という多次元時系列データの特徴量エンジニアリングに役立ちそうなライブラリを紹介します。

以下の記事を参考にさせて頂きました。
- [時系列データから自動で特徴抽出するライブラリ tsfresh](https://qiita.com/yuko1658/items/871df86f99a9134cc9ef)
- [tsfreshで時系列データの統計的処理を簡単に](https://qiita.com/KI1208/items/e27f827bac8654c9adcc)

# 1. tsfreshのインストール

僕は pip 経由でインストールしました。pip を新しめにしておかないと pip から install できなかったので、pip を upgrade しといて下さい。


```terminal
pip install --upgrade pip
```
で pip を upgrade して、

```terminal
pip install tsfresh
```
で tsfresh を install します。

※この際に、pandas が対応していないなどの警告が出た場合は pandas のバージョンを、

```terminal
pip install pandas==0.21 
```
みたいに変更してください。

僕の使用しているバージョンはこちらです。
- pip: 20.0.2
- pandas: 0.21.0
- tsfresh: 0.14.1


# 2. 擬似的な時系列データを用意する

多次元の時系列データを見つけるのが面倒だったので、tsfresh からダウンロードできるデータセットを今回は擬似的に変形して使用します。
(既に自分のデータがある方は読み飛ばしてください。)

先に言っておくと、この擬似データでは手順は掴めますが、出てくる結果は全く面白くないので、自分のデータがある方はそれを使用することをオススメします。

[UEA & UCR Time Series Classification Repository](http://timeseriesclassification.com/dataset.php)
に面白そうな時系列データがたくさんありそうですね...

まず、データを読み込みます。


```python:In[1]

import pandas as pd
import numpy as np

from tsfresh.examples.har_dataset import download_har_dataset, load_har_dataset

download_har_dataset()
df = load_har_dataset()
print(df.shape)
df.head()
```
で確認すると、このデータは7352サンプルポイント、128変数(128次元)あることが分かります。


次に、100サンプルポイント、50変数だけ切り取ります。

```python:In[2]
df = df.iloc[0:100, 0:50]
print(df.shape)
df.head()
```

今回は、
「被験者が5人いて、体につけたセンサーから10変数分の時系列データを取得して、その被験者が子供か大人かを分類する。」
みたいなシチュエーションを想定してください。

今回の目的は特徴量エンジニアリングの流れを通してみることなので、
値の対応はめちゃくちゃです。もちろんこのデータで分類できるわけもありません。


```python:In[3]
# id:5個、変数10個

#10変数ごとに、各個体（被験者）に割り当てる。
df_s1 = df.iloc[:,0:10].copy()
df_s2 = df.iloc[:,10:20].copy()
df_s3 = df.iloc[:,20:30].copy()
df_s4 = df.iloc[:,30:40].copy()
df_s5 = df.iloc[:,40:50].copy()

#それぞれの個体idを値とする列を作成する。
df_s1['id'] = 'sub1'
df_s2['id'] = 'sub2'
df_s3['id'] = 'sub3'
df_s4['id'] = 'sub4'
df_s5['id'] = 'sub5'

#それぞれのカラムの変数名を書き直す。
columns = ['x1', 'x2', 'x3', 'x4', 'x5', 'x6', 'x7', 'x8', 'x9', 'x10', 'id']
df_s1.columns = columns
df_s2.columns = columns
df_s3.columns = columns
df_s4.columns = columns
df_s5.columns = columns

df_s1.head()
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/541751/87bd4219-280b-c9f3-ce8d-9723c1e8bbe7.png)

こんな感じのデータフレームが5個体ぶんある状態です。

※よくあるミスですが、データフレームを分割抽出する際に、.copy() を使わないと、もとの値ごと変更されてしまいます(参照渡し)。この記事が参考になるかと思います。
[pandasで=でコピーした場合とcopy()を使った場合の挙動の違い](https://qiita.com/analytics-hiro/items/e8fdf8652869817851cd)


# 3. 時系列データを対応する形のデータフレームに変形する



[公式ドキュメント](https://tsfresh.readthedocs.io/en/latest/text/data_formats.html)によると、この記事での主人公的な関数になる **extract_features()** には、引数として渡す際の形式が指定されています。
データ型は pandas の dataframe オブジェクト型なのですが、その形式が3種類あります。

- Flat DataFrame (平坦なデータフレーム)
- Stacked DataFrame　(データフレームを縦積みしたもの)
- Dictionary of flat DataFrames (平坦なデータフレームを辞書にしたもの)

以上がその3種類です。今回は、1番目のこの形式に整形します。

```
id	time	x	y
A	t1	x(A, t1)	y(A, t1)
A	t2	x(A, t2)	y(A, t2)
A	t3	x(A, t3)	y(A, t3)
B	t1	x(B, t1)	y(B, t1)
B	t2	x(B, t2)	y(B, t2)
B	t3	x(B, t3)	y(B, t3)
```


先ほどの続きで、

```python:In[4]
df = pd.concat([df_s1, df_s2, df_s3, df_s4, df_s5], axis=0)
print(df['id'].nunique())
df.head()
```
として、連結すると、

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/541751/ff4945b4-db63-44f1-916c-c6f7fc6b3382.png)

となります。
ユニークな id の数が5個なので、問題なく連結されています。

これで、extract_features() 関数に渡せる形になりました。

# 4. extract_features() で特徴量を抽出する

先ほどのデータフレームに対して、

```python:In[5]
from tsfresh import extract_features

df_features = extract_features(df, column_id='id')
df_features.head()
```
を適用すると、
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/541751/86e5438c-c54b-f073-ba3c-ca954e9a93f9.png)

といった感じで特徴量が計算されます。
この特徴量は一つの変数に対して754個あります。
かなり多いと思いますが、時系列データ特有の扱い辛さが解消されたと思えば有用だと思います。

各特徴量が何を意味しているかは、
ドキュメンテーションの [Overview on extracted features](https://tsfresh.readthedocs.io/en/latest/text/list_of_features.html) に説明があります。パラメーターを要する特徴量(統計量)は複数のパラメーターで計算してくれているみたいです。

# 5. select_features()で特徴量をフィルタリングする

ドキュメンテーションでは、上記のように特徴量を抽出(extract)したあとは、[特徴量をフィルタリング]
(https://tsfresh.readthedocs.io/en/latest/text/feature_filtering.html) することが推奨されています。
その特徴量フィルタリングのための関数が **select_features()** です。
この関数は、その特徴量をもとにした場合に統計的な有意差がありそうな特徴量のみになるように、統計的仮説検定を用いて特徴量を選択してくれます。

その前に、

```python:In[6]
from tsfresh.utilities.dataframe_functions import impute

df_features = impute(df_features)
```
をすることで、先ほど求めた特徴量においてNaNや無限大などの扱えない値を補完します。

また、select_features() は従属変数 y をもとに特徴量を絞るので、そのデータを擬似的に用意します。
この y は[ソースコード](https://tsfresh.readthedocs.io/en/latest/_modules/tsfresh/feature_selection/selection.html)にもあるように、pandas.Seriesまたはnumpy.arrayである必要があります。
今回は、

```python:In[7]
from tsfresh import select_features

X = df_features
y = [0, 0, 0, 1, 1]
y = np.array(y)
```

と言った感じで特徴量抽出した X のデータフレームと、適当に割り振った y の numpy 配列を用意します。

そして、

```python:In[8]
X_selected = select_features(X, y)
print(X_selected.shape)
X_selected
```
をすると、

で〜ん。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/541751/48cdca8b-e115-1fc3-c9b0-5e48a09954ff.png)

見事に何も出てきません。
そりゃそうです、データが適当なやつなので...
ちゃんとしたデータですると、うまく特徴量選択できるはずです。
(自分のデータで特徴量選択ができた方(できなかった方)、はぜひコメントに一言頂けると泣いて喜びます。)

# おまけ. 主成分分析(PCA)で特徴量を次元圧縮してみた

これで終わるのもなんなので、
今回は統計的仮説検定を用いた特徴量選択の代わりに、PCAで次元圧縮してみました。
pca というインスタンスを生成して、それに fit して、そのインスタンスを使って transform する感じですね。
※ n_components の数は、個体数(今回は5)と特徴量数(今回は7540)の小さい方以下である必要があります。

```python:In[9]
from sklearn.decomposition import PCA
pca = PCA(n_components=4)
pca.fit(df_features)
X_PCA = pca.transform(df_features)
```

これをデータフレームに変形して表示すると、

```python:In[10]
X_PCA = pd.DataFrame( X_PCA )
X_PCA.head()
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/541751/7ddd5bdf-40eb-b23f-6849-9f34783c38b7.png)


となります。

さらに、
[Python: scikit-learn で主成分分析 (PCA) してみる](https://blog.amedama.jp/entry/2017/04/02/130530)
を参考にして、寄与率・累積寄与率を求めてみると、

```python:In[11]
print('各次元の寄与率: {0}'.format(pca.explained_variance_ratio_))
print('累積寄与率: {0}'.format(sum(pca.explained_variance_ratio_)))
```

```text:out[11]
各次元の寄与率: [0.30121012 0.28833114 0.22187195 0.1885868 ]
累積寄与率: 0.9999999999999999
```
となります。
もしかすると、
特徴量選択の代わりに PCA で特徴量圧縮してもいいのかもしれないですね。


パワープレーではありますが、この一連の流れを使えば時系列ならではの煩わしさからは逃れられそうですね。

そこからは、
[【Kaggle】Baseline modelの構築、Pipeline処理](https://qiita.com/yuki_edy/items/0923a69855a8f0042070)
みたいにベースラインモデルを作ってみたりなどなど、いろいろな手が打てますね。

次は実際の多次元時系列データで試してみようと思います。