読者です 読者をやめる 読者になる 読者になる

データサイエンティスト(仮)

元素粒子理論屋。いっぱしのデータ分析屋になるために修行中です。

久々にC++を触ってみたら自明な入力しかできないのはくやしい

C++

自明な振り返り

自明な出力しかできませんでした。
tekenuko.hatenablog.com

動機

くやしいので非自明なことを言いたい。

自明(?)なコード

自明なことを入力すると自明だと言われます。

// trivial.cpp

# include <iostream>
# include <string>
using namespace std;

int main(){
  string str;
  cout << "自明な言葉はなーんだ?" << endl;

  // 空白があっても一つのカタマリと認識させるのに使用
  getline(cin, str);

  if(str == "Hello world!"){
    cout << "自明ですね" << endl;
  } else{
    cout << "自明じゃないですね" << endl;
  }

  return 0;
}

自明な実行

trivialという実行ファイルを作って実行します。自明なクイズが始まるので、自明な入力をすると...

$ g++ -o trivial trivial.cpp 
$ ./trivial

# 出力
自明な言葉はなーんだ?  //自明なクイズが始まる
Hello world! // コマンドラインで自明に入力する
自明ですね //自明だったらしいです

非自明なことをいう

とりあえず非自明なことを言ってみます。

$ g++ -o trivial trivial.cpp 
$ ./trivial

# 出力
自明な言葉はなーんだ?  //自明なクイズが始まる
ほとんどいたるところ連続 // コマンドラインで入力する
自明じゃないですね //自明じゃなかったらしいです

非自明だったようです。よかったです。

久々にC++を触ってみたら自明なことしかできない

C++

自明な背景

3年ぶりくらいにC++触ってみたら自明なことしかできなかったので、自明なことをしたメモをします。

自明なコード

実行すると自明な言葉を発します。

// helloworld.cpp

#include <iostream>
using namespace std;

int main(){
  cout << "Hello world!" << endl;
  return 0;
}

自明な実行

自分はコマンドライン上でコンパイルして実行するという自明なことしかできない。helloworidという実行ファイルを使って実行するという自明な営みをすると自明なことをします。

$ g++ -o helloworld helloworld.cpp 
$ ./helloworld 

# 出力
Hello world!

自明なNext Step

今の自分のC++力は自明なコードを書けるレベルとホモトピー同値です。
修行します。

Pythonでデータ分析:機械学習の自動化

Python auto-sklearn

導入

何か問題を解決するにあたって機械学習を活用する場合、膨大なアルゴリズム、そのアルゴリズムに付随する多くのハイパーパラメータが存在します。分析の要件が「とにかく精度、中身は問わない」だった場合、何とかして効率的にモデルとパラメータを知りたい、という状況が起こりえます*1

そのような要件をサポートしてくれるPythonのライブラリに、auto-sklearnというものがあります。
github.com

これは、大雑把にいうと、与えられたデータを見て、scikit-learnの中から

  • 良さそうな前処理の選定
  • 良さそうなアルゴリズムの選定
  • 良さそうなハイパーパラメータの選定

を行い、これらのアンサンブルを取る、といったことをします。

f:id:tekenuko:20161002143359p:plain
まず、Meta Learningという、データの数、種類などから、どの手法がどのデータに向いているかをモデル化し、ベイズ的最適化により最適な組み合わせを探る、そしてアンサンブルをとる、という流れです。より詳しくは、以下の論文をご参照ください。
What is auto-sklearn? — AutoSklearn 0.0.2 documentation

インストール

pip経由でインストールします。まずは、auto-sklearnを使うために必要なライブラリをインストール後、auto-sklearnをインストールします。

pip install -r https://raw.githubusercontent.com/automl/auto-sklearn/master/requirements.txt
pip install auto-sklearn

注意

auto-sklearnはUbuntuで開発されたらしく、Linuxでは動作するがMacやWindowでは動かないかもしれない、とのこと。実際、自分の環境(MacOS)で最初はうまくいきませんでした。自分の場合は、gccのバージョンが古い(4.8以前)ことが原因だったようで、バージョンを4.8以降にしたら導入できました。

コードを試したいときは

例えば、auto-sklearnのページにコード例があります。これをJupyter notebookか何かにコピーして動かしてみると良いかと思います*2
What is auto-sklearn? — AutoSklearn 0.0.2 documentation

注意

やんちゃに動かすと大量のログを吐きながら膨大な計算をします*3。少ないデータ量でも2〜3時間はかかります。オプションでモデル探索の時間などをいじれるので、すぐに結果を見たいならオプションの指定を適切に行ったほうが良いでしょう。
APIs — AutoSklearn 0.0.2 documentation

所感

ローカルPCで動かすには、やっていることと処理時間の両方の観点でオーバースペック気味かなという感想です。使いどころとしては、ハイスペックな分析環境があり、かつどうしても精度がほしいといった場合かなと思います。

*1:個人としてはこういった中身の理解が介在しないことはしたくないですが、そういう要件の場合もありますので仕方ない。

*2:何かデモを載せようかと思いましたが、注意に書いた状況のためやめにしました(汗)。

*3:今もPCがものすごいうなりながら計算してます(汗)。

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

トポロジカルデータアナリシス 不定期

導入

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

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

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

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

問題点

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

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

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

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

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

チェック複体

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

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

f:id:tekenuko:20160928223420j:plain

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

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

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

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

アルファ複体

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

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

f:id:tekenuko:20160928225105j:plain

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

これまで、TDAで登場する単体的複体についてざっと紹介しました。しかし、なぜこれらの単体的複体で近似すればよいかについては何も言っていませんでした。実は、点の周りに描く領域を凸集合に(球は凸集合)し、k+1個の球が重なっているときにk単体が要素になる(チェック複体とアルファ複体)ような抽象単体複体は、もとの見立てた図と定性的な性質が同じである*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次元三角形、...と呼ぶことにします。これは、N次元空間内のk次元三角形という拡張ができます。このk次元三角形をk単体といいます。もう少し正確には、N次元空間内にk + 1個の点を配置し、線形独立なベクトルを与え、その凸包を考える、という定式化をしますが、詳しくは参考書を参考にしていただけばと思います。

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

f:id:tekenuko:20160926221730p:plain

単体的複体(複体)

ゆるゆる定義

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

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

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

http://infoshako.sk.tsukuba.ac.jp/~hachi/math/csc2/image/complex.gif
参考(http://infoshako.sk.tsukuba.ac.jp/~hachi/math/csc2/complex.html

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

(注意)単体的複体は図形ではない

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

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

では、我々が想像する図形はどう定義されているのでしょうか。これは、単体的複体に含まれるすべての要素(つまり、単体)の和集合を取ったものと定義されます。単体的複体をK、その要素を\sigmaとすると
|K| = \cup _{\sigma \in K}\sigma
が図形を表します。各単体が我々が住んでいるような空間(ユークリッド空間)にあるときは、|K|Kが定める多面体、先程簡単に紹介した抽象単体複体の世界では|K|K幾何学的実現といいます。つまり、何かパーツ(単体)を組み合わせた図形を表現する際は、単体的複体を定義するだけでは不十分で、そこから幾何学的実現をとってやる必要がある、ということです。ややこしいですね。

Next Step

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

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

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

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

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

Python 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')

f:id:tekenuko:20160924214300p:plain

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

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

他のプロット

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

sinhrks.hatenablog.com


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

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

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

R アソシエーション分析

動機

最近こんな本を購入しました。
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)

f:id:tekenuko:20160924005258p:plain

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

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")

f:id:tekenuko:20160924011554p:plain

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

バブルチャート

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

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

f:id:tekenuko:20160924011607p:plain

横軸がアソシエーションの左辺(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"))

f:id:tekenuko:20160924012017p:plain

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

おわりに

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

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

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

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

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

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