Pandasで「ビジネス活用事例で学ぶデータサイエンス入門」を勉強する（第５章）

はじめに

マーケティング寄りのデータ分析の知識を補うため、以下の本で勉強を始めた。事例ベースな内容で、とても読みやすい。 Pandasも習得したいので、Pandasに翻訳しながら読み進めている。今回は第５章を勉強した。

ビジネス活用事例で学ぶ データサイエンス入門

ちなみに、CSVファイルなどのデータは、本のサポートページ（SBクリエイティブ:ビジネス活用事例で学ぶ データサイエンス入門）で配布されている。

翻訳したコード

第５章では、Web系のデータ分析では定番？のバナー広告のA/Bテストを扱っている。作業はJupyter Notebook上で行った。 カイ二乗検定が使われているが、Pandasにはカイ二乗検定はない様子なので、SciPyのstats.chi2_contingencyを使う。

%matplotlib inline # Jupyter Notebookでmatplotlibの図が表示されないときにつけるおまじない
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
from scipy.stats import chi2_contingency

# データを読み込む。
ab_test_imp = pd.read_csv('section5-ab_test_imp.csv')
ab_test_goal = pd.read_csv('section5-ab_test_goal.csv')

# transition_idが同じものをくっつける。
ab_test_imp_goal = pd.merge(ab_test_imp, ab_test_goal, on=['transaction_id'], how='outer', suffixes=['','.g'])

# user_id.gがNaNだった場合に「0」それ以外では「1」として、クリックされたかを判定するフラグとする。
ab_test_imp_goal['is_goal'] = 1 - ab_test_imp_goal['user_id.g'].apply(np.isnan).astype(int)

# 集計してクリック率を計算する。
sum_goal = ab_test_imp_goal.groupby('test_case')['is_goal'].sum()
sz_user = ab_test_imp.groupby('test_case')['user_id'].size()
sum_goal / sz_user
# ※ Jupyter Notebook上には以下の用にクリック率が表示される。
#    Aが約8%で、Bが約11%となった。
# test_case
# A    0.080256
# B    0.115460
# dtype: float64

# カイ二乗検定により、クリック率に統計的な差があるかを確認する。
cr = ab_test_imp_goal.pivot_table(index='test_case',columns='is_goal', values='user_id', aggfunc='count')
# ※ crは以下のようなクロス集計結果となる
# is_goal    0        1
# test_case
# A          40592    3542
# B          38734    5056
#
chi2, p, dof, expected = chi2_contingency(cr.as_matrix())
print(chi2)
# 308.375052893
print(p)
# 4.93413963379e-69

chi2の値が本と同様の値となり、p値も本と同様にとても小さな値となった。 p値が小さいということで 「バナーAとバナーBでは、クリックされなかった（0）とクリックされた（1）の割合に差がない」という帰無仮説が棄却され、広告をバナーAとバナーBに分けたことで、クリック率が変わったと判断できる。 本では、この後、テストケースごとのクリック率を計算したりするが、今回は省略。

おわりに

Rのchisq.test関数と勝手が違ったので、手間取ってしまった。本ではカイ二乗検定の詳しい説明がないので、気になる方は別な統計本を買って勉強する必要がある。「統計的に有意な差」と「ビジネス的に意味があるか」という２つの観点が出てきたりして、ちょっと、この章は読みづらい感がある。

Pandasで「ビジネス活用事例で学ぶデータサイエンス入門」を勉強する（第４章）

はじめに

マーケティング寄りのデータ分析の知識を補うため、以下の本で勉強を始めた。事例ベースな内容で、とても読みやすい。 Pandasも習得したいので、Pandasに翻訳しながら読み進めている。今回は第４章を勉強した。

ビジネス活用事例で学ぶ データサイエンス入門

ちなみに、CSVファイルなどのデータは、本のサポートページ（SBクリエイティブ:ビジネス活用事例で学ぶ データサイエンス入門）で配布されている。

翻訳したコード

第４章では、クロス集計を扱っている。作業はJupyter Notebook上で行った。Pandasなコードは以下のとおり。

%matplotlib inline # Jupyter Notebookでmatplotlibの図が表示されないときにつけるおまじない
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np

# データを読み込む。
dau = pd.read_csv('section4-dau.csv')
user_info = pd.read_csv('section4-user_info.csv')

# user_idとapp_nameが同じものをくっつける。
dau_user_info = pd.merge(dau, user_info, on=['user_id', 'app_name'])

# 月次で集計するため、日付を月までにする。
dau_user_info['log_month'] = dau_user_info['log_date'].str.slice(0, 7)

# 性別で集計する。
dau_user_info.pivot_table(index='log_month',columns='gender', values='user_id', aggfunc='count')

gender	F	M
log_month
2013-08	47343	46842
2013-09	38027	38148

# 年代で集計する（表は省略）。
dau_user_info.pivot_table(index='log_month',columns='generation', values='user_id', aggfunc='count')

# 性別と年代で集計する（表は省略）。
dau_user_info.pivot_table(index='log_month',columns=['gender', 'generation'], values='user_id', aggfunc='count')

# デバイスで集計する（表は省略）。
dau_user_info.pivot_table(index='log_month',columns='device_type', values='user_id', aggfunc='count')

# 集計結果を可視化する。
dau_user_info.groupby(['log_date', 'device_type'])['app_name'].count().groupby('device_type').plot(legend=True,rot=45)

デバイスでわけて、日付ごとのユーザー数を集計し、可視化したものが以下のとおり。 本と同様に、9月になってから、Androidユーザーが少なくなっている。

おわりに

pivot_tableで多重なクロス集計もできるので便利。今回も手こずったのは、plotの部分で、x軸のラベルを日付だけにしたかったけど諦め。細かいことは、matplotlib使えば良いんだろうけど、勉強中なので、しばらくはPandasのplotにこだわりたい。

Pandasで「ビジネス活用事例で学ぶデータサイエンス入門」を勉強する（第３章）

はじめに

マーケティング寄りのデータ分析の知識を補うため、勉強を開始した。「チュートリアル的な事例ベースの教材がないかな〜」と色々と探していたところ、ぴったりの良い本が見つかった。第１章と第２章には、データ分析がどういう仕事か書かれている。事例は第３章からとなる。

ビジネス活用事例で学ぶ データサイエンス入門

この本のみならず、多くの（マーケティング寄りの）データ分析の本では、R言語が使われている。PythonのPandasも習得したかったので、Pandasで翻訳しながら読み進めることにした。ちなみに、CSVファイルなどのデータは、本のサポートページ（SBクリエイティブ:ビジネス活用事例で学ぶ データサイエンス入門）で配布されている。

翻訳したコード

第３章では、売上を可視化・比較するため、ヒストグラムを作成する。作業はJupyter Notebook上で行った。Pandasのコードは以下のとおり。

%matplotlib inline # Jupyter Notebookでmatplotlibの図が表示されないときにつけるおまじない
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np

# データを読み込む。
dau = pd.read_csv('section3-dau.csv')
dpu = pd.read_csv('section3-dpu.csv')
install = pd.read_csv('section3-install.csv')

# DAUデータとInstallデータで、user_idとapp_nameが同じものをくっつける。
dau_install = pd.merge(dau, install, on=['user_id', 'app_name'])

# さらにDPUデータをくっつける。
dau_install_payment = pd.merge(dau_install, dpu, on=['log_date', 'app_name', 'user_id'], how='outer')

# 未課金ユーザの課金額をNAから0とする。
dau_install_payment = dau_install_payment.fillna(0)

# 月次で集計するため、日付を月までにする。
dau_install_payment['log_month'] = dau_install_payment['log_date'].str.slice(0, 7)
dau_install_payment['install_month'] = dau_install_payment['install_date'].str.slice(0, 7)

# 利用月・利用開始月・ユーザーIDでまとめる。
mau_payment = dau_install_payment.groupby(['log_month', 'user_id', 'install_month']).sum().reset_index()

# 利用月と利用開始月を確認することで、新規ユーザーか既存ユーザーかを判定する。
mau_payment['user_type'] = np.where(mau_payment['install_month']==mau_payment['log_month'], 'install', 'existing')

# 課金額を集計して可視化する。
mau_payment_summery = mau_payment.groupby(['log_month', 'user_type']).sum().reset_index()
mau_payment_summery.pivot('log_month', 'user_type').plot.bar(y='payment', stacked=True)

可視化した先月と今月の売上状況は以下のとおり。本では、この後、新規ユーザーの売上比較も行うが、今回は省略。

おわりに

Numpy/Scipy脳の私には、Pandasは、ちょっととっつきにくい。一番手こずったのは、plotの部分で、なかなか思ったような図が出せなかった。このへんは「慣れ」な気もしなくもない。逆にR言語に慣れてる人は、関数名の違いさえ把握してしまえば、Pandasを使いこなすのは余裕なのでは？

PyCallを使えばRubyでもKerasでDeep Learningができる

mrknさんが開発しているPyCallを使うと、RubyからPythonオブジェクトを操作できる。 Rubyから、Pythonの機械学習・統計分析のツールを利用することを目的としており、 ネット上にもnumpyやscikit-learnを実行する例があがっている。

Rubyist Magazine - PyCall があれば Ruby で機械学習ができる

このPyCallで、Kerasを叩くことができれば、RubyでもDeep Learningできると思い試してみた。 まずはインストールから。

$ gem install --pre pycall

ちなみに、実行環境をまとめると、Ruby 2.4.0、PyCall 0.1.0.alpha.20170317、Python 3.6.0、Theano 0.9.0、Keras 2.0.2である。 試したのは、MNISTの手書き数字画像を、畳み込みニューラルネットで認識するサンプルコードである。

keras/mnist_cnn.py at master · fchollet/keras · GitHub

これを、細かいところは適宜はぶきながら、Ruby+PyCallで移植すると次のようになる。

require 'pycall/import'
include PyCall::Import

# Kerasの必要なものをimportする．
pyimport 'keras'
pyfrom 'keras.datasets', import: 'mnist'
pyfrom 'keras.models', import: 'Sequential'
pyfrom 'keras.layers', import: ['Dense', 'Dropout', 'Flatten']
pyfrom 'keras.layers', import: ['Conv2D', 'MaxPooling2D']

# MNISTは28x28の大きさの手書き数字画像で、各画像が10個のクラスにわけられている．
nb_classes = 10
img_rows = 28
img_cols = 28

# MNISTデータセットを読み込む．初回実行時はダウンロードするところから始まる．
(x_train, y_train), (x_test, y_test) = mnist.load_data.()

# データのreshapeは、元のコードでは...
#   pyfrom 'keras', import: 'backend'
#   backend.image_data_format.()
# の結果で処理を分けている．
# 試してみたところ "channels_last" だったので、
# そちらの処理を移植した．
x_train = x_train.reshape.(x_train.shape[0], img_rows, img_cols, 1)
x_test = x_test.reshape.(x_test.shape[0], img_rows, img_cols, 1)
# 型をfloat32にして、要素を[0.0,1.0]にする．
x_train = x_train.astype.('float32')
x_test = x_test.astype.('float32')
x_train /= 255
x_test /= 255

# ラベル情報をクラスベクトル形式にする．
y_train = keras.utils.to_categorical.(y_train, nb_classes)
y_test = keras.utils.to_categorical.(y_test, nb_classes)

# ネットワークを定義する．
model = Sequential.()
model.add.(Conv2D.(32, kernel_size: [3, 3], activation: 'relu', 
                   input_shape: [img_rows, img_cols, 1]))
model.add.(Conv2D.(64, kernel_size: [3, 3], activation: 'relu'))
model.add.(MaxPooling2D.(pool_size: [2, 2]))
model.add.(Dropout.(0.25))
model.add.(Flatten.())
model.add.(Dense.(128, activation: 'relu'))
model.add.(Dropout.(0.5))
model.add.(Dense.(nb_classes, activation: 'softmax'))

# ネットワークをコンパイルする．初回実行時はそれなりに時間がかかる．
model.compile.(loss: keras.losses.categorical_crossentropy,
              optimizer: keras.optimizers.Adadelta.(),
              metrics: ['accuracy'])

# ネットワークを学習する．
model.fit.(x_train, y_train,
          batch_size: 128,
          epochs: 10,
          verbose: 1,
          validation_data: [x_test, y_test])

# 分類性能を評価する．
score = model.evaluate.(x_test, y_test, verbose: 0)
print(sprintf("Test loss: %.6f\n", score[0]))
print(sprintf("Test accuracy: %.6f\n", score[1]))

これを実行すると、問題なくネットワークの学習が動き始める！

$ ruby keras_test.rb
Using Theano backend.
Using cuDNN version 5005 on context None
Mapped name None to device cuda: GeForce GTX 1080 (0000:06:00.0)
Train on 60000 samples, validate on 10000 samples
Epoch 1/10
60000/60000 [==============================] - 27s - loss: 0.3227 - acc: 0.9030 - val_loss: 0.0730 - val_acc: 0.9760
...
Epoch 10/10
60000/60000 [==============================] - 26s - loss: 0.0410 - acc: 0.9881 - val_loss: 0.0298 - val_acc: 0.9890
Test loss: 0.029782
Test accuracy: 0.989000

Kerasは、ブロックをつなげる感じで、簡単にネットワーク構造を定義できる。 Pythonに詳しくないRubyエンジニアが、Deep Learningを試してみるには、PyCall+Kerasが最高な気がする。 また、Pythonで学習したネットワークを、Ruby+PyCallで読み込んで使うという形にすれば、 容易にDeep Learningな機能をRailsアプリに組み込める、という感じで夢が広がる。

scikit-learnで近似最近傍探索したいときはLSHForestがある

2017/8/13追記: LSHForestはパフォーマンスがよろしくないため、0.19からDEPRECATEDとなった。0.21から削除されるようなので、使用しないほうが良い。

scikit-learnでは、ver. 0.16から近似最近傍探索手法のLSHForestが実装されている。LSHForestは、ハッシングによる近似最近傍探索の代表的な手法であるLocality Sensitive Hashing（LSH）をベースにした、木構造の近似最近傍探索手法である。LSHは特徴ベクトルをRandom Projectionと閾値処理で0と1の短いバイナリベクトルに変換し、これをキーとしてハッシュテーブルを作る。LSHForestでは、LSHによって得られれたバイナリベクトルから木構造（あるビットが0か1で二分木が作れる、これをLSHTreeと呼ぶ）を複数個つくる。検索クエリが与えられると、それぞれのLSHTreeから候補を割り出して、それら候補を距離で並べることで検索結果を得る。ちなみに、コサイン距離（1-コサイン類似度）による検索しか行えない。この点では、FLANNの方が柔軟である。 それでは、検索インデックスを作成するスクリプトは、次のようになる。これをgen_idx.pyとする。

#!/usr/bin/env python
# -*- coding: utf-8 -*-

from sklearn.datasets import load_svmlight_file
from sklearn.neighbors import LSHForest
from sklearn.externals import joblib

def main():
  # MNIST（手書き数字画像）データセットを読み込む。
  # ※検索対象の例として使う。
  targets, _ = load_svmlight_file('mnist.scale')

  # LSHはアルゴリズム中でRandom Projectionを用いるので、
  # 再現性を確保したい場合は、random_stateに何か与えると良い。
  search_idx = LSHForest(random_state=1984)
  
  # 検索インデックスを作成する。
  search_idx.fit(targets)

  # 検索インデックスを保存する。
  joblib.dump(search_idx, 'lshtree.pkl.cmp', compress=True)

if __name__ == '__main__':
  main()

検索インデックスを読み込み、検索するスクリプトは次のようになる。これをsearch.pyとする。

#!/usr/bin/env python
# -*- coding: utf-8 -*-

from sklearn.datasets import load_svmlight_file
from sklearn.neighbors import LSHForest
from sklearn.externals import joblib

def main():
  # 検索クエリとして、MNISTのテストデータを読み込む。
  queries, q_labels = load_svmlight_file('mnist.scale.t')

  # 検索対象データのラベルを読み込む。
  # ※あとで検索結果を確認したいためで、検索には必要ない。
  _, t_labels = load_svmlight_file('mnist.scale')
  
  # 検索インデックスを読み込む。
  search_idx = joblib.load('lshtree.pkl.cmp')

  # 各クエリの5-近傍を検索する。
  dists, ids = search_idx.kneighbors(queries, n_neighbors=5)

  # 検索クエリの一番目のデータで、検索結果を確認する。
  print(q_labels[0])
  print(t_labels[ids[0,:]])
  print(ids[0,:])
  print(dists[0,:])

if __name__ == '__main__':
  main()

実行してみると、検索クエリと同じ「7」のMNISTデータが検索できているのがわかる。

$ ./gen_idx.py
$ ./search.py
7.0
[ 7.  7.  7.  7.  7.]
[15260 16186 14563  9724 31073]
[ 0.08516171  0.09843745  0.10159849  0.10445509  0.10882645]

LSHForestにも、LSHTreeの数などのパラメータがあるが、論文中で書かれている値であったりと、デフォルトで問題なさそう。 内部のLSHとしては、32ビットのバイナリベクトルに変換する様子。 FLANNと違って、検索時に、元の検索対象データを必要としないのが良い。 ちなみに、検索インデックスに新たにデータを加えるときは、partial_fitメソッドを使う。