◇◆機械学習モデリング前編 Jane Street Market Predictionで市場予測をしよう！◆◇

# モチベーション

本職でずっと投資の研究をしているのですが、アイデアや戦略は結構思いついたりします。しかし、経験者であれば分かる通り、上がるか下がるかを予測することは非常に難易度が高いです。

また、アイデア自体は思いついても正直枯渇します。スポンジのように知識を吸収しては、絞り切ってを繰り返すと思います。

なので、ある程度、市場予測の一連の作業を楽に出来ないか考えた所、学習コストはかかるものの、モデリングを鍛えるという結論に到達しました。

本記事では、Kaggleのnotebookを題材にモデリングについて解説していきたいと思います。

ただ、もちろん時系列データの扱いは非常に繊細なので、検証段階で将来の情報リーク等によって良い結果を出すこと等もあり、そういった時系列データの取り扱いの細かい部分は [◆アルゴリズム取引　～投資の定量分析の犯しやすいミスとその解決策～◆](https://qiita.com/kanawoinvest/items/7d830b1ce2ba3dcccb8f) で解説してますので、そちらをご参照くださると助かります。

# 方針

Kaggleの[Jane Street Market Prediction](https://www.kaggle.com/c/jane-street-market-prediction)からnotebookを参考にしながら分析手法やモデリング手法について学んでいきたいと思います。

本記事で採用するnotebookは[Jane_street_Extensive_EDA](https://www.kaggle.com/muhammadmelsherbini/jane-street-extensive-eda)で、こちらでは本コンペのデータの理解を深めるためにデータを細かく見ます。

ですので、**モデリング編では、データ探索をして、データへの理解度を深めてから、様々な機械学習モデルを用いて市場予測していく方針です。**

# データの説明

*参考:https://www.kaggle.com/c/jane-street-market-prediction/data*

このデータセットには、実際の株式市場データを表す匿名化された一連の特徴量 `feature_ {0 ... 129}` が含まれています。データセットの各行は取引機会を表しており、`action` を予測します。1は取引を行い、0は取引を渡します。各取引には関連する`weight`と`resp`があり、これらを合わせて取引の収益を表します。`date` 列は取引の日を表す整数であり、`ts_id` は時間の順序を表します。匿名化された特徴量に加えて、**features.csv**では特徴量に関するメタデータが提供されます。

訓練データ**train.csv**には、`resp`値と、さまざまな期間のリターンを表す他のいくつかの`resp_ {1,2,3,4}`値が提供されます。これらの変数はテストセットに含まれていません。`weight = 0`の取引は、完全を期すために意図的にデータセットに含まれていますが、そのような取引はスコア評価に貢献しません。

これは、時系列APIに依存して、モデルが時間的に前向きにならないようにするコード競争です。 APIを使用するには、評価ページの指示に従ってください。ノートブックを送信すると、目に見えないテストで再実行されます。

- コンペのモデルトレーニングフェーズ中、この目に見えないテストセットは、約100万行の履歴データで構成されます。
- ライブ予測フェーズでは、テストセットは定期的に更新されるライブ市場データを使用します。

コンペの第2（予測）フェーズでは、ノートブックの時間制限は、テストセットで提示された取引の数に比例することに注意してください。詳細については、コード要件を参照してください。

# 成績評価基準

この大会は、効用スコア(*utility score*)で評価されます。テストセットの各行は、アクション値を予測する取引機会を表します。1は取引を行い、0は取引を見送ります。各トレード`j`には、リターンを表す重み(`weight`)と各アセットのリターン(`resp`)が関連付けられています。

各`date`に対して、

```math
\displaystyle p_i = \sum_{j}(weight_{ij} * resp_{ij} * action_{ij})
```

```math
\displaystyle t = \frac{\sum p_i}{\sqrt{\sum p_i^2}} * \sqrt{\frac{250}{|i|}}
```

ただし、`i`はテスト期間における独自な日付の個数です。効用スコアは次のように定義します。

```math
u = \min( \max(t, 0), 6)\sum p_i
```

**解説**: ある銘柄`j`を考えますと株価からリターンを計算出来ます。ここで、株価のリターンをそれぞれ$resp_1, resp_2, ..., resp_i$となっています。また、重みは本コンペで定義された重みですが、ポートフォリオ理論的には銘柄`i`に資産をどの程度分配したかを表します。なので、本コンペが`weight = 0`と書いてある銘柄に関してはそもそも評価に貢献しません。

なので、`weight * resp`で個別銘柄の戦略による損益で、`action`はそれを取引するかどうかです。$p_i$では全銘柄の投資戦略の損益を表していると言えます。また、$t$に関しては、投資戦略のシャープレシオを表していまして、$\sqrt(250/i)$は投資戦略によって選んだ銘柄の個数を調整するための掛け目です。

最後に$u$ですが、こちらが評価関数になります。意味はシャープレシオ*リターンで、シャープレシオには、上限値と下限値を設けております。

前置きが長くなりましたが、さっそくデータについて見ていきましょう！

# notebook解説

## インポート

```python
import os
import gc
import sys

import numpy as np
import pandas as pd

import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
import seaborn as sns
import plotly.express as px
import plotly.graph_objects as go
%matplotlib inline
sns.set_style('darkgrid')
pd.set_option('display.max_columns', 140)
# from sklearn.preprocessing import StandardScaler as scale
# from sklearn.decomposition import PCA
# from sklearn.neural_network import MLPClassifier as mlp
# from sklearn.linear_model import LogisticRegression as log
# from sklearn.model_selection import train_test_split as split 
# from sklearn.model_selection import GridSearchCV as Grid
# from sklearn.naive_bayes import GaussianNB as GNB
# from sklearn.svm import SVC
# from sklearn.linear_model import LinearRegression as linear
# import xgboost as xgb
# from sklearn.decomposition import IncrementalPCA as ipca
# from sklearn.metrics import (roc_auc_score, precision_score, recall_score, f1_score,
#                              confusion_matrix, accuracy_score, roc_curve, auc)
```

コメントアウトしている行に関しましては、基本的にモデリングに使う行ですが、今回はデータセットの解説から入るため、こちらは一旦不要です。

## データの読み込み

```python
df = pd.read_csv('../input/jane-street-market-prediction/train.csv')
```

```python
# 元データのサイズを保存
org_len = len(df)
```

データの大きさを確認

```python
df.info()
```

```
<class 'pandas.core.frame.DataFrame'>
RangeIndex: 2390491 entries, 0 to 2390490
Columns: 138 entries, date to ts_id
dtypes: float64(135), int64(3)
memory usage: 2.5 GB
```

今回の訓練用のデータセットは138列、2,390,491行あります。また、メモリは2.5GB使うため非常に重いです。このままデータを用いるとコンピュータに負担がかかるのでメモリを節約しましょう。

まずはいくつかの変数の型を変更します。

```python
df.feature_0 = df.feature_0.astype(np.int8)
df.date= df.date.astype(np.int16)
df.ts_id = df.ts_id.astype(np.int32)
```

次のコードは、データフレームをループしながら、-0.0001から0.0001の間のデータの量が非常に少ない場合（0.1％の非常に控えめなスレッショルド）にのみfloat64列をfloat32に変更して、小さな値の列の精度を損なうことを回避します。

```python
for i in df:
    if df[i].dtype == np.float64:
        if (((df[i] < .0001) & (df[i] > -.0001)).mean()) > .001:
            print(i)
```

```
weight
resp_1
resp_2
resp_3
resp_4
resp
```

ここで、これらのカラムに対してメモリ削減のために、型を変更、および、ガベージコレクション(`gc`)でメモリ管理をします。

```
for i in df:
    if df[i].dtype == np.float64:
        if (((df[i] < .0001) & (df[i] > -.0001)).mean()) < .001:
            df[i] = df[i].astype(np.float32)
            gc.collect();
```

では、もう一度、このデータセットのメモリを見てみましょう。

```python
df.info()
```

```
<class 'pandas.core.frame.DataFrame'>
RangeIndex: 2390491 entries, 0 to 2390490
Columns: 138 entries, date to ts_id
dtypes: float32(129), float64(6), int16(1), int32(1), int8(1)
memory usage: 1.3 GB
```

なんと、2.5GBのメモリを消費するデータセットを1.3GBまでに落とすことができました。これで計算も速くなりますし、メモリ不足の問題も多少は改善されると思います。

次に、データセットは時系列データなので、時系列にちゃんとなっているように並び替えましょう。

```python
df.sort_values(by= ['date','ts_id'],inplace=True)
```

次に`action`を定義します。これは、銘柄を買うか売るかでしたね。

```python
df['action'] = np.where(df['resp'] > 0,1,0)
df.action = df.action.astype('category')
```


ここで`np.where`が出てきましたが、これは`np.where(condition, x, y)`で第一引数に条件文を入れて、それが真なら第二引数の値、偽であれば第三引数の値になります。

これはカテゴリカル変数として定義しておきます。つまり、1か0には数値的な意味はなく、買ったか買っていないかを分類するための数字でしかないということです。

# EDA(探索的データ解析)

EDAはデータセットを調査して、データの中身について理解を深めることです。これによって、データを用いて予測する際に、どの変数を用いるか、もしくはどのような変数を生成すれば良いのか等について知見を深めることが出来ます。名前はかっこいいですね。枝と読みたい方もいらっしゃるでしょうけど、残念ながら巷ではイーディーエーと読みます。

## 一般的なデータ分析

### `Resp`データ分析

```python
fig = plt.figure(figsize=(16,6))
ax = plt.subplot(1,1,1)
df.groupby('date')[['resp_1', 'resp_2', 'resp_3', 'resp_4', 'resp']].sum().cumsum().plot(ax=ax)
plt.title('Cumulative Sum of Different RESP\'s',fontsize=18)
plt.xlabel('Date',fontsize=14)
plt.axvspan(0,92,linestyle=':',linewidth=2,label='first 92 days',color='darkorange',alpha=.2)
plt.legend(fontsize=12,ncol=3,loc=2);
```

![__results___22_0.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/729216/5d429d94-fd3b-34f7-2a0c-2dcfd5db29ae.png)

>最初の92日間でより多くの利益があったことに気付くことができます。これは、多く分析者に、この時点より前に観測を取り除くことをお勧めするという洞察につながります。
また、`resp_4`の累積合計が最大であるのに対し、`resp_1`の累積合計は最小であることがわかります。

では、`Resp`の日毎の平均をプロットしてみましょう。

```python
fig = px.line(df.groupby('date')[['resp_1', 'resp_2', 'resp_3', 'resp_4','resp']].mean(),
              x= df.groupby('date')[['resp_1', 'resp_2', 'resp_3', 'resp_4','resp']].mean().index,
              y= ['resp_1', 'resp_2', 'resp_3', 'resp_4','resp'],
              title= '\naverage Resp per day')
fig.layout.xaxis.title = 'Day' 
fig.layout.yaxis.title = 'Avg Resp'
fig.show()
```

![newplot.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/729216/da29f952-95f7-f004-070f-95f28c5a4fe4.png)


```python
fig,((ax11,ax12,ax13),(ax21,ax22,ax23),(ax31,ax32,ax33),(ax41,ax42,ax43),(ax51,ax52,ax53)) = plt.subplots(5,3,figsize=(18,24))
plt.subplots_adjust(hspace=0.35)
ax11.hist(df.resp,bins=150, color='darkblue',alpha=.6)
ax11.axvline(df.resp.mean()+df.resp.std(),color='darkorange',linestyle=':',linewidth=2)
ax11.axvline(df.resp.mean()-df.resp.std(),color='darkorange',linestyle=':',linewidth=2)
df.resp.plot.hist(bins= 150,ax=ax12,color='darkblue',alpha=.6)
ax12.axvline(df.resp.mean()+df.resp.std(),color='darkorange',linestyle=':',linewidth=2)
ax12.axvline(df.resp.mean()-df.resp.std(),color='darkorange',linestyle=':',linewidth=2)
ax12.set_xlim(-.08,.08)
ax13.hist(df.resp,bins=150, color='darkblue',alpha=.6)
ax13.set_yscale('log')
skew= round(df.resp.skew(),4)
kurt= round(df.resp.kurtosis())
std1= round((((df.resp.mean()-df.resp.std()) < df.resp ) & (df.resp < (df.resp.mean()+df.resp.std()))).mean()*100,2)
props = dict(boxstyle='round', facecolor='white', alpha=0.5)
ax11.text(.02,.96,'μ = {}\nstd = {}\nskewness = {}\nkurtosis = {}\n% values in 1 std = {}%'.format(round(df.resp.mean(),4),round(df.resp.std(),4),skew,kurt,std1),
         transform=ax11.transAxes, verticalalignment='top',bbox=props,fontsize=10)
ax11.set_title('Resp Hist Normal scale',fontsize=14)
ax12.set_title('Resp Hist normal scale zoomed',fontsize=14)
ax13.set_title('Resp Hist with freq on a log scale',fontsize=14);
ax11.set_xlabel('')
ax11.set_ylabel('')
ax12.set_xlabel('')
ax12.set_ylabel('')
ax13.set_xlabel('')
ax13.set_ylabel('')
ax21.hist(df.resp_1,bins=150,color='darkblue',alpha=.6)
ax21.axvline(df.resp_1.mean()+df.resp_1.std(),color='darkorange',linestyle=':',linewidth=2)
ax21.axvline(df.resp_1.mean()-df.resp_1.std(),color='darkorange',linestyle=':',linewidth=2)
df.resp_1.plot.hist(bins= 150,ax=ax22,color='darkblue',alpha=.6)
ax22.axvline(df.resp_1.mean()+df.resp_1.std(),color='darkorange',linestyle=':',linewidth=2)
ax22.axvline(df.resp_1.mean()-df.resp_1.std(),color='darkorange',linestyle=':',linewidth=2)
ax22.set_xlim(-.08,.08)
ax23.hist(df.resp_1,bins=150,color='darkblue',alpha=.6)
ax23.set_yscale('log')
skew= round(df.resp_1.skew(),4)
kurt= round(df.resp_1.kurtosis())
std1= round((((df.resp_1.mean()-df.resp_1.std()) < df.resp_1 ) & (df.resp_1 < (df.resp_1.mean()+df.resp_1.std()))).mean()*100,2)
ax21.text(.02,.96,'μ = {}\nstd = {}\nskewness = {}\nkurtosis = {}\n% values in 1 std = {}%'.format(round(df.resp_1.mean(),4),round(df.resp_1.std(),4),skew,kurt,std1),
         transform=ax21.transAxes, verticalalignment='top',bbox=props,fontsize=10)
ax21.set_title('Resp_1 Hist Normal scale',fontsize=14)
ax22.set_title('Resp_1 Hist normal scale zoomed',fontsize=14)
ax23.set_title('Resp_1 Hist with freq on a log scale',fontsize=14);
ax21.set_xlabel('')
ax21.set_ylabel('')
ax22.set_xlabel('')
ax22.set_ylabel('')
ax23.set_xlabel('')
ax23.set_ylabel('')
ax31.hist(df.resp_2,bins=150,color='darkblue',alpha=.6)
ax31.axvline(df.resp_2.mean()+df.resp_2.std(),color='darkorange',linestyle=':',linewidth=2)
ax31.axvline(df.resp_2.mean()-df.resp_2.std(),color='darkorange',linestyle=':',linewidth=2)
df.resp_2.plot.hist(bins= 150,ax=ax32,color='darkblue',alpha=.6)
ax32.axvline(df.resp_2.mean()+df.resp_2.std(),color='darkorange',linestyle=':',linewidth=2)
ax32.axvline(df.resp_2.mean()-df.resp_2.std(),color='darkorange',linestyle=':',linewidth=2)
ax32.set_xlim(-.08,.08)
ax33.hist(df.resp_2,bins=150, color='darkblue',alpha=.6)
ax33.set_yscale('log')
skew= round(df.resp_2.skew(),4)
kurt= round(df.resp_2.kurtosis())
std1= round((((df.resp_2.mean()-df.resp_2.std()) < df.resp_2 ) & (df.resp_2 < (df.resp_2.mean()+df.resp_2.std()))).mean()*100,2)
ax31.text(.02,.96,'μ = {}\nstd = {}\nskewness = {}\nkurtosis = {}\n% values in 1 std = {}%'.format(round(df.resp_2.mean(),4),round(df.resp_2.std(),4),skew,kurt,std1),
         transform=ax31.transAxes, verticalalignment='top',bbox=props,fontsize=10)
ax31.set_title('Resp_2 Hist Normal scale',fontsize=14)
ax32.set_title('Resp_2 Hist normal scale zoomed',fontsize=14)
ax33.set_title('Resp_2 Hist with freq on a log scale',fontsize=14);
ax31.set_xlabel('')
ax31.set_ylabel('')
ax32.set_xlabel('')
ax32.set_ylabel('')
ax33.set_xlabel('')
ax33.set_ylabel('')
ax41.hist(df.resp_3, color='darkblue',alpha=.6,bins=150)
ax41.axvline(df.resp_3.mean()+df.resp_3.std(),color='darkorange',linestyle=':',linewidth=2)
ax41.axvline(df.resp_3.mean()-df.resp_3.std(),color='darkorange',linestyle=':',linewidth=2)
df.resp_3.plot.hist(bins=150, color='darkblue',alpha=.6,ax=ax42)
ax42.axvline(df.resp_3.mean()+df.resp_3.std(),color='darkorange',linestyle=':',linewidth=2)
ax42.axvline(df.resp_3.mean()-df.resp_3.std(),color='darkorange',linestyle=':',linewidth=2)
ax42.set_xlim(-.08,.08)
ax43.hist(df.resp_3, color='darkblue',alpha=.6,bins=150)
ax43.set_yscale('log')
skew= round(df.resp_3.skew(),4)
kurt= round(df.resp_3.kurtosis())
std1= round((((df.resp_3.mean()-df.resp_3.std()) < df.resp_3 ) & (df.resp_3 < (df.resp_3.mean()+df.resp_3.std()))).mean()*100,2)
ax41.text(.02,.96,'μ = {}\nstd = {}\nskewness = {}\nkurtosis = {}\n% values in 1 std = {}%'.format(round(df.resp_3.mean(),4),round(df.resp_3.std(),4),skew,kurt,std1),
         transform=ax41.transAxes, verticalalignment='top',bbox=props,fontsize=10)
ax41.set_title('Resp_3 Hist Normal scale',fontsize=14)
ax42.set_title('Resp_3 Hist normal scale zoomed',fontsize=14)
ax43.set_title('Resp_3 Hist with freq on a log scale',fontsize=14);
ax41.set_xlabel('')
ax41.set_ylabel('')
ax42.set_xlabel('')
ax42.set_ylabel('')
ax43.set_xlabel('')
ax43.set_ylabel('')
ax51.hist(df.resp_4,bins=150, color='darkblue',alpha=.6)
ax51.axvline(df.resp_4.mean()+df.resp_4.std(),color='darkorange',linestyle=':',linewidth=2)
ax51.axvline(df.resp_4.mean()-df.resp_4.std(),color='darkorange',linestyle=':',linewidth=2)
df.resp_4.plot.hist(bins= 150,color='darkblue',alpha=.6,ax=ax52)
ax52.axvline(df.resp_4.mean()+df.resp_4.std(),color='darkorange',linestyle=':',linewidth=2)
ax52.axvline(df.resp_4.mean()-df.resp_4.std(),color='darkorange',linestyle=':',linewidth=2)
ax52.set_xlim(-.08,.08)
ax53.hist(df.resp_4,bins=150,color='darkblue',alpha=.6)
ax53.set_yscale('log')
skew= round(df.resp_4.skew(),4)
kurt= round(df.resp_4.kurtosis())
std1= round((((df.resp_4.mean()-df.resp_4.std()) < df.resp_4 ) & (df.resp_4 < (df.resp_4.mean()+df.resp_4.std()))).mean()*100,2)
ax51.text(.02,.96,'μ = {}\nstd = {}\nskewness = {}\nkurtosis = {}\n% values in 1 std = {}%'.format(round(df.resp_4.mean(),4),round(df.resp_4.std(),4),skew,kurt,std1),
         transform=ax51.transAxes, verticalalignment='top',bbox=props,fontsize=10)
ax51.set_title('Resp_4 Hist Normal scale',fontsize=14)
ax52.set_title('Resp_4 Hist normal scale zoomed',fontsize=14)
ax53.set_title('Resp_4 Hist with freq on a log scale',fontsize=14)
ax51.set_xlabel('')
ax51.set_ylabel('')
ax52.set_xlabel('')
ax52.set_ylabel('')
ax53.set_xlabel('')
ax53.set_ylabel('')
fig.suptitle('RESPs Historgrams on Different Scales',fontsize=18,y=.92);
```

![__results___26_0.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/729216/b24fa4ca-2a4c-8ea4-7d6f-f2566be6de99.png)


次にペアプロットです。

```python
sns.pairplot(df[['resp_1', 'resp_2', 'resp_3', 'resp_4', 'resp']],corner=True);
```

![__results___27_0.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/729216/d7e94aa0-84cc-937d-4d55-94520092e0f3.png)


>この図に示されている関係と、この前の図に示されている標準偏差と分布から、`Resp`は`Resp_4`との関連性が高く、`Resp_1`と`Resp_2`は相互に関連性が高いことがわかります。より長い期間の投資はより多くのリターンとより高いリスクに繋がります。

**投資期間に関する基礎知識**


投資期間は、特定の目標のために投資を行うことが期待される期間です。投資は一般に、株式（リスクが高い）と債券（リスクが低い）の2つの主要なカテゴリに分類されます。投資期間(以下、タイム・ホライズン)が長ければ長いほど、投資家はより積極的またはリスクの高いポートフォリオを構築できます。タイム・ホライズンが短ければ短いほど、投資家はより保守的またはリスクが少なくなります。

### 日付(`date`)

今回は日付がちょうど500日分あるため、これは大体２年分のデータであると判断して良いです。理由は、株式市場の営業日が年間大体250日であるからです。$t$の定義式の掛け目も一年間=250営業日と定義しているように思われます。

では、日別に`resp`がある銘柄数について見てみましょう。株式市場では、もちろん全銘柄に対してリターンは存在しますが、本コンペでは、ある程度投資対象を絞っているため、全銘柄について調査しなくて済みます。

```python
fig = px.area(data_frame= df.groupby('date')[['resp']].count(),title='Number of operation per day')
fig.update_traces( showlegend = False)
fig.layout.xaxis.title = 'Day' 
fig.layout.yaxis.title = 'Number of operations'
fig.show()
```

![newplot (1).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/729216/705d7750-46e5-fb0c-dd20-838016181bab.png)

日別合計リターン

```python
fig = px.area(data_frame= df.groupby('date')[['resp']].sum(),title='Resp sum of operation per day')
fig.update_traces( showlegend = False)
fig.layout.xaxis.title = 'Day' 
fig.layout.yaxis.title = 'Resp sum'
fig.show()
```

![newplot (2).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/729216/60b19b8a-68b1-bb6e-a23d-864eab394868.png)

>`resp`の分散が大きいことが伺えます。

それぞれの平均に対して20営業日（1ヶ月=約20営業日を示唆）の移動標準偏差を作成します。

```python
date_df = df.groupby('date')[['resp']].mean()
std20 = []
for i in range(len(date_df)):
    if i <20:
        std20.append(np.nan)
    else:
        moving_std = date_df['resp'][i-20:i].std()
        std20.append(moving_std)
date_df['moving_std'] = std20
date_df.tail(2)
```

||resp|moving_std|
|---|---|---|
|date|||	
|498|0.001554|0.001962
|499|0.000299|0.001969


```python
fig = px.line(data_frame=date_df,y=['resp','moving_std'],title='Average Resp & 20 day moving standard deviation')
fig.layout.xaxis.title = 'Day' 
fig.layout.yaxis.title = 'Avg Resp'
fig.show()
```
![newplot (3).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/729216/3bda0cfe-7dc9-e36e-81b0-b9d76a6f574e.png)

ここで、毎日の各`resp`の標準偏差を確認します。

```python
fig, (ax1,ax2) = plt.subplots(2,1,figsize=(14,12))
df.groupby('date')[['resp_1', 'resp_2', 'resp_3', 'resp_4']].std().plot(ax=ax1,color=['steelblue','darkorange','red','green'],alpha=.8)
df.groupby('date')[['resp_1', 'resp_2', 'resp_3', 'resp_4']].std().plot.kde(ax=ax2)
fig.suptitle('Resp\'s Std',fontsize=18,y=.96)
ax2.set_xlabel('')
ax1.set_xlabel('')
ax2.set_title('kde of each resp std', fontsize=14)
ax1.set_title('std of Resp\'s for each trading day',fontsize=14);
```

![__results___41_0.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/729216/cc4920ab-1d90-98e8-8cbf-4b27e7f4dc37.png)

>前に述べたように、`resp`の標準偏差は、長期的な投資とともに増加するようです。
>また、多くのKagglerが80日後に何らかの取引モデルの調整が行われた可能性があると述べたため、最初の100日間で偏差がやや高かったことにも気付くことができます。

### 重み(`weight`)

```python
fig = plt.figure(figsize=(18,7))
grid =  gridspec.GridSpec(2,3,figure=fig,hspace=.3,wspace=.2)
ax1 = fig.add_subplot(grid[0, 0])
ax2 = fig.add_subplot(grid[0, 1])
ax3 = fig.add_subplot(grid[1, 0])
ax4 = fig.add_subplot(grid[1, 1])
ax5 = fig.add_subplot(grid[:, 2])
sns.boxplot(x = df.weight,width=.5,ax=ax1)
ax2.hist(df.weight, color='#404788ff',alpha=.6, bins= list([-.05] + list(10**np.arange(-2,2.24,.05))))
ax2.set_xscale('symlog')
ax2.set_xlim(-.05,227)
sns.boxplot(x = df.weight[df.weight != 0],width=.5,ax=ax3)
ax1.set_title('Weights including zero weights',fontsize=14)
ax3.set_title('Weights not including zero weights',fontsize=14)
ax2.set_title('Weights including zero weights (log)',fontsize=14)
ax4.set_title('Weights not including zero weights (log)',fontsize=14)
props = dict(boxstyle='round', facecolor='white', alpha=0.4)
ax1.text(.2,.9,'μ = {}    std = {}\nmin = {}    max = {}'.format(round(df.weight.mean(),3),round(df.weight.std(),3),round(df.weight.min(),3),round(df.weight.max(),3)),
         transform=ax1.transAxes, verticalalignment='top',bbox=props,fontsize=12)
ax3.text(.2,.9,'μ = {}        std = {}\nmin = {}    max = {}'.format(round(df.weight[df.weight != 0].mean(),3),round(df.weight[df.weight != 0].std(),3),
                                                              round(df.weight[df.weight != 0].min(),3),round(df.weight[df.weight != 0].max(),3)),
         transform=ax3.transAxes, verticalalignment='top',bbox=props,fontsize=12)
ax4.hist(df.weight[df.weight !=0],color='#404788ff',alpha=.6,bins=10**np.arange(-2.16,2.24,.05))
ax4.set_xscale('log')
ax4.set_xticks((.01,.03,.1,.3,1,3,10,30,100))
ax4.set_xticklabels((.01,.03,.1,.3,1,3,10,30,100))
ax5.pie(((df.weight==0).mean(),(1-(df.weight==0).mean())),startangle=300,wedgeprops=dict(width=0.5),
        labels=('Zeros\n{}%'.format(round((df.weight==0).mean()*100,2)),'Nonzeros\n{}%'.format(round((1-(df.weight==0).mean())*100,2))),
        textprops={'fontsize': 12},colors=['#404788ff','#55c667ff'])
ax5.set_title('Zeros vs non-zero weights',fontsize=14)
ax1.set_xlabel('')
ax2.set_xlabel('')
ax3.set_xlabel('')
ax2.set_ylabel('')
ax5.set_ylabel('')
ax4.set_xlabel('');
```

![__results___44_0.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/729216/aaa3e807-8d21-9454-1b5c-a02af90e25b9.png)

```python
fig = plt.figure(figsize=(15,10))
fig.suptitle('Nonzero weights histogram in different scales',fontsize=18)
ax1 = plt.subplot(3,1,1)
ax1.hist(df.weight[df.weight !=0],color='darkblue',alpha=.7, bins=10**np.arange(-2.16,2.23,.05))
plt.xscale('log')
plt.xticks((.01,.03,.1,.3,1,3,10,30,100),(.01,.03,.1,.3,1,3,10,30,100))
ax2 = plt.subplot(3,1,2)
sns.distplot(df.weight[df.weight != 0], color='darkblue', bins=400, ax=ax2) 
ax3 = plt.subplot(3,1,3)
ax3.hist(df.weight[(df.weight !=0) & (df.weight < 3.197 )],color='darkblue',alpha=.7, bins=200)
ax3.set_xlim(0,3.3)
ax2.set_xlabel('') 
ax1.set_title('All values (log-scale)',fontsize=14)
ax2.set_title('kde of the distribution',fontsize=14)
ax3.set_title('75% of the Values',fontsize=14)
plt.subplots_adjust(hspace=.4);
```

![__results___45_0.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/729216/409d28d4-6999-0638-f346-a75026e3dd2d.png)

> 重みの値が0に集中していて、かなり分散が激しいですね。特に右端の外れ値に関してより詳しく見ていきましょう。

```python
fig, (ax1,ax2) = plt.subplots(2,1,figsize=(16,8))
fig.suptitle('Weight outliers',fontsize=18)
sns.boxplot(df.weight,width=.5, ax=ax1)
ax1.axvline(np.percentile(df.weight,95), color= 'green',label='95.0%',linestyle=':',linewidth=3)
ax1.axvline(np.percentile(df.weight,99), color= 'darkblue',label='99.0%',linestyle=':',linewidth=3)
ax1.axvline(np.percentile(df.weight,99.9), color= 'darkorange',label='99.9%',linestyle=':',linewidth=3)
ax1.axvline(np.percentile(df.weight,99.99), color= 'magenta',label='99.99%',linestyle=':',linewidth=3)
ax1.legend(fontsize=13)
sns.boxplot(df.weight[df.weight !=0],width=.5, ax=ax2)
ax2.axvline(np.percentile(df.weight[df.weight !=0],95), color= 'green',label='95.0%',linestyle=':',linewidth=3)
ax2.axvline(np.percentile(df.weight[df.weight !=0],99), color= 'darkblue',label='99.0%',linestyle=':',linewidth=3)
ax2.axvline(np.percentile(df.weight[df.weight !=0],99.9), color= 'darkorange',label='99.9%',linestyle=':',linewidth=3)
ax2.axvline(np.percentile(df.weight[df.weight !=0],99.99), color= 'magenta',label='99.99%',linestyle=':',linewidth=3)
ax2.legend(fontsize=13)
ax1.set_title('All weights', fontsize= 14)
ax2.set_title('Non-zero weights', fontsize= 14)
ax1.set_xlabel('')
ax2.set_xlabel('');
```

![__results___47_0.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/729216/da3e4ea3-1cbb-82c0-2e7b-459013b5523d.png)

```python
sns.scatterplot(data=df, x='resp',y='weight', color= 'blue', alpha=.3)
plt.title('Resp vs Weight\ncorrelation={}'.format(round(df.weight.corr(df.resp),4)));
```

![__results___48_0.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/729216/6923df20-df5e-6380-5852-34ab0e72eea4.png)

>重みは`Resp`と線形に相関していないことがわかりますが、重みが大きいほど`Resp`値が低いことにのみ関連していることは明らかです。

## 特徴量データ分析

### 特徴量csvを読み込む

```python
df_f =  pd.read_csv('../input/jane-street-market-prediction/features.csv')
df_f.head(5)
```
||feature|tag_0|tag_1|tag_2|tag_3|tag_4|tag_5|tag_6|tag_7|tag_8|tag_9|tag_10|tag_11|tag_12|tag_13|tag_14|tag_15|tag_16|tag_17|tag_18|tag_19|tag_20|tag_21|tag_22|tag_23|tag_24|tag_25|tag_26|tag_27|tag_28
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
0|feature_0|False|False|False|False|False|False|False|False|False|False|False|False|False|False|False|False|False|False|False|False|False|False|False|False|False|False|False|False|False
1|feature_1|False|False|False|False|False|False|True|True|False|False|False|False|False|False|False|False|False|False|False|False|False|False|False|False|False|False|False|False|False
2|feature_2|False|False|False|False|False|False|True|True|False|True|False|False|False|False|False|False|False|False|False|False|False|False|False|False|False|False|False|False|False
3|feature_3|False|False|False|False|False|False|True|False|True|False|False|False|False|False|False|False|False|False|False|False|False|False|False|False|False|False|False|False|False
4|feature_4|False|False|False|False|False|False|True|False|True|True|False|False|False|False|False|False|False|False|False|False|False|False|False|False|False|False|False|False|False

```python
fig = px.bar(df_f.set_index('feature').T.sum(), title='Number of tags for each feature')
fig.layout.xaxis.tickangle = 300
fig.update_traces( showlegend = False)
fig.layout.xaxis. dtick = 5
fig.layout.xaxis.title = ''
fig.layout.yaxis.title = ''
fig.show()
```

![newplot (4).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/729216/ff4cb9ba-bc66-7c1c-cfa6-230e69c9e50e.png)

**欠損値**

```python
fig = px.bar(x = df.isnull().sum().index,y= df.isnull().sum().values,title= 'Number of Null values')
fig.layout.xaxis.tickangle = 300
fig.layout.xaxis. dtick = 5
fig.layout.yaxis. dtick = 100000
fig.layout.xaxis.title = ''
fig.layout.yaxis.title = ''
fig.layout.xaxis.showgrid = True
fig.show()
```

![newplot (6).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/729216/2f6b61d3-098e-3eec-b7f8-ce6c5d51da57.png)


欠損値が多いカラムが目立ちますね。データの10%以上が欠損しているカラムを見つけましょう。

```python
nulls = df.isnull().sum()
nulls_list = list(nulls[nulls >(0.1 * len(df))].index)
nulls_list
```

```
['feature_7',
 'feature_8',
 'feature_17',
 'feature_18',
 'feature_27',
 'feature_28',
 'feature_72',
 'feature_78',
 'feature_84',
 'feature_90',
 'feature_96',
 'feature_102',
 'feature_108',
 'feature_114']
```

null値の数にはある種のパターンがあるため、これらの特徴量間の関係を調べます。※参照notebookでは、色が表示されますが、こちらでは、普通のテーブルを使用します。

```python
df[['resp','resp_1','resp_2','resp_3','resp_4','weight']+nulls_list].corr().style.background_gradient(cmap='coolwarm')
```
||resp|resp_1|resp_2|resp_3|resp_4|weight|feature_7|feature_8|feature_17|feature_18|feature_27|feature_28|feature_72|feature_78|feature_84|feature_90|feature_96|feature_102|feature_108|feature_114
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
resp|1.000000|0.452159|0.595274|0.815972|0.956197|-0.006948|0.021207|0.014660|-0.028277|-0.017448|0.048293|0.032502|-0.002034|0.000072|0.011500|0.003604|0.013469|0.004613|0.013209|0.004385
resp_1|0.452159|1.000000|0.890214|0.678381|0.358949|-0.016091|0.033919|0.026142|0.004333|0.006632|0.036125|0.025671|-0.000889|0.001371|0.008330|0.001973|0.010467|0.002994|0.010255|0.002969
resp_2|0.595274|0.890214|1.000000|0.823383|0.469230|-0.012955|0.030800|0.025124|-0.005811|0.000202|0.039829|0.028633|-0.000642|0.001266|0.008487|0.001454|0.010247|0.002223|0.010270|0.002225
resp_3|0.815972|0.678381|0.823383|1.000000|0.805952|-0.009229|0.019317|0.014859|-0.002078|0.000566|0.023373|0.016951|-0.004490|-0.003221|0.007704|0.000655|0.011542|0.003766|0.010231|0.002438
resp_4|0.956197|0.358949|0.469230|0.805952|1.000000|-0.005817|0.015059|0.009246|-0.018851|-0.012306|0.033533|0.022307|-0.003965|-0.002240|0.009916|0.003087|0.012989|0.005355|0.011978|0.004442
weight|-0.006948|-0.016091|-0.012955|-0.009229|-0.005817|1.000000|-0.016164|-0.000819|-0.011541|-0.002568|-0.010643|0.001118|0.002058|0.000367|-0.041597|0.052486|-0.044057|0.056261|-0.046157|0.053893
feature_7|0.021207|0.033919|0.030800|0.019317|0.015059|-0.016164|1.000000|0.898832|0.402987|0.374951|0.487964|0.430844|-0.009719|-0.006903|-0.015400|-0.054359|-0.009083|-0.051083|-0.001724|-0.046770
feature_8|0.014660|0.026142|0.025124|0.014859|0.009246|-0.000819|0.898832|1.000000|0.373778|0.425576|0.424587|0.463257|-0.003916|0.001188|-0.002604|-0.038638|0.000781|-0.039802|0.010609|-0.034254
feature_17|-0.028277|0.004333|-0.005811|-0.002078|-0.018851|-0.011541|0.402987|0.373778|1.000000|0.906956|-0.501928|-0.464884|-0.070201|-0.069432|-0.021351|-0.053571|0.039123|0.007819|0.016894|-0.020790
feature_18|-0.017448|0.006632|0.000202|0.000566|-0.012306|-0.002568|0.374951|0.425576|0.906956|1.000000|-0.446295|-0.503707|-0.078615|-0.074634|-0.021682|-0.054821|0.040938|0.006212|0.017579|-0.023075
feature_27|0.048293|0.036125|0.039829|0.023373|0.033533|-0.010643|0.487964|0.424587|-0.501928|-0.446295|1.000000|0.900539|0.059322|0.060716|0.020922|-0.004243|-0.029925|-0.059175|0.000032|-0.028572
feature_28|0.032502|0.025671|0.028633|0.016951|0.022307|0.001118|0.430844|0.463257|-0.464884|-0.503707|0.900539|1.000000|0.074591|0.075029|0.030016|0.005652|-0.027775|-0.055675|0.006896|-0.021918
feature_72|-0.002034|-0.000889|-0.000642|-0.004490|-0.003965|0.002058|-0.009719|-0.003916|-0.070201|-0.078615|0.059322|0.074591|1.000000|0.773391|0.268906|0.209429|-0.208541|-0.144612|0.029765|0.032632
feature_78|0.000072|0.001371|0.001266|-0.003221|-0.002240|0.000367|-0.006903|0.001188|-0.069432|-0.074634|0.060716|0.075029|0.773391|1.000000|0.221750|0.234574|-0.149852|-0.183428|0.035854|0.024444
feature_84|0.011500|0.008330|0.008487|0.007704|0.009916|-0.041597|-0.015400|-0.002604|-0.021351|-0.021682|0.020922|0.030016|0.268906|0.221750|1.000000|0.645486|0.873262|0.534054|0.951272|0.608618
feature_90|0.003604|0.001973|0.001454|0.000655|0.003087|0.052486|-0.054359|-0.038638|-0.053571|-0.054821|-0.004243|0.005652|0.209429|0.234574|0.645486|1.000000|0.547641|0.891866|0.628328|0.971499
feature_96|0.013469|0.010467|0.010247|0.011542|0.012989|-0.044057|-0.009083|0.000781|0.039123|0.040938|-0.029925|-0.027775|-0.208541|-0.149852|0.873262|0.547641|1.000000|0.633677|0.956765|0.609289
feature_102|0.004613|0.002994|0.002223|0.003766|0.005355|0.056261|-0.051083|-0.039802|0.007819|0.006212|-0.059175|-0.055675|-0.144612|-0.183428|0.534054|0.891866|0.633677|1.000000|0.618610|0.971214
feature_108|0.013209|0.010255|0.010270|0.010231|0.011978|-0.046157|-0.001724|0.010609|0.016894|0.017579|0.000032|0.006896|0.029765|0.035854|0.951272|0.628328|0.956765|0.618610|1.000000|0.645604
feature_114|0.004385|0.002969|0.002225|0.002438|0.004442|0.053893|-0.046770|-0.034254|-0.020790|-0.023075|-0.028572|-0.021918|0.032632|0.024444|0.608618|0.971499|0.609289|0.971214|0.645604|1.000000

これらのカラムのnull値の数は膨大であるため（25万を超える）、これらの特徴量と`Resp`および重みとの相関関係がないため、null値が10％を超える特徴量は取り除きます。

```python
df.drop(columns=nulls_list,inplace=True)
```

残りの欠損値については、最初に変動係数に注目します。

```python
(df.iloc[:,7:-2].std() / df.iloc[:,7:-2].mean()).head(5)
```

```
feature_0    101.635945
feature_1      6.652154
feature_2      6.905722
feature_3    217.145702
feature_4    421.125070
dtype: float64
```

> 平均値がゼロに近いため、変動係数は信頼できないようです


これで、特徴量の分布を俯瞰することができます。

```python
df.iloc[:,7:-2].hist(bins=100,figsize=(20,74),layout=(29,4));
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/729216/ae4391f1-3df2-66ae-2be0-d7a39988badd.png)

次のコードは、特徴量の分布の包括的で確実な理解を得るために、平均のプロットで水平箱ひげ図のグリッドを作成します。

>データが強く中央に偏っているため、0.1％から99.9％のひげを使用して極端な外れ値を取り除いたことに注意してください。

```python
fig = plt.figure(figsize=(20,80))
fig.suptitle('Features Box plot with 0.1% 99.9% whiskers',fontsize=22, y=.89)
grid =  gridspec.GridSpec(29,4,figure=fig,hspace=.5,wspace=.05)
featstr = [i for i in df.columns[7:-2]]
counter = 0
for i in range(29):
    for j in range(4):
        subf = fig.add_subplot(grid[i, j]);
        sns.boxplot(x= df[featstr[counter]],saturation=.5,color= 'blue', ax= subf,width=.5,whis=(.1,99.9));
        subf.axvline(df[featstr[counter]].mean(),color= 'darkorange', label='Mean', linestyle=':',linewidth=3)
        subf.set_xlabel('')
        subf.set_title('{}'.format(featstr[counter]),fontsize=16)
        counter += 1
        gc.collect()
plt.show();
```

![__results___68_0.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/729216/735003d7-73bd-8d5a-066e-01b03070ee49.png)


各特徴量の分布に影響を与える外れ値がたくさんあることがわかります。

また、値の大部分は平均を中心としているため、平均を使用してnull値を補完します。

```python
df.fillna(df.mean(axis=0),inplace=True)
```

**特徴量の成長**

ここでは、日別の平均のプロットの累積和を見て、特徴量の増加について着目します。

```python
df.groupby('date')[featstr].mean().cumsum().plot(layout=(29,4),subplots=True,figsize=(20,82),xlabel='')
fig = plt.gcf()
fig.text(0.5, 0.19, 'Date',ha='center', fontsize = 24)
fig.suptitle('Cumulative features means per day',fontsize=24,y=.886);
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/729216/d76ecbda-1142-40e5-af05-c00a47a3b466.png)


多くの特徴量の日別平均の累積和は直線的に増加しているように見えますが、81、82、83などの一部の特徴量は実際には減少しており、特徴量3のように変動する特徴量もあります。

### 特徴量同士の相関

相関行列のデータフレームを作成します。

```python
corr = df.iloc[:,7:-2].corr()
corr.style.background_gradient(cmap='coolwarm')
```

この相関行列は非常に大きい(130✖︎130)ため、ここでは表示を割愛します。

**相関行列のヒートマップ**

```python
fig = plt.figure(figsize=(18,12))
ax = plt.subplot(1,1,1)
sns.heatmap(corr,ax= ax, cmap='coolwarm');
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/729216/77b6a51e-bdfc-52b8-1632-0b96e853d203.png)

特徴量間に多重共線性がたくさんあるように見えます。また、特徴量空間にカップルのパターンがあるように見えます。このパターンは、`feature_41`などの一部のフィーチャで壊れています。

```python
featstr2 = [ i for i in featstr if i not in ['feature_41','feature_64']]
len(featstr)
```

```
116
```

```python
fig = plt.figure(figsize=(22,44))
grid =  gridspec.GridSpec(12,5,figure=fig,hspace=.5,wspace=.2)
counter = 1
for i in range(12):
    for j in range(5):
        if counter == 113:
            break
        subf = fig.add_subplot(grid[i, j]);
        sns.scatterplot(x= df[featstr2[counter]], y = df[featstr2[counter+1]], ax= subf);
        cor = round(df[featstr2[counter]].corr(df[featstr2[counter+1]]) * 100,2)
        subf.set_xlabel('')
        subf.set_ylabel('')
        subf.set_title('{} & {}\nCorrelation = {}%'.format(featstr2[counter],featstr2[counter+1],cor),fontsize=14)
        counter += 2
        gc.collect();  
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/729216/343f7b9e-eb23-863c-e86f-ea19615f8b5c.png)

>想定通り、これらは財務関連の特徴量であるため、多くの機能は互いに強く相関しています。

次に、相関性の高い特徴量のグループを調査します。

まずは`[feature_19, feature_20, feature_21, feature_22, feature_23, feature_24, feature_25, feature_26, feature_29]`から始めます。なぜなら、多重共線生が見えるからです。

```python
plt.figure(figsize=(12,6)) 
sns.heatmap(df[featstr2[15:23]].corr(),center=0,cmap='coolwarm',annot=True,cbar=False);
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/729216/e774eee3-e564-aa10-c301-3c820a77705e.png)

```python
sns.pairplot(df[featstr2[15:23]],corner=True);
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/729216/5017c6ab-5653-263e-415f-2eccc6b78726.png)

>ピアソンの相関係数がこれらの特徴間で非常に高いという事実にもかかわらず、関係は完全に線形ではありませんが、外れ値が散布図の形状に影響を与えることにも注意できます。

では、他のグループについても見てみましょう。

```python
plt.figure(figsize=(12,6)) 
sns.heatmap(df[featstr2[23:31]].corr(),center=0,cmap='coolwarm',annot=True,cbar=False);
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/729216/8278c5f9-a708-21a6-e075-9c2fdcee6511.png)

```python
sns.pairplot(df[featstr2[23:31]],corner=True);
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/729216/8caa2019-21f9-9abc-cb85-a578b6f1dc08.png)

>要点は、他のクラスターとほぼ同じように見えますが、これらの特徴量クラスターの両方が互いに負の相関関係にあることにも言及する価値があります。

```python
plt.figure(figsize=(18,6)) 
sns.heatmap(df[featstr2[15:31]].corr(),center=0,cmap='coolwarm',annot=True,cbar=False);
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/729216/a8f5da32-4302-ac30-0c46-d0929ff5ced0.png)

**外れ値**

まずは特徴量の平均値について見ましょう。

```python
fig = px.bar(df[featstr].mean(), title='Features mean values')
fig.layout.xaxis.tickangle = 300
fig.update_traces(showlegend = False)
fig.layout.xaxis. dtick = 5
fig.layout.xaxis.title = ''
fig.layout.yaxis.title = ''
fig.show()
```

![newplot (7).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/729216/39b6fb37-e6f7-7d30-6ae8-cc3db603a832.png)

次に最大値について見ましょう。

```python
fig = px.bar(df[featstr].max(), title='Features Max Values')
fig.layout.xaxis.tickangle = 300
fig.update_traces(showlegend = False)
fig.layout.xaxis. dtick = 5
fig.layout.xaxis.title = ''
fig.layout.yaxis.title = ''
fig.show()
```

![newplot (8).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/729216/ec73a887-efb5-bbd9-bdaf-6c2a17e5a5a0.png)

次は最小値です。

```python
fig = px.bar(df[featstr].min(), title='Features Min Values')
fig.layout.xaxis.tickangle = 300
fig.update_traces(showlegend = False)
fig.layout.xaxis. dtick = 5
fig.layout.xaxis.title = ''
fig.layout.yaxis.title = ''
fig.show()
```

![newplot (9).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/729216/8fdce1c5-7e58-ddf2-1aca-7dce1c081cd4.png)

```python
fig, (ax1,ax2,ax3)= plt.subplots(3,1,figsize=(10,12))
plt.subplots_adjust(hspace=.3)
sns.distplot(df[featstr].max(),ax= ax1 )
sns.distplot(df[featstr].min(),ax= ax2)
sns.distplot(df[featstr].mean(),ax= ax3)
fig.suptitle('distribution of mean max and min for features',fontsize=16)
ax1.set_title('distribution  of features max values',fontsize=14)
ax1.text(.82,.56,'std = {}'.format(round(df[featstr].max().std(),2)),transform=ax1.transAxes, verticalalignment='top',bbox=props,fontsize=12)
ax2.set_title('distribution  of features min values',fontsize=14)
ax2.text(.82,.56,'std = {}'.format(round(df[featstr].min().std(),2)),transform=ax2.transAxes, verticalalignment='top',bbox=props,fontsize=12)
ax3.set_title('distribution  of features mean values',fontsize=14)
ax3.text(.82,.56,'std = {}'.format(round(df[featstr].mean().std(),2)),transform=ax3.transAxes, verticalalignment='top',bbox=props,fontsize=12);
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/729216/9d96c95b-3120-ec79-c805-dfa55c2ea171.png)

>平均値は互いにそれほど異ならないという事実にもかかわらず、最小値と最大値は非常に偏った分布で非常にずれていることがわかります。

**外れ値へのより統計的な探索**

```python
for i in featstr[1:]:
    print('{}\n0.1%:99.9% are between: {}\nmax: {}\nmin: {}\n75% are under: {}'.format(i,np.percentile(df[i],(.1,99.9)), df[i].max(),df[i].min(),np.percentile(df[i],75)),
         '\n===============================')
```

```
feature_1
0.1%:99.9% are between: [-3.1720264  16.70666595]
max: 74.42988586425781
min: -3.1720263957977295
75% are under: 1.5784167051315308 
===============================
feature_2
0.1%:99.9% are between: [-3.09318233 16.21772709]
max: 148.07630920410156
min: -3.093182325363159
75% are under: 1.52639901638031 
===============================
feature_3
0.1%:99.9% are between: [-8.56044363  8.42387932]
max: 25.87166404724121
min: -25.424652099609375
75% are under: 1.0469255447387695 
===============================
feature_4
0.1%:99.9% are between: [-7.86522797  7.90286264]
...
```

```python
df[(df.feature_56== df.feature_56.max())|(df.feature_57== df.feature_57.max())|(df.feature_58== df.feature_58.max()) | (df.feature_59== df.feature_59.max())]
```

||date|weight|resp_1|resp_2|resp_3|resp_4|resp|feature_0|feature_1|feature_2|feature_3|feature_4|feature_5|feature_6|feature_9|feature_10|feature_11|feature_12|feature_13|feature_14|feature_15|feature_16|feature_19|feature_20|feature_21|feature_22|feature_23|feature_24|feature_25|feature_26|feature_29|feature_30|feature_31|feature_32|feature_33|feature_34|feature_35|feature_36|feature_37|feature_38|feature_39|feature_40|feature_41|feature_42|feature_43|feature_44|feature_45|feature_46|feature_47|feature_48|feature_49|feature_50|feature_51|feature_52|feature_53|feature_54|feature_55|feature_56|feature_57|feature_58|feature_59|feature_60|feature_61|feature_62|feature_63|feature_64|feature_65|feature_66|feature_67|feature_68|feature_69|feature_70|feature_71|feature_73|feature_74|feature_75|feature_76|feature_77|feature_79|feature_80|feature_81|feature_82|feature_83|feature_85|feature_86|feature_87|feature_88|feature_89|feature_91|feature_92|feature_93|feature_94|feature_95|feature_97|feature_98|feature_99|feature_100|feature_101|feature_103|feature_104|feature_105|feature_106|feature_107|feature_109|feature_110|feature_111|feature_112|feature_113|feature_115|feature_116|feature_117|feature_118|feature_119|feature_120|feature_121|feature_122|feature_123|feature_124|feature_125|feature_126|feature_127|feature_128|feature_129|ts_id|action
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
235877|38|3.339560|-0.001017|-0.007426|-0.024829|-0.046755|-0.033357|-1|5.399380|5.108447|3.613957|3.282637|1.874917|1.627905|3.537130|2.426174|2.572712|3.605390|-3.484599|-2.825172|2.857596|3.233086|-3.352329|-3.912990|-12.886112|-10.339006|-4.377385|-5.277115|-3.872348|-4.283693|3.506302|4.035473|3.110046|4.478292|3.713556|3.731210|4.706809|4.334348|-3.058521|-3.192860|3.088272|4.539015|-0.136437|3.195502|0.780103|1.167637|-0.327845|0.913224|1.520857|1.316261|1.053037|0.099054|0.855302|1.775148|3.746835|0.985716|388.749817|66.636330|660.351135|800.703613|350.911865|-0.671831|-0.569680|-1.486158|-1.611610|-2.942103|-1.607305|-1.071010|-1.565248|-1.505246|2.665544|-0.095570|1.875405|0.898643|0.415992|1.163417|0.861380|0.671799|0.222680|0.190999|0.681430|0.211753|0.399736|1.027329|0.849244|1.615301|0.705474|-1.423562|-0.041847|-0.107219|-0.006064|-0.310370|-2.582340|-1.338859|0.884709|-0.057046|-1.194013|-1.427599|-0.940190|-0.260405|-0.789533|-1.781693|-2.874468|-0.429607|0.747440|0.427909|-0.787168|-1.908855|-1.880670|-0.350742|-0.874343|-2.260690|-1.953951|6.032958|4.065476|3.634447|1.015100|5.885513|4.415474|8.390346|4.077566|5.459621|2.131100|235877|0
282230|44|10.859674|-0.001306|-0.012651|-0.016530|-0.052089|-0.048993|1|7.414509|9.399608|5.123063|6.146243|4.196547|4.889341|11.098932|11.976443|2.739184|4.887201|6.628507|7.351760|8.330588|11.728739|9.664564|13.641744|5.919429|5.534737|6.929136|9.939451|9.581731|13.009244|-2.449863|-3.702712|-1.050895|-2.082486|-1.886110|-2.196391|-3.542542|-3.657356|5.737075|7.694194|-3.350399|-6.203264|-0.113139|1.447529|-0.254701|7.538714|1.115668|2.226885|3.178381|4.724526|6.178386|1.928744|3.180385|-0.291821|3.230325|2.360054|746.739868|254.094452|816.712097|138.990021|55.821522|-1.349714|-1.129695|-0.094330|-0.107096|-2.041466|-0.597949|-0.441602|-1.366051|-1.314805|16.274912|13.676999|19.838009|-5.815048|-2.006737|-3.691047|-5.506912|0.018003|-2.335785|-0.879096|-2.030294|-2.191662|0.140086|-1.238212|-0.142781|-2.321156|-1.641860|-1.393412|-1.515613|-1.216623|-1.746285|-1.086886|-2.817912|9.008026|2.096552|3.640483|6.876660|-1.182548|3.427998|0.083712|1.071575|3.375455|-2.786946|6.497206|0.910018|1.932831|5.401005|-1.713007|2.378500|-0.578150|0.342619|1.677919|-2.013862|9.998623|10.710449|3.118533|1.845357|6.339649|7.125723|4.243608|2.327895|3.403054|1.688110|282230|0
286937|44|0.219883|0.009060|0.000275|0.008128|0.006852|-0.004134|1|2.464552|3.314075|6.717828|8.427411|5.999205|7.184940|7.984714|8.634279|10.861753|16.884825|16.154589|19.937834|22.595144|31.915413|10.789083|15.700766|14.886054|16.765507|17.146122|24.403795|24.710970|34.675587|-7.076792|-10.310453|-5.595434|-9.575839|-9.691500|-12.738782|-16.241947|-20.111063|7.307025|10.078041|-1.336552|-2.614510|0.127066|0.966415|-0.375578|8.362889|1.284996|-1.139717|-1.475705|-1.193366|-0.914533|-2.279577|-1.075379|-1.819851|8.140661|8.705054|161.595047|2228.305664|175.472198|234.087051|983.675720|-0.342254|-0.294564|1.821556|1.884842|-0.513369|0.651058|0.537046|-2.615676|-2.487624|4.619804|4.233167|5.308126|-5.576128|-2.581536|-4.365026|-6.566184|-1.339482|-2.204539|-1.126629|-2.361625|-2.628877|-0.362707|0.566276|-0.091297|0.440614|1.080030|-0.755546|-0.407814|-1.183697|-0.848132|-0.313204|-2.134031|8.850863|2.711628|4.811763|8.958502|-0.366285|3.318763|0.329840|1.538568|4.348969|-1.816757|6.563427|1.253030|2.842146|7.301911|-0.718124|2.376652|-0.260213|0.996117|2.531418|-1.418708|2.053362|1.714995|5.849802|5.440816|1.884865|1.741266|2.375170|0.947722|3.170853|1.950660|286937|0
834854|168|0.596468|0.002326|0.002032|0.005509|0.007654|0.003616|-1|-1.436731|-2.133809|4.648743|2.969781|0.907301|0.450538|1.191757|-0.143514|9.682172|8.077853|18.023613|10.920301|25.319948|17.384563|-1.470224|-1.229787|-3.497435|-1.578663|-4.015155|-3.872986|-5.276647|-4.489394|4.190367|3.187598|10.996456|9.465394|19.159149|12.718694|31.042587|19.925367|-0.962800|-0.650292|4.028402|4.142321|0.611631|1.429632|2.867671|-2.501698|-0.290131|-0.013775|-0.325107|-0.721573|-0.889294|-0.459434|-0.440759|-0.153293|7.292922|6.501546|23.773369|1444.840088|87.525017|303.269165|1372.886841|1.325848|1.173081|-0.216417|-0.236258|1.529365|2.043775|1.743800|3.031896|2.989350|-0.761621|-2.308117|-1.328938|0.220087|0.031546|0.051950|0.083291|0.902025|0.278071|0.052404|0.093597|0.077576|1.635452|-0.001564|-1.001347|-1.168230|-0.163837|-0.383313|0.311700|-1.112664|-0.479528|0.244936|0.945099|-0.420525|-1.712330|-0.716374|-0.219099|-0.496815|-0.363410|-0.959098|-0.536406|0.115152|0.687828|-1.948212|-1.659458|-1.141388|-1.382763|-0.595652|-1.113544|-1.988687|-1.147521|-0.065230|0.619337|0.335134|0.268782|7.619550|7.574471|0.335153|0.244880|0.339170|0.232385|8.562437|7.662621|834854|1

>データセットには極端な外れ値があると推測できます。また、一部の外れ値には、特徴間の高い多重共線性の結果である隣接列からの大きな値が付随していることにも言及する価値があります。

次に、特徴量データの99.9％を超えるこれらの極端な外れ値を取り除いてみましょう。

データセットのループ中にデータが削除されないようにするために、すべての特徴量に対して99.9％をマークするためのリストを作成します。また、後で調べるために「.1％のマークを使用して」という負の外れ値のリストを作成します。

```python
n999 = [ np.percentile(df[i],99.9) for i in featstr[1:]]
n001 = [ np.percentile(df[i],.1) for i in featstr[1:]]
```

```python
for i, j in enumerate(featstr[1:]):
    df = df[df[j] < n999[i]]
    gc.collect()
```

どの程度のデータを削除したか見ましょう。

```python
str(round(((org_len - len(df))/org_len)*100,2))+'%'
```

```
'4.01%'
```

```python
fig = px.bar(df[featstr].max(), title='Features Max Values')
fig.layout.xaxis.tickangle = 300
fig.update_traces(showlegend = False)
fig.layout.xaxis. dtick = 5
fig.layout.xaxis.title = ''
fig.layout.yaxis.title = ''
fig.show()
```

![newplot (10).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/729216/3c945d73-8490-651a-0396-a9d2b892a3fe.png)

次に、カスタマイズされた.1％と99.9％のひげを使用して箱ひげ図グリッドを再度作成します。

```python
fig = plt.figure(figsize=(20,80))
fig.suptitle('Features Box plot with 0.1% 99.9% whiskers',fontsize=22, y=.89)
grid =  gridspec.GridSpec(29,4,figure=fig,hspace=.5,wspace=.05)
counter = 0
for i in range(29):
    for j in range(4):
        subf = fig.add_subplot(grid[i, j]);
        sns.boxplot(x= df[featstr[counter]],saturation=.5,color= 'blue', ax= subf,width=.5,whis=(.1,99.9));
        subf.set_xlabel('')
        subf.set_title('{}'.format(featstr[counter]),fontsize=16)
        counter += 1
        gc.collect();
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/729216/5bb8aa92-0c54-6cac-61f0-e4d6dac9aa3e.png)

>正の外れ値のみを削除したため、データセットにはまだ多くの明らかな外れ値、特に負の値があります。そのため、以前は対称的な外れ値を持っていた特徴量に、ある種の左に歪んだ外れ値があります。
>両方の箱ひげ図を比較すると（正の外れ値を削除する前と後）、以前は対称的な外れ値を持っていた`feature_3`から`feature_40`までの特徴量が、トリミング後に極端な負の外れ値を持つように変更されていることがわかります。

手動でこれらの特徴量をトリミングします。

```python
for i,j in zip(featstr[1:][2:34],n001[2:34]):
    df = df[df[i] > j]
    gc.collect();
```

次のコードは、これまでのクリーニングのために失われた観測値の比率を示しています。

```python
str(round(((org_len - len(df))/org_len)*100,2))+'%'
```

```
'4.89%'
```

外れ値を削除し、`resp`値を考慮した後の特徴量のkde(カーネル密度推定)とヒストグラムを見ましょう。

```python
fig = plt.figure(figsize=(20,80))
fig.suptitle('KDE plot of Features',fontsize=24,transform =fig.transFigure, y=.89)
grid =  gridspec.GridSpec(29,4,figure=fig,hspace=.5,wspace=.01)
counter = 0
for i in range(29):
    for j in range(4):
        subf = fig.add_subplot(grid[i, j]);
        sns.distplot(df[df.action==0][featstr[counter]],bins= 100,label='Negative',
                     color='darkorange', kde_kws={'linewidth':4},ax=subf)
        sns.distplot(df[df.action!=0][featstr[counter]],bins= 100,label='Positive',
                     color='blue', kde_kws={'alpha':.9,'linewidth':2},hist_kws={'alpha':.3},ax=subf)
        subf.axvline(np.percentile(df[featstr[counter]],99.5),color= 'darkblue', label='99.5%', linestyle=':',linewidth=2)
        subf.axvline(np.percentile(df[featstr[counter]],.5),color= 'red', label='0.5%', linestyle=':',linewidth=2)
        subf.legend().set_visible(False)
        subf.set_xlabel('')
        subf.set_title('{}'.format(featstr[counter]),fontsize=16)
        kurt=round(df[featstr[counter]].kurt(),2)
        skew=round(df[featstr[counter]].skew(),2)
        subf.text(.6,.92,'Kurt = {:.2f}\nSkew = {:.2f}'.format(kurt ,skew),
         transform=subf.transAxes, verticalalignment='top',bbox=props,fontsize=10)
        counter += 1
        gc.collect();
handles, labels = subf.get_legend_handles_labels()
fig.legend(handles, labels,ncol=4, bbox_to_anchor=(0.86, 0.893),fontsize=10,
           title= 'Resp',title_fontsize=14,bbox_transform =fig.transFigure);
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/729216/defedc3a-ae5e-8d1e-3670-85a30c1c1d89.png)

>特徴量の分布に`resp`値を追加すると、次のようになります。
> - 特徴のヒストグラムは、外れ値がはるかに少なくなり、より正式な分布になりました。
> - また、1、2、85、87、88、91などの一部の機能には多くの負の外れ値があることがわかります。
> - 49、50、51、55、56、57、58、59などの一部の機能には、まだ多くの正の外れ値があります。
> - 特徴量の分布は、`resp`値の影響を受けません。

**特徴量と`resp`の相関**

まず、`resp`と各特徴量の関係についてpandasの相関系列を作成します。

```python
respcorr =  pd.Series([ df.resp.corr(df[i]) for i in featstr],index=featstr)
```

```python
fig = px.bar(respcorr,color = respcorr, color_continuous_scale=['red','blue'], title= 'Features Correlation with Resp')
fig.layout.xaxis.tickangle = 300
fig.layout.xaxis. dtick = 5
fig.layout.xaxis.title = ''
fig.layout.yaxis.title = 'pearson correlation'
fig.update(layout_coloraxis_showscale=False)
fig.show();
```

![newplot (11).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/729216/5f6c995f-8d50-1de4-b1ac-e651237969a2.png)

特徴量が`resp`とそこまで高い相関がないことがわかります。

**特徴量と重み(`weight`)**

ここでも、pandasの一連の重みと特徴量の相関を見ますが、0より大きい重みのみを考慮します。

```python
wecorr = pd.Series([df[df.weight != 0].weight.corr(df[df.weight != 0][i]) for i in featstr],index=featstr)
```

```python
wecorr.head(10)
```

```
feature_0    -0.004811
feature_1    -0.039106
feature_2     0.061847
feature_3     0.000140
feature_4     0.005604
feature_5    -0.001310
feature_6     0.001633
feature_9    -0.096143
feature_10    0.012905
feature_11   -0.041127
dtype: float64
```

```python
fig = px.bar(wecorr,title= 'Features Correlation with Weight (not including zero weights)')
fig.layout.xaxis.tickangle = 300
fig.layout.xaxis. dtick = 5
fig.layout.xaxis.title = ''
fig.layout.yaxis.title = 'pearson correlation'
fig.update(layout_coloraxis_showscale=False)
fig.update_layout(showlegend=False)
fig.show()
```

![newplot (12).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/729216/be3ca923-65b5-c20b-79b1-1ee6198811ea.png)

散布図を使用して、相関係数の最大値および最小値を調べます。

まずは、重みに対して最も高い相関を持つ特徴量(feature_51)について見ます。

```python
fig = plt.figure(figsize=(8,6))
sns.scatterplot(df[df.weight != 0].weight,df[df.weight != 0].feature_51, color = 'darkblue', alpha=.3)
plt.xlabel('Weight',fontsize=14)
plt.ylabel('Featre_51',fontsize=14)
plt.title('Feature_51 vs Weight\nCorrelation = {}%'.format(round(df[df.weight != 0].weight.corr(df[df.weight != 0].feature_51),4)*100),fontsize=16);
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/729216/808965f1-16ba-2a3f-5326-9cd1a53b320f.png)

重み(`weight`)はfeature_51と散布図からも高い相関があることがわかります。

重みに最も低い相関を持つ特徴量(feature_126)について見ましょう。

```python
fig = plt.figure(figsize=(8,6))
sns.scatterplot(df[df.weight != 0].weight,df[df.weight != 0].feature_126, color = 'darkblue', alpha=.3)
plt.xlabel('Weight',fontsize=14)
plt.ylabel('Featre_126',fontsize=14)
plt.title('Feature_126 vs Weight\nCorrelation{}%'.format(round(df[df.weight != 0].weight.corr(df[df.weight != 0].feature_126),4)*100),fontsize=16);
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/729216/46b72ae9-8951-711c-70af-814bcd54a23a.png)

>重みとfeature_126の間にはある種の負の相関関係がありますが、関係は弱いようです

**Feature 0**

feature_0の一意の値を見つけましょう。

```python
plt.figure(figsize=(7,5)) 
df.feature_0.value_counts().plot.bar(color='darkblue',alpha=.6,width=.5)
plt.title('Feature_0',fontsize=18) 
plt.xticks(rotation=0,fontsize=14);
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/729216/1cdb44e0-000f-133a-12a9-0872efcf3d6e.png)

>Feature_0はバイナリな特徴量のようです

`Resp`を考慮する。

```python
plt.figure(figsize=(8,6)) 
sns.countplot(data=df, x='feature_0', hue='action',palette='viridis')
plt.xticks(fontsize=14)
plt.yticks(fontsize=13)
plt.xlabel('Feature 0',fontsize=12)
plt.title('Feature 0 and Resp', fontsize=18)
plt.ylabel('')
plt.xlim(-1,2)
h, l = plt.gca().get_legend_handles_labels()
plt.legend(h,['Negative','Positive'],ncol=1, fontsize=12, loc=3,title= 'Resp',title_fontsize=14);
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/729216/c991704f-ff86-7cb4-b540-f72d0678f8a1.png)

>feature_0とrespが負または正である間に明確な関係がないようです