Obey Your MATHEMATICS.

機械学習関連の純粋数学や実験など

つくってあそぼう!ヒッグス粒子発見器 -データ前処理編-

こんにちは。

皆さん、ヒッグス粒子をご存知でしょうか。


ヒッグス粒子 - Wikipedia


2013年に、スイスにあるCERN加速器実験で発見した、と発表されたばかりの新しい粒子です。*1

なんでも、”質量はどこから生まれるのか”みたいな問いに答えを与えるような粒子だそうで、なにやら凄いみたいです。*2



この記事のシリーズではUCI Machine Learning Repositoryにある、Higgs粒子発見器を作るためのデータセット


UCI Machine Learning Repository: HIGGS Data Set


を用いて、Higgs粒子発見器(!)を作って行こうと思います。


第一回目の今回は、前処理に重点を置いていきます。

§1. データセットの概要

データの概要ですが、とりあえず引用させてもらいます:

The data has been produced using Monte Carlo simulations. The first 21 features (columns 2-22) are kinematic properties measured by the particle detectors in the accelerator. The last seven features are functions of the first 21 features; these are high-level features derived by physicists to help discriminate between the two classes. There is an interest in using deep learning methods to obviate the need for physicists to manually develop such features. Benchmark results using Bayesian Decision Trees from a standard physics package and 5-layer neural networks are presented in the original paper. The last 500,000 examples are used as a test set.*3

つまり、

・各サンプルは、ヒッグス粒子の理論値のモンテカルロシミュレーション or ノイズデータ

・サンプル数は11,000,000個。

・目的はノイズ信号とヒッグス粒子の信号を区別する識別器を作ること。

・変数は全部で28個

・最初の21個は加速器により検出される数値

・残りの7個は最初の21個の変数の関数の値で、物理学者によって作られた分類を助けるためのデータ

・最後の50万個のサンプルをテストに使いましょう。

と言った感じです。


さて、

https://archive.ics.uci.edu/ml/machine-learning-databases/00280/

ここから圧縮されたデータをダウンロードして解凍しますと、8GB程の.csvが入っています。*4

これを使って前処理をしてきましょう。

§2. データセットの分割

今回は1100万個のデータ全てではなくて、マシンのキャパの関係上250万個のデータを使っていきます。
本来の設定と同じ数である50万個のデータでテストは行うのでまあ良しとしましょう。

まずデータを読み込み、200万行のトレーニング用のcsvと50万行のテスト用csvを別々に作ります。

import pandas as pd
import numpy as np
import csv as csv
all_df = pd.read_csv("HIGGS.csv", header=0, nrows=2500000)
train_df = all_df[:2000000]
test_df = all_df[2000000:]
train_df.to_csv('Higgs_train.csv')
test_df.to_csv('Higgs_test.csv')

次に、特徴選択をします。上述の通りこのデータセットには

21個のlow-level feature

7個のhigh-level feature

が存在します。またこのデータセットを提供した人の論文に

[1402.4735] Searching for Exotic Particles in High-Energy Physics with Deep Learning

・全ての特徴を使った分類器
・low-levelのみの分類器
・high-levelのみの分類器

それぞれのベンチマークが載っています。今回はhigh-level featureのみを使った分類機を作ろうと思いますので、それ以外の列を落とした.csvを作ります。

#読み込み
train_df = pd.read_csv("Higgs_train.csv", header=None)
test_df = pd.read_csv("Higgs_test.csv", header=None)

#第0列は無駄なIDなので削除 / 1はHiggs or noiseなので残す/ 2~22はlow featureなので削除 / 23~29はhigh featureなので残す
drop_n_list = range(2,23)
drop_n_list.insert(0, 0)
train_df =train_df.drop(drop_n_list,axis=1)
test_df =test_df.drop(drop_n_list,axis=1)

#保存
train_df.to_csv('high_feature_train.csv')
test_df.to_csv('high_feature_test.csv')


これで”前処理するための”準備が整いました。

§3. 欠損値の有無の確認

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

import pandas as pd
import numpy as np
import csv as csv
import matplotlib.pyplot as plt
import seaborn as sns

train_df = pd.read_csv("high_feature_train.csv", header=0, index_col=0, names=['is_Higgs','m_jj', 'm_jjj', 'm_lv', 'm_jlv', 'm_bb', 'm_wbb', 'm_wwbb'])

次に欠損値を確認& "is_Higgs"をカテゴリー変数に変換します。

#欠損値確認
len(train_df) - train_df.count()

#is_Higgs --> 0 or 1 のカテゴリー変数に変換。
train_df["is_Higgs"] = train_df["is_Higgs"].astype('category', categories=[1,0], ordered=False)

f:id:mathetake:20170113210121p:plain

と言う出力が得られると思います。欠損値なしです。

テストデータに関して同様の操作をし、欠損値がないことも確認できます。

§4. データの概要&視覚化

欠損値が無いことが確認できたので、基本統計量を見たり視覚化したりしましょう。

まず、基本統計量の確認から:

#基本統計量。ラベルは無視。
train_df.describe()

#ラベル毎の頻度
train_df['is_Higgs'].value_counts()

f:id:mathetake:20170113210851p:plain

f:id:mathetake:20170113211327p:plain

こんな感じの出力が得られると思います。5%程のラベル間のサンプル数の差がありますが(不均衡データ)、今回は無視してそのまま進みます。


次に散布図行列を描いてみましょう。is_Higgs =0 or 1で色分けします。

#散布図行列
sns.pairplot(train_df, hue='is_Higgs')

f:id:mathetake:20170113211912p:plain


最後に箱ひげ図を見てみます。

# 箱ひげ図
fig, axes = plt.subplots(2, 1)
for i, (n, g) in enumerate(train_df.groupby('is_Higgs')):
    sns.boxplot(data=g.iloc[:, 1:], ax=axes[i])
    axes[i].set_ylabel(n)

f:id:mathetake:20170113212406p:plain

これを見るとなんとなーーく外れ値がちらほらありそうな気がします。

外れ値検出しましょう。

§5. 外れ値検出

今回はIsolationForest を使って外れ値検出します。

箱ひげ図を見た限り、is_Higgsのクラス毎に外れ値のスケールが結構変わってきそうなので別々に処理します。

#  Scikit-learnのIsolationForestを使います。
from sklearn.ensemble import IsolationForest

#is_Higgs =1 or 0 でわけでarrayで格納
train_Higgs = train_df.loc[(train_df.is_Higgs == 1 )]
train_array_Higgs =train_Higgs.values

train_noise = train_df.loc[(train_df.is_Higgs == 0 )]
train_array_noise =train_noise.values 

#それぞれ分類機を作って、結果を'is_inlier'と言うkey名でデータフレームで保存
clf_Higgs = IsolationForest(n_estimators=50,max_samples=10000,contamination=0.00002)
clf_Higgs.fit(train_array_Higgs)
train_Higgs_SVM_result = clf_Higgs.predict(train_array_Higgs)
train_Higgs_SVM_result = np.reshape(train_Higgs_SVM_result,(-1,1))
train_Higgs_SVM_result = pd.DataFrame(train_Higgs_SVM_result, columns=['is_inlier'])

clf_noise = IsolationForest(n_estimators=50,max_samples=10000,contamination=0.00002)
clf_noise.fit(train_array_noise)
train_noise_SVM_result = clf_noise.predict(train_array_noise)
train_noise_SVM_result = np.reshape(train_noise_SVM_result,(-1,1))
train_noise_SVM_result = pd.DataFrame(train_noise_SVM_result, columns=['is_inlier'])

#train_Higgs/noiseに上の結果を結合 (indexのリセットをしないと結合がうまくいかないので注意。)
train_Higgs = pd.concat([train_Higgs.reset_index( drop = True ), train_Higgs_SVM_result.reset_index( drop = True )], axis=1)
train_noise= pd.concat([train_noise.reset_index( drop = True ), train_noise_SVM_result.reset_index( drop = True )], axis=1)

IsolationForestのパラメータですが、分類器の数はデフォルトの50で、各分類機に投入されるサンプルの数は1万、contaminationは0.0002としています。
contaminationは外れ値の割合のことですが、ざっくり0.02%としました。箱ひげ図を見た感じで決めました。


それぞれのクラスでどんな値が異常値として検出されたのか、散布図行列で見てみましょう。

#それぞれどんな値がoutlierになったか可視化しよう。重要そうな4つの成分だけ見る。
sns.pairplot(train_Higgs[['m_jj','m_jjj','m_jlv','m_bb','is_inlier']], hue='is_inlier')
plt.savefig("Higgs_outlier.png")
sns.pairplot(train_noise[['m_jj','m_jjj','m_jlv','m_bb','is_inlier']], hue='is_inlier')
plt.savefig("noise_outlier.png")

まずHiggs粒子のクラスタ:
f:id:mathetake:20170113220056p:plain

次にnoiseのクラスタ;
f:id:mathetake:20170113220134p:plain


青い点が検出された外れ値になります。なんか上手い具合に外れたものが検出されている(ような)気がします。
検出された外れ値を除外&その数を計算してみましょう:

#それぞれoutlierの行('is_inlier' == -1)を落とす。
n = train_Higgs.shape[0]
train_Higgs = train_Higgs.loc[(train_Higgs.is_inlier == 1 )].reset_index( drop = True )
print "Higgs outliers :",  n - train_Higgs.shape[0]

n = train_noise.shape[0]
train_noise = train_noise.loc[(train_noise.is_inlier == 1 )].reset_index( drop = True )
print "Noise outliers : ", n - train_noise.shape[0]

f:id:mathetake:20170113220745p:plain

それぞれ20個程度の外れ値が検出された事が分かります。

除外して得られたクラス毎のデータを結合して、§4のように可視化していきましょう。

#除外後のデータを1つのDataFrameに結合
train_df = pd.concat([train_noise, train_Higgs], ignore_index=True).drop("is_inlier", axis=1)

#除外後の基本統計量
train_df.describe()

#散布図
sns.pairplot(train_df, hue='is_Higgs')
plt.savefig("00002RF_pair_plot.png")

#箱ひげ図
fig, axes = plt.subplots(2, 1)
for i, (n, g) in enumerate(train_df.groupby('is_Higgs')):
    sns.boxplot(data=g.iloc[:, 1:], ax=axes[i])
    axes[i].set_ylabel(n)
plt.savefig("00002RF_box_plot.png") 

◯基本統計量
f:id:mathetake:20170113221028p:plain
こんな感じです。除外前とほとんど変わりませんが、'm_jj'と'm_jjj'のmaxがだいぶ低くなっている事が分かります。

◯散布図行列
f:id:mathetake:20170113221412p:plain

◯箱ひげ図
f:id:mathetake:20170113221459p:plain
箱ひげ図を除外前と見比べるとかなりいい感じに外れ値が除外されているのが見て取れると思います。
(若干'm_jj'の外れ値っぽいのが残っていますが、2個/10万個 程度なので無視しましょう…。)


最後に外れ値を除外したDataFrameを.csvで保存して前処理完了とします。

train_df.to_csv("outlier_removed_high_feature_train.csv")


以上で

つくってあそぼう!ヒッグス粒子発見器 -データ前処理編-

はおわりです。

次回、実際に上で作ったデータセットを元に分類機を作っていきます。それでは。

*1:正確には実験結果がヒッグス粒子である事を強く示唆していると発表した

*2:僕は物理畑の人間ではないので、詳しくは知りません。

*3:UCI Machine Learning Repository: HIGGS Data Set

*4:2.6Gあります。