トポロジカルデータアナリシス：データと図形を結びつける

導入

不定期でトポロジカルデータアナリシス(TDA)に関する紹介をします。今回の内容は、データから図形を見立てる方法と、TDAにて登場する単体的複体の紹介をします。

データを図形に読み替える

以下の図の左側のような点の集まりを、データの集まりだとします。次に、点の周りにある半径の球を描きます*1。すると、球どうしで重なるものがでます。重なったものは、おおざっぱに「つながっている」ち見なします。すると、データは右側のような図に見立てられます。

この見立てた図から何か情報が引き出せれば、それはデータの持つ大事な性質と関係しているかもしれません。

問題点

情報の引き出し方はいずれ紹介するとして、実際に引き出すにあたって以下の問題があります。

1. 見立てた図形をコンピュータで解析することは難しい

見立てた図は人間が「何となくこんな感じ」と考えたもので、細かくみると変てこりんな形をしています。そのため、こういった複雑な形のうまく計算するように実装するのは大変です。そのため、実際の計算では工夫が必要になります。

実際に解析するのは抽象単体複体

実際に計算するために有効なのが抽象単体複体です。大雑把な流れとしては、点の周りにある領域を考え、つながり方から抽象単体複体、およびその幾何学的実現（図形）*2を構成します。これは、点と線や三角形といった図形の組み合わせとなっており、コンピュータで計算することが比較的容易です*3。つまり、見立てた図形を何らかのルールによって三角形的な図形で分割して近似する、ということを（大雑把には）やっています。何らかのルールの部分は幾つか種類がありえますが、以下で代表的な3つの方法を紹介します。

チェック複体

最も素直な構成の仕方です。

チェック複体を作る流れ

1. 点の周りにある半径の球を描く
2. 球が接していたら「つながっている」とみなし、抽象単体複体（チェック複体）を構成する
   • 抽象単体複体を構成するルール例（ざっくり）
     • 2つの球が重なる：2つの点でラベルされる1単体が抽象単体複体の要素
     • 3つの球が重なる：3つの点でラベルされる2単体が抽象単体複体の要素

ヴィートリス・リップス複体

チェック複体は、最も素直な構成の仕方ですが、計算量が多いという欠点があります。それは、何個の球が重なっているかをすべて洗い出す必要があるからです。考える点の数が増えると、考慮する場合の数も爆発的に増えてしまいます。そこで、チェック複体の簡易版として、以下のヴィートリス・リップス複体というものも考えられています。

ヴィートリス・リップス複体を作る流れ

1. 点の周りにある半径の球を描く
2. 球が接していたら「つながっている」とみなし、抽象単体複体（ヴィートリス・リップス複体）を構成する
   • チェック複体との違い
     • 個の球の「任意の」2点が重なっていれば、単体が抽象単体複体の要素
     • チェック複体は、個の球すべてが重なっている領域がないとダメ

注意点

あくまでチェック複体の簡易版であり、もとの調べたい図形と性質が変わる場合があります。以下に、その具体例を載せます。左側は穴が一つ空いている図と見なせそうですが、そこからチェック複体の幾何学的実現を構成したもの（右側）は、穴の空いていない図になります。TDAでは穴の数が本質的な量であると見なしており、両者の性質は異なることになります。

アルファ複体

ちょっと構成法が複雑ですが、計算量がチェック複体ほど多くならなそうな単体的複体です*4。

アルファ複体を作る流れ

1. 点と点の間を垂直2等分線で分割する。ただし線を引くのは他の垂直2等分線とぶつからない部分だけ
2. 各点がボロノイ領域というもので分割される*5
3. さらに点の周りにある半径の球を描き、ボロノイ領域と球の共通部分を考える
4. 上で考えた共通部分どうしが接していたら「つながっている」とみなし、抽象単体複体（アルファ複体）を構成する
   • ルールはチェック複体と同じ

なぜ抽象単体複体を調べればよいのか

これまで、TDAで登場する単体的複体についてざっと紹介しました。しかし、なぜこれらの単体的複体で近似すればよいかについては何も言っていませんでした。実は、点の周りに描く領域を凸集合に（球は凸集合）し、個の球が重なっているときに単体が要素になる（チェック複体とアルファ複体）ような抽象単体複体は、もとの見立てた図と定性的な性質が同じである*6ことが保証されています。保証する定理は、脈体定理といいます。この脈体定理があるので、TDAで知りたいような「定性的な量」は

1. いったん抽象単体複体の幾何学的実現の世界にいき、そこで必要な量を求める
2. 脈対定理より、抽象単体複体で求めた量が見立てた図で必要だった量と等しい

という手順で求めればよいということになります。ただし、ヴィートリス・リップス複体に関しては、脈対定理の前提を満たしていません。なので必ずしももとの見立てた図とヴィートリス・リップス複体の幾何学的実現の性質は一致しないことに注意です。

脈体定理に関しては、例えば
www.kyoritsu-pub.co.jp
に概要が載っています。証明までは載っていないので、気になる方は位相幾何学の本を読むとよいと思います*7。

残された疑問点

今回、データから図形を見立てる方法と、トポロジカルデータアナリシス（TDA）にて登場する単体的複体の紹介をしました。しかしながら、以下のような疑問点がまだ残っています。

• 見立てた図形から、何がわかるのか
• 点の周りに描く円の半径を変えたら、見立てられる図形も変わるけど、それは何を意味しているのか

一つ目を知るには、ホモロジー群、二つ目を知るにはパーシステントホモロジー群という数学的な対象を知る必要があります。これらをつぶさに紹介するのは非常に大変なため、詳しい説明は省略します*8。一つ目は、大雑把には図形の穴の数を抽出するイメージです*9。
tekenuko.hatenablog.com

二つ目は、TDAにてノイズ込みのデータからどうやって情報を抽出するかの部分と関係するものです。次回は、その工夫の仕方をざっくり紹介しようと考えています。

*1:必ずしも球である必要はないです。

*2:幾何学的実現は、抽象単体複体の要素の和集合（つまり、くっつけれる単体を全部くっつけた図形)です。

*3:ここの言い回しは、実際のアルゴリズムやその他の手法などを自分が理解していないため、不正確かもしれません。

*4:これは、要素として現れうる単体の次元が、アルファ複体の方が小さいことが理由の一つだと自分は思っています。実際の計算量も含めて細かいところは理解していないものがあるので、理解が深まったら都度紹介していってもアリかなと思っています。

*5:ボロノイ領域が接していた場合に「つながっている」とみなし、抽象単体複体を構成することもできます。この抽象単体複体をドロネー複体といいます。

*6:ホモトピー同値という性質です。

*7:自分は証明したことありません。。

*8:自分もきちんと理解しているかは疑問です（汗）。

*9:穴の数がわかったから何なんだ、という疑問は大いにあると思います。このあたりもそのうちうまく紹介したいものです。。

トポロジカルデータアナリシス : 単体的複体

導入

不定期でトポロジカルデータアナリシス(TDA)に関する紹介をします。今回は、私たちが普段から何気なく慣れ親しんできた図形に関する話題です。ここでは、単体的複体をいう対象を考えます。

参考

このブログでは、数式を交えた説明は基本的にしません。きちんと理解したい方は、例えば
www.kyoritsu-pub.co.jp
があります。ただし、数学科の2〜3年生くらいの知識がある程度無いと読み進めるのは難しいです。

もう少し数学科以外に向けた本としては
www.seshop.com
があります。素粒子論をやりたい学生が、何故か（ほぼ）皆が申し合わせたように学部生あたりのときに読んでいる本です。微分積分学や線形代数を理解していれば読むことは可能です。

単体

ざっくり定義

単体的複体を考えるために、その構成要素である単体から考えています。
単体とは、点（0次元）や辺（1次元）、三角形（2次元）、四面体（3次元）といった図形を一般化したものです。これらをそれぞれ0次元三角形、1次元三角形、...と呼ぶことにします。これは、次元空間内の次元三角形という拡張ができます。この次元三角形を単体といいます。もう少し正確には、次元空間内に個の点を配置し、線形独立なベクトルを与え、その凸包を考える、という定式化をしますが、詳しくは参考書を参考にしていただけばと思います。

単体はより次元の低い単体を一部に含みます。それを「面」といいます。例えば、四面体（3単体）は三角形（2単体）*1や辺（1単体）、点（0単体）を一部に含んでいます。

単体的複体(複体)

ゆるゆる定義

単体的複体は、単体をあるルールによって集めた集合です。あるルールというのは、大雑把には以下の2つです。

• ある単体的複体に含まれる単体の面（部分集合的なもの）も単体的複体のメンバー
• ２つの単体が共通部分をもつときは、それぞれから見てきちんと面になっている部分で共有している

以下に、単体的複体を構成できる図形とそうでない例をあげます*2。

参考（http://infoshako.sk.tsukuba.ac.jp/~hachi/math/csc2/complex.html）

左側は三角形同士が点や辺の部分で貼り合わされていますが、右側ではそうではないものがあります。例えば、辺と点がくっついているもの（右側の図形のうち、左側）に着目します。その共通部分は点*3になっていますが、片方の三角形はそれが「面」（つまり部分集合）になっていません。つまり、2つ目の条件で、右側のようなはり合わせ方を禁止するようなルールを考えているわけです。

（注意）単体的複体は図形ではない

単体的複体は、単体をある「組み合わせ的な性質」で集めた集合であって、図形そのものではありません。そのため、ある特定の次元にある図形という具体的な対象を忘れ、集合の要素と「組み合わせ的性質」の2つだけで抽象化することが可能です。この組み合わせ的性質のみに着目したものを抽象単体複体をいいます。なんてまどろっこしいことを…と感じる方もいると思いますが、このような抽象化がデータをうまく図形に近似して特徴を取り出す際に役に立ちます。詳しい定義などは、参考書を見てください。

図形は単体的複体の幾何学的実現

では、我々が想像する図形はどう定義されているのでしょうか。これは、単体的複体に含まれるすべての要素（つまり、単体）の和集合を取ったものと定義されます。単体的複体を、その要素をとすると

が図形を表します。各単体が我々が住んでいるような空間（ユークリッド空間）にあるときは、をが定める多面体、先程簡単に紹介した抽象単体複体の世界ではをの幾何学的実現といいます。つまり、何かパーツ（単体）を組み合わせた図形を表現する際は、単体的複体を定義するだけでは不十分で、そこから幾何学的実現をとってやる必要がある、ということです。ややこしいですね。

Next Step

この単体的複体は、データという点の集まりからどうやって情報を引き出すか、を考える際に使います。次回はどのような単体的複体がトポロジカルデータアナリシスで使われるかを簡単に紹介していきます。本当は、同値とかホモトピー同値とかいう「同じ」って何ですか、ってところを考えておいたほうがよいのですが、まずは解析するのに必要なキーワード的知識を紹介していくことにフォーカスしていこうと思っています。

*1:おそらくこれが最も素朴な面のイメージだと思います。

*2:回りくどい言い方ですが、単体的複体と図形は異なるものなので、言い分けています。

*3:この言い方は実は正確ではないです。片方の三角形については面の定義を満たしていないからです。一点で交わっているという言い方のほうが良いかもしれません。

複数の棒グラフを表示させるのはpandasが便利

経緯

ある対象に対して、複数のアプローチの結果を可視化したいとき、棒グラフで並べて比較する方法があります。これをmatplotlib.pyplot.bar()で描いていましたが、棒の太さやら目盛の調整が大変でした。matplotlibは柔軟な可視化ができる反面、匠の技が要求されることもあるため、使いこなすのは結構難しいです。より簡単にできる方法はないかと考えていたところ、PandasのDataFrameにmatplotlibのラッパーがあったことを思い出しました。というわけで、復習がてら載せておきます。

可視化する対象

幾つかのモデルの回帰係数の値を可視化します。使用データは、ボストン近郊の住宅情報、モデルは線形回帰、Lasso、Ridge、Elastic Netの4種類です。Lasso、Ridge、Elastic Netは罰則付き回帰と呼ばれる手法で、目的関数に罰則項というモデルの複雑さを調整する項を入れることで、過学習を緩和させることが目的の一つです。手法によって、特定の重要でない変数に対する回帰係数が0になったり、係数が全体的に縮小する、などバリエーションがあります。今回は、代表的な3手法の推定結果にどのような違いが出るかを可視化してみます。

可視化

ライブラリのインポート

必要なライブラリたちをインポートしておきます。matplotlibを呼び出していますが、これはグラフの整形に使っています。

# データ加工
import numpy as np
from pandas import DataFrame
# サンプルデータ
from sklearn.datasets import load_boston
# 学習用、検証用データに分割
from sklearn.cross_validation import train_test_split
# 比較対象：線形回帰
from sklearn.linear_model import LinearRegression
# 罰則付き回帰、正則化項をクロスバリデーションで求めるライブラリを使用
from sklearn.linear_model import LassoCV
from sklearn.linear_model import RidgeCV
from sklearn.linear_model import ElasticNetCV
# MSE、R^2
from sklearn.metrics import mean_squared_error
from sklearn.metrics import r2_score
# 可視化
import matplotlib
import matplotlib.pyplot as plt
%matplotlib inline
# きれいに書くおまじない
plt.style.use('ggplot')
font = {'family': 'meiryo'}
matplotlib.rc('font', **font)

サンプルデータ

サンプルデータは、ボストン近郊の住宅情報です。そろそろ使いすぎてデータマートを作る操作が滑らかになってきました。

# データのロード、データマート作成
boston = load_boston()
df = DataFrame(boston.data, columns = boston.feature_names)
df['MEDV'] = boston.target
# 説明変数、目的変数
X = df.iloc[:, :-1].values
y = df.loc[:, 'MEDV'].values
# 学習用、検証用データに分割
(X_train, X_test, y_train, y_test) = train_test_split(X, y, test_size = 0.3,  random_state = 666)

モデル作成、プロット

matploblibの場合、プロットしたい対象を一つずつ追加していくという方法が取られます。PandasのDataFrameのplotメソッドは、DataFrameのデータを可視化することができます。棒グラフは、plot(kind = 'bar')と指定するかplot.bar()とすると指定できます。DataFrameに対してこのメソッドを使用すると、indexが横軸で、columnsごとにデータを表示します。今回、4つのモデルの係数に関する棒グラフを表示したいので、DataFrameにはindexには各変数を、columnsには各モデルをラベルするようにデータを格納します。

コード例と結果は以下となります。

# モデルインスタンス生成
Linear = LinearRegression()
Lasso = LassoCV()
Ridge = RidgeCV()
ElasticNet = ElasticNetCV()
# 学習
Linear.fit(X_train, y_train)
Lasso.fit(X_train, y_train)
Ridge.fit(X_train, y_train)
ElasticNet.fit(X_train, y_train)
# 予測値
y_train_Linear_pred = Linear.predict(X_train)
y_test_Linear_pred = Linear.predict(X_test)
y_train_Lasso_pred = Lasso.predict(X_train)
y_test_Lasso_pred = Lasso.predict(X_test)
y_train_Ridge_pred = Ridge.predict(X_train)
y_test_Ridge_pred = Ridge.predict(X_test)
y_train_ElasticNet_pred = ElasticNet.predict(X_train)
y_test_ElasticNet_pred = ElasticNet.predict(X_test)
# データフレームに格納
df_coef = DataFrame(data = np.transpose(np.array([Linear.coef_, Lasso.coef_, Ridge.coef_, ElasticNet.coef_])),  
                    index = df.iloc[:, :-1].columns, 
                    columns = ['Linear MSE train: %.3f, test: %.3f, R^2 train: %.3f, test: %.3f' 
                               %(mean_squared_error(y_train, y_train_Linear_pred),  mean_squared_error(y_test, y_test_Linear_pred),
                                 r2_score(y_train, y_train_Linear_pred), r2_score(y_test, y_test_Linear_pred)),
                               'Lasso MSE train: %.3f, test: %.3f, R^2 train: %.3f, test: %.3f' 
                               %(mean_squared_error(y_train, y_train_Lasso_pred),  mean_squared_error(y_test, y_test_Lasso_pred),
                                 r2_score(y_train, y_train_Lasso_pred), r2_score(y_test, y_test_Lasso_pred)), 
                               'Ridge MSE train: %.3f, test: %.3f, R^2 train: %.3f, test: %.3f' 
                               %(mean_squared_error(y_train, y_train_Ridge_pred),  mean_squared_error(y_test, y_test_Ridge_pred),
                                 r2_score(y_train, y_train_Ridge_pred), r2_score(y_test, y_test_Ridge_pred)), 
                               'ElasticNet MSE train: %.3f, test: %.3f, R^2 train: %.3f, test: %.3f' 
                               %(mean_squared_error(y_train, y_train_ElasticNet_pred),  mean_squared_error(y_test, y_test_ElasticNet_pred),
                                 r2_score(y_train, y_train_ElasticNet_pred), r2_score(y_test, y_test_ElasticNet_pred)) ])
# プロット
df_coef.plot.bar(figsize = (15, 8), width = 1)
plt.title('罰則付き回帰で推定された回帰係数の比較', size = 20)
plt.xlabel('特徴量')
plt.ylabel('回帰係数の大きさ')
plt.legend(loc = 'best')

いい感じに横に並べることができました。ちなみに、棒の太さはwidthでいじることができ、デフォルトは0.8です。これを大きくしすぎると、棒グラフどうしがかぶってしまって何のことやら…という状態になるので、振る舞いを見てちょうどよい値を選択すると良いと思います。

回帰係数の結果ですが、Lassoは重要でない変数が0になる、Ridgeは線形回帰と少し似ているが、全体的に係数のサイズが小さめ、Elastic Netはそれらの折衷案、のような振る舞いをしています。また、MSEとを見ると、LassoとElastic Netは学習データの精度を少し犠牲にして汎化性能が向上しているようには見えます。ただし、全体的に線形回帰と比較して精度が落ちているため、これらの間で何を採用するかは悩ましいところです。

他のプロット

pandasのplotだけでも工夫するといろいろといじれます。例えば、以下が参考になると思います。

sinhrks.hatenablog.com

matplotlibを使いこなして職人になるのもよいですが、時間の制約があった場合はこういった便利なラッパーを使うとよい場面があるかもしれません*1。

*1:自分は割りと原理主義者なのでmatplotlibを調整できるようになりたいですが（笑）。

アソシエーション分析のちょい復習

動機

最近こんな本を購入しました。
www.shoeisha.co.jp
第I部の「データアナリティクスの基礎」は、データ分析のプロジェクトの全体像に触れていて、データ分析のプロジェクトをどう進めていくかについて整理するのには良いのかな、と思いました*1。後半部分の分析内容は結構あっさりしていて、すでに専門的なスキルを持っている方向けではなさそうです。

この書籍の第8章で、Rを使ったアソシエーション分析を取り扱っています。これは、マーケティングの分野でよく使用される手法です。顧客の購買履歴から、「商品Aを買った人は商品Bを買う確率が高い」といった関係性を見つけ出します。自分は、業務では予測分析がメインなので、そういえば久しくこういった手法に触れてなかったな…と感じました。というわけで、復習の一環で簡単にデモしてみようと思います*2。

アソシエーション分析の説明は、例えば以下のサイトが参考になるかと思います。
sinhrks.hatenablog.com
sinhrks.hatenablog.com
実務では、今回使用するサンプルデータのような形式に簡単にまとめられない場合があるので、前処理も含めた解説は非常にためになります。

本当はPythonでやろうと思ったけど、自分の環境ではうまく構築できなかった

Pythonでは、Orangeというものでアソシエーション分析が可能ということで、pipでOrange（何も指定しないとVer. 2.7.8）をインストールしようとしました。しかし、setup.pyのビルドでコケてしまい、結局解消せず*3。大本のページ（Orange Data Mining）を見ると、最新版はVer. 3.3.7とのこと。pip経由だと、Orange3と指定するとインストールできそうなので実行。こちらはインストールに成功しましたが、アソシエーション分析に必要なメソッドが無い。OrangeでPythonからAprioriを動かす - Qiitaを見ると、どうやら本当に無い模様。

というわけで、いったんPythonで進めるのは保留にして、Rでデモをやるだけにします。

Rの{arules} + {arulesVIz}を使って分析

Rでアソシエーション分析を行う場合、例えば{arules}というパッケージがあります。今回は、このパッケージを使用していきます。

{arules}のインストール

if(!require(arules)){
  install.packages("arules", dependencies = TRUE)
}
require(arules)

データセットの読み込み

店舗の30日間のPOSデータを使用します。9835トランザクション、196カテゴリからなるデータで、「milk（牛乳）」、「root vegitables（根菜類）」といった商品の併売情報が入っています。

data(Groceries)
Groceries

# 出力
transactions in sparse format with
 9835 transactions (rows) and
 169 items (columns)

先頭の数行を確認してみます。併売情報を出力するのに、inspect()を用います。

inspect(head(Groceries))

# 出力
  items                     
1 {citrus fruit,            
   semi-finished bread,     
   margarine,               
   ready soups}             
2 {tropical fruit,          
   yogurt,                  
   coffee}                  
3 {whole milk}              
4 {pip fruit,               
   yogurt,                  
   cream cheese ,           
   meat spreads}            
5 {other vegetables,        
   whole milk,              
   condensed milk,          
   long life bakery product}
6 {whole milk,              
   butter,                  
   yogurt,                  
   rice,                    
   abrasive cleaner}

このようなデータを作ったり、アクセスしたりするのは一筋縄ではいかないのですが、ここではうまくトランザクションデータが手に入っているという状況で先に進みます。

summaryで基本統計量の確認

summary(Groceries)

# 出力
transactions as itemMatrix in sparse format with
 9835 rows (elements/itemsets/transactions) and
 169 columns (items) and a density of 0.02609146 

most frequent items:
      whole milk other vegetables       rolls/buns             soda           yogurt 
            2513             1903             1809             1715             1372 
         (Other) 
           34055 

element (itemset/transaction) length distribution:
sizes
   1    2    3    4    5    6    7    8    9   10   11   12   13   14   15   16   17 
2159 1643 1299 1005  855  645  545  438  350  246  182  117   78   77   55   46   29 
  18   19   20   21   22   23   24   26   27   28   29   32 
  14   14    9   11    4    6    1    1    1    1    3    1 

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  1.000   2.000   3.000   4.409   6.000  32.000 

includes extended item information - examples:
       labels  level2           level1
1 frankfurter sausage meat and sausage
2     sausage sausage meat and sausage
3  liver loaf sausage meat and sausage

出力結果の意味は以下です。

• most frequent items : 頻度（=商品が売れた数）
• element (itemset/transaction) length distribution : 1トランザクションで購入した商品の数
• Min, 1st Qu., Median, Mean, 3rd Qu., Max : 最小値、第一四分位数、中央値、平均値、第三四分位数、最大値

商品ごとのヒストグラム

先程の要約だと、商品ごとの頻度の情報が十分ではありませんでした。ここでは、購入された商品に対して相対出現頻度によるヒストグラムを出力します。これは、itemFrequencyPlot()で実行できます。全商品を可視化すると複雑になるので、引数のsupportで支持度4%以上の商品のみに限定しています。

# 目盛りを無駄にメイリオに設定
quartzFonts(Hir = c("HiraMaruPro-W4", "HiraMaruProN-W4", 
                    "Meiryo", "HiraMinProN-W3"))
itemFrequencyPlot(Groceries, support = 0.04, cex.names = 0.7, 
                  col = "red", family = "Hir", font = 3)

アソシエーション分析の実行

apriori()で実行できます。support（支持度）とconfidence（確信度）は、それぞれ手で0.5%以上、1%以上をしきい値としました。

rules <- apriori(Groceries, parameter = list(support = 0.005, confidence = 0.01))

# 表示結果
Apriori

Parameter specification:
 confidence minval smax arem  aval originalSupport support minlen maxlen target   ext
       0.01    0.1    1 none FALSE            TRUE   0.005      1     10  rules FALSE

Algorithmic control:
 filter tree heap memopt load sort verbose
    0.1 TRUE TRUE  FALSE TRUE    2    TRUE

Absolute minimum support count: 49 

set item appearances ...[0 item(s)] done [0.00s].
set transactions ...[169 item(s), 9835 transaction(s)] done [0.00s].
sorting and recoding items ... [120 item(s)] done [0.00s].
creating transaction tree ... done [0.00s].
checking subsets of size 1 2 3 4 done [0.00s].
writing ... [2138 rule(s)] done [0.00s].
creating S4 object  ... done [0.00s].

この条件だと併売に関する2138個のルールが出来上がります。

特定の商品に注目した場合

アソシエーションのルールは、「X ならば Y」という形になっています。作成したルールを用いて、Yを固定したときにXがどのようなものがあるかを確認してみます。

# Y = "onions"（たまねぎ）にしてみる
onionRules <- subset(rules, subset = rhs %in% "onions")
# リフト値が高い順に表示
inspect(sort(onionRules, by = "lift"))

# 出力
     lhs                                   rhs      support     confidence lift     
1304 {root vegetables,other vegetables} => {onions} 0.005693950 0.12017167 3.8750440
1307 {other vegetables,whole milk}      => {onions} 0.006609049 0.08831522 2.8478038
316  {root vegetables}                  => {onions} 0.009456024 0.08675373 2.7974523
324  {other vegetables}                 => {onions} 0.014234875 0.07356805 2.3722681
308  {whipped/sour cream}               => {onions} 0.005083884 0.07092199 2.2869434
310  {citrus fruit}                     => {onions} 0.005592272 0.06756757 2.1787771
314  {tropical fruit}                   => {onions} 0.005693950 0.05426357 1.7497776
312  {bottled water}                    => {onions} 0.005897306 0.05335787 1.7205725
320  {yogurt}                           => {onions} 0.007219115 0.05174927 1.6687019
326  {whole milk}                       => {onions} 0.012099644 0.04735376 1.5269647
322  {rolls/buns}                       => {onions} 0.006812405 0.03703704 1.1942927
47   {}                                 => {onions} 0.031011693 0.03101169 1.0000000
318  {soda}                             => {onions} 0.005287239 0.03032070 0.9777183

これより「root vegetables（根菜類）, other vegetables（他の野菜）」やroot vegetables（根菜類）, whole vegetables（牛乳）」といった条件がonions（たまねぎ）に結びつくものがリフト値が高い（単に玉ねぎを買うよりよりもXを買ってさらにYを買う確率が高い）ことがわかります。

結果の可視化

{arulesViz}のインストール

if(!require(arulesViz)){
  install.packages("arulesViz", dependencies = TRUE)
}
require(arulesViz)

（図が単純になるように）ルールを再設定

rules <- apriori(Groceries, parameter = list(support = 0.004, 
                                             confidence = 0.5))

以下では、有用そうな図をいくつか紹介していきます。

散布図

散布図は、method = "scatterplot"で指定します。

graphics::plot(rules, method = "scatterplot")

横軸が支持度、縦軸が確信度で、点の色が赤に近いほどリフト値が高くなっています。これより、例えばリフト値が高い領域がわかり、パラメータ設定の参考にすることができます。

バブルチャート

method = "grouped"と指定することで、XとYの関係を直観的（丸の大きさと色）で表現できます。

graphics::plot(rules, method = "grouped")

横軸がアソシエーションの左辺（X）、縦軸が右辺（Y）、支持度の大きさは丸のサイズ、リフト値の大きさは色の濃さで表されています。

アソシエーショングラフ

有向グラフ形式（矢印で関係が表現されているグラフ）で表現することも可能です。表示を簡単にするため、リフト値の上位10項目のみに着目します。method = "graph"で有向グラフを表現できます。

# リフト値の高い順に10項目
rules_high_lift <- head(sort(rules, by = "lift"), 10)
graphics::plot(rules_high_lift, method = "graph", control = list(type = "items"))

herbsを購入した人は、root vegetablesを購入する可能性が高い（?）などが表示されます。
実はWindowとMacでこの表示が異なるそうで、Windowsはマウス操作で移動ができるそうです。

おわりに

アソシエーション分析を用いると、複雑なトランザクションデータだけを見ている場合にはわからなかったインサイトが得られる可能性があります。しかしながら、ルールがわかってもそれがビジネス戦略とあっているか*4はまた別の問題です*5。アソシエーション分析で課題解決を図る際は、きちんと目的を明確にし、分析が目的に沿っているかを精査しておく必要があります。自分が業務で携わるときは、このあたりを肝に命じて進めていこうと思います。

*1:一応自分はコンサルタントなので、たびたび全体感を整理しとくのは業務でも有用かなと思っています。

*2:というわけで非常にクオリティが低い内容です。。

*3:自分はPython3.5を使っているのが原因の一つっぽいですが、setup.pyをPython3に合うように書き換えても成功には至らず。。

*4:売上を上げたいなど。

*5:ルールが利益に結びつかない可能性があります。

FastBDTの計算時間が速いかを確認してみる

はじめに

最近、ブースティング系のアルゴリズムでXGboostより速いものが実装されているようです。
github.com

論文は以下になります。
[1609.06119] FastBDT: A speed-optimized and cache-friendly implementation of stochastic gradient-boosted decision trees for multivariate classification

論文を見ると、Stochastic Gradient Boosted Decision Tree (SGBDT)という手法が用いられているようです*1。用途としては、信号（例えば、ラベル1）と背景事象（1以外）を分類する場合が想定されています。このアルゴリズムは、素粒子実験の一つである、Belle II実験の解析で使用されているとのことです。ばくっと探していたらカールスルーエの学生の修論があったので、載せておきます。
https://ekp-invenio.physik.uni-karlsruhe.de/record/48602/files/EKP-2015-00001.pdf

Belle II実験

Belle II実験に関しては、修士課程のときの自分の研究内容と深く関わっているので少し補足を。
Belle II実験とは、日本のKEKにある電子-陽電子加速器を用いた実験です。この実験の目的の一つ*2は、B中間子とよばれる、ボトムクォークを含む複合粒子をたくさん生成し、その「稀な」崩壊現象の精密測定から新しい理論の探索を行うことです*3。標準模型を超える物理があった場合、直接高エネルギー加速器で直接生成されるだけでなく、間接的に様々な物理量に寄与する可能性があります*4。その寄与は、標準模型の予測値からも「ずれ」という形で表れます。その「ずれ」を実験で精密に検証することで、新しい物理の姿を間接的に見ることができます。

間接的に新しい物理の効果を見るにあたって、「稀な」崩壊現象を見ることが非常に大事です。それは、比較的起こりやすい崩壊現象よりも、新しい物理の効果を捉えやすいためです。標準模型では、例えばGlashow‐Iliopoulos‐Maiani機構（GIM機構）といった、ある種類の崩壊を抑制する機構が備わっているため、素朴な次元解析での評価よりも崩壊率が小さいものがあります。これらは標準模型では「稀な」崩壊と予測されます。しかし、新しい物理は必ずしもGIM機構のようなものが備わっているわけではないため*5、大きな補正が「稀な」崩壊現象に加わる可能性があります。そのため、比較的起こりやすい崩壊現象とくらべて変化を捉えやすくなります。

B中間子に着目する理由は、その「稀さ」が絶妙だからです。B中間子の物理現象は、抑制されているカテゴリーの中では比較的大きな寄与を持っています*6。そのため、稀ではあるが、標準模型が正しいと仮定した場合の予測値が実験で検知できる程度には寄与がある、といった現象になり、「標準模型からのずれ」の検証が他の中間子よりも容易になっています。

インストールから使用可能になるまで

つい素粒子の話になると長くなってしまうので、話を戻します。Mac OSの場合、サイトから直接ダウンロードし、makeすれば使えるようになります。cmakeは予めインストールしておきます。

$ git clone https://github.com/thomaskeck/FastBDT
$ cd FastBDT/
$ cmake . 
$ make
$ make install

上記の操作をすれば、実行ファイルが生成されるので、FastBDTのディレクトリで

$ ./FastBDF [あといろいろ必要なオプション指定]

と唱えると実行できます。

Pythonから使用したい場合

Pythonで動かす方法もあります。ただし、まだ完全に設定などを理解していないので、現状を備忘録的にまとめておきます。

修正点

まず、インポートする元は ~/FastBDT/python/FastBDT.pyにあります。ただし、このスクリプトで指定しているライブラリの名前が自分の環境では異なっていたので、少し修正が必要です。修正方法は2通りあります。

• [19行目]で指定しているライブラリをlibFastBDT_CInterface.soからlibFastBDT_CInterface.dylibに書き換える
• [14行目-16行目]のコメントアウトを外し、[19行目]をコメントアウト

こうすれば、インポートしたときに正しく読み込みができるようになりました。

インポート

以下のようにパスを上手く指定してインポートします。

import sys
sys.path.append('/Users/[ユーザ名]/FastBDT/python')
import FastBDT

インポートまでの一連の流れですが、不慣れなこともあり、かなり洗練されていない感じになっています（というよりセオリーを知らない。）が、ご了承ください（汗）。

モデルはClassifierで作成することができます。モデルインスタンス作成後は、scikit-learnのようにfitやpredictといったメソッドが使えます。

計算時間の比較

学習にかかった時間を計測します。比較対象として、パーセプトロン、ランダムフォレスト、XGboostを選択します。

サンプルデータ

FastBDTは'1'とそれ以外を識別する問題にのみ対応しているようです*7。そのため、簡単のため、人工的に2値データを生成して、それで比較することにしました。クラス分類のテスト用の人工データは、scikit-learn.datasetのmake_classificationで作成することが可能です。参考にしたサイトは以下です。
overlap.hatenablog.jp

説明変数が5つ（n_features = 5）、サンプル数100万（n_samples = 1000000）の2値データ（n_classes = 2）を生成します。他のオプションは適当です。

from sklearn.datasets import make_classification
# データ作成
dat = make_classification(
    n_samples=1000000, n_features=5, n_informative=2, n_redundant=2,
    n_repeated=0, n_classes=2, n_clusters_per_class=2, weights=None,
    flip_y=0.01, class_sep=1.0, hypercube=True, shift=0.0, scale=1.0,
    shuffle=True, random_state=None)

データマートを作成し、さらに学習用、検証用データを作成します。

# ライブラリ
import pandas as pd
from pandas import DataFrame
# DataFrameに格納
df = DataFrame(dat[0], columns = ['Var1', 'Var2', 'Var3', 'Var4', 'Var5'])
df['Class'] = dat[1]

# 説明変数、目的変数
X = df.iloc[:, :-1].values
y = df.loc[:, 'Class'].values
# 学習用、検証用データ作成
from sklearn.cross_validation import train_test_split
(X_train, X_test, y_train, y_test) = train_test_split(X, y, test_size = 0.3, random_state = 666)

計算時間、誤分類数、Accuracyを算出

ばくっと計算する関数を作っておきます。

def CalcTimes(mod, X_train, X_test, y_train, y_test):
    import numpy as np
    from sklearn.metrics import accuracy_score
    
    start = time.time()
    # モデルインスタンス作成、学習
    mod = mod
    mod.fit(X_train, y_train)
    elapsed_time = time.time() - start
    print('elapsed_time: %.2f' %elapsed_time + ' [sec]')
    
    # 誤分類数、Accuracy
    y_train_pred = np.round(mod.predict(X_train))
    y_test_pred = np.round(mod.predict(X_test))
    print('Misclassified samples: %d' % (y_test != np.round(y_test_pred)).sum())
    print('Accuracy: %.2f' % accuracy_score(y_test, y_test_pred))

FastBDT

from FastBDT import Classifier
mod = Classifier()
CalcTimes(mod, X_train, X_test, y_train, y_test)

# 出力
elapsed_time: 6.50 [sec]
Misclassified samples: 34948
Accuracy: 0.88

パーセプトロン

from sklearn.linear_model import Perceptron
mod = Perceptron()
CalcTimes(mod, X_train, X_test, y_train, y_test)

# 出力
elapsed_time: 0.56 [sec]
Misclassified samples: 66686
Accuracy: 0.78

ランダムフォレスト

from sklearn.ensemble import RandomForestClassifier
mod = RandomForestClassifier()
CalcTimes(mod, X_train, X_test, y_train, y_test)

# 出力
elapsed_time: 44.74 [sec]
Misclassified samples: 39241
Accuracy: 0.87

XGboost

from xgboost import XGBClassifier
mod = XGBClassifier()
CalcTimes(mod, X_train, X_test, y_train, y_test)

# 出力
elapsed_time: 45.03 [sec]
Misclassified samples: 34131
Accuracy: 0.89

比較結果

今回のデータセットに関しては、Macbook Air使用で、FastBDTはランダムフォレストやXGboostと比較して6倍以上速い結果になりました。しかも、Accuracyが大きく損なわれずにです。条件を変えると倍率は変わると思われますが、概ね似たカテゴリーのアルゴリズムよりは高速なんだと思います。もっと開発が進んで使いやすくなったら、非常に強力なアルゴリズムなのではないかと期待されます。

*1:細かい部分で工夫があるようですが、まだキャッチアップできていません。

*2:他にもK中間子の稀崩壊探索なども目的です。

*3:前身として、Belle実験というものが同じKEKで行われていました。Belle実験では、B中間子の崩壊現象を用いた標準模型の検証が主目的でした。

*4:量子効果（摂動論の高次効果）で寄与します。

*5:4次元以外の「余剰な」次元を持っている理論がその例です。超対称性を持った理論でSuper GIM機構というものがある場合は、余剰次元をもつ理論よりは寄与が抑制されます。

*6:色々省略しますが、これはトップクォークという素粒子がとてつもなく重いことが起源です。

*7:使ってみた結果、こう認識しています。間違いがあればご指摘ください。