酒と泪とRubyとRailsと

Ruby on Rails と Objective-C は酒の肴です!

【書評: データサイエンティスト養成読本】で学ぶデータサイエンティストの生態系

統計解析の勉強をしていく中で、どうやったら実際の業務の中に活かしていけるかを模索していて、この本にいきあたりました。人々が生み出すデータがますます増加していく中で、そのデータをビジネスにどう活用していくかをエンジニアの視点で考えられる素晴らしい本だと思います。

今回はこの本の中で特に参考になった点を中心にピックアップしていきます。


本書のサンプルデータ

サポートページ:データサイエンティスト養成読本 [ビッグデータ時代のビジネスを支えるデータ分析力が身につく!]
インターネットの力、マジですごい!サンプルデータはほんとうに有難い!

データサイエンティストの仕事術

Cross Industry Standard Process for Data Mining(CRISP-DM)

データマイニングプロジェクトにおける標準的な6つの手順。

1) Business Understanging
ビジネスの全体像の理解(ヒアリング)と課題設定

2) Data Understanding
データ収集と理解に努めて、仮説構築の下地を作る

3) Data Preparation
数理モデルの作成のために、データマート(データウェアハウス)を構築。
欠損値や外れ値を除去していく。

4) Modeling
仮説にもどついて、数理モデルを構築する

5) Evaluation
数理モデルの評価。ビジネスリスクの評価

6) Deployment
数理モデルに基づく、ビジネス施策の実施

必要なスキルセット

# ハードスキル
1) Data Understanding 〜 Data Preparation
RDBMS(SQL)、Hadoop(HDFS)、MapReduce、Hiveなど

2) Modeling 〜 Devaluation
統計解析、機械学習の知識、R, Pythonなど

# ソフトスキル
1) Business Understanding 〜 Data Understanding
ビジネスにおける業界、業務の知識。質問力、理解力が必要

2) Deployment
分析結果の伝達力、説得力、プロジェクトの推進力が必要

根本としては好奇心。「問題を核心まで掘り下げて、そこで見つけた疑問を明確で検証可能な一連の仮説に落としこみたいという欲求」が必要。

データサイエンスよりも大切なこと

1) 統計的な正しさよりもビジネスの成功
2) 組織の構造や実務の改善効果を重視する
3) 意味のある分析のためには泥臭いクリーニングが必須
4) データは適切な保持が大切

パラメトリックとノンパラメトリック

@TJO_datasciさんからご指摘頂きました。パラメトリックとノンパラメトリックの違いの端的なまとめは次の通りです。

* パラメトリック検定:正規分布に限らず厳密な分布の形の仮定を必要とする
* ノンパラメトリック検定:分布の形の仮定が不要である

引用元: 「パラメトリック検定」と「ノンパラメトリック検定」の違いについて出典を明示して書いておく

僕自身幾つか文献を見ましたが、まだまだ理解できていません。引き続き勉強して理解できたら、また自分の言葉で書き直したいと思います!

データ分析実践入門

事前準備: RStudio

この章では主にRを使います。RStudioについては次の章で説明します。

ロジスティック回帰

Rにデフォルトで入っている、タイタニック号の乗客の生存情報(Titanic)にロジスティック回帰モデルを使う。乗客の各属性から、生存の要因を探す。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
## データ整形
z <- data.frame(Titanic)
titanic.data <- data.frame(
    Class=rep(z$Class, z$Freq),
    Sex=rep(z$Sex, z$Freq),
    Age=rep(z$Age, z$Freq),
    Survived=rep(z$Survived, z$Freq)
  )
head(titanic.data)
# Class  Sex   Age Survived
# 1   3rd Male Child       No
# 2   3rd Male Child       No
# 3   3rd Male Child       No
# 4   3rd Male Child       No
# 5   3rd Male Child       No
# 6   3rd Male Child       No

## モデル構築
titanic.logit <- glm(Survived~., data=titanic.data, family=binomial)
summary(titanic.logit)
# Call:
#   glm(formula = Survived ~ ., family = binomial, data = titanic.data)
# 
# Deviance Residuals: 
#   Min       1Q   Median       3Q      Max  
# -2.0812  -0.7149  -0.6656   0.6858   2.1278  
# 
# Coefficients:
#   Estimate Std. Error z value Pr(>|z|)    
# (Intercept)   0.6853     0.2730   2.510   0.0121 *  
#   Class2nd     -1.0181     0.1960  -5.194 2.05e-07 ***
#   Class3rd     -1.7778     0.1716 -10.362  < 2e-16 ***
#   ClassCrew    -0.8577     0.1573  -5.451 5.00e-08 ***
#   SexFemale     2.4201     0.1404  17.236  < 2e-16 ***
#   AgeAdult     -1.0615     0.2440  -4.350 1.36e-05 ***
#   ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# (Dispersion parameter for binomial family taken to be 1)
# 
# Null deviance: 2769.5  on 2200  degrees of freedom
# Residual deviance: 2210.1  on 2195  degrees of freedom
# AIC: 2222.1
# 
# Number of Fisher Scoring iterations: 4

オッズ比の計算

1
2
3
4
5
6
7
8
9
10
install.packages('epicalc')
library('epicalc')
logistic.display(titanic.logit, simplified=T)
# OR lower95ci  upper95ci     Pr(>|Z|)
# Class2nd   0.3612825 0.2460447  0.5304933 2.053331e-07
# Class3rd   0.1690159 0.1207510  0.2365727 3.691891e-25
# ClassCrew  0.4241466 0.3115940  0.5773549 5.004592e-08
# SexFemale 11.2465380 8.5408719 14.8093331 1.431830e-66
# AgeAdult   0.3459219 0.2144193  0.5580746 1.360490e-05
#=> 例えば、女性のほうが11.2465380倍生存確率が高いということがわかる

決定木モデル

1
2
3
4
5
6
library(rpart)
install.packages('partykit')
library(partykit)

titanic.rp <- rpart(Survived~., data=titanic.data)
plot(as.party(titanic.rp), tp_args=T)

スクリーンショット 2014-03-25 22.45.03

主成分分析

1
2
state.pca <- prcomp(state.x77[,1:6], scale=T)
biplot(state.pca)

スクリーンショット 2014-03-25 22.48.23

非階層クラスタリング: k-means

k-meansの手順。

1) k個のクラスタの中心の初期値を決める
2) 各データをk個クラスタ中心との距離を求め、最も近いクラスタに分類
3) 形成されたクラスタの中心を求める
4) クラスタの中心が変化しない時点まで(2), (3)を繰り返す

サンプルソースはこちら。

1
2
3
4
5
6
7
8
9
10
# k-meansの実行(クラスタ数を3つに設定)
state.km <- kmeans(scale(state.x77[,1:6]),3)

# 主成分分析の結果にクラスタ情報を付与
state.pca.df <- data.frame(state.pca$x)
state.pca.df$name <- rownames(state.pca.df)
state.pca.df$cluster <- as.factor(state.km$cluster)

# 描画
ggplot(state.pca.df, aes(x=PC1, y=PC2, label=name, col=cluster)) + geom_text() + theme_bw(16)

スクリーンショット 2014-03-25 23.01.10

レーダーチャート

1
2
3
4
5
6
7
8
9
10
11
12
13
install.packages('fmsb')
library('fmsb')
df <- as.data.frame(scale(state.km$centers))
dfmax <- apply(df, 2, max) + 1
dfmin <- apply(df, 2, min) - 1
df <- rbind(dfmax, dfmin, df)

# 描画
radarchart(df, seg=5, plty=1, pcol=rainbow(3))
legend("topright", legend=1:3, col=rainbow(3), lty=1)
# クラスタ1 => 平均寿命、高卒率、収入が多い
# クラスタ2 => 殺人件数と非識字率が高い
# クラスタ3 => 人口が多い

スクリーンショット 2014-03-26 7.44.11

RStudioのすすめ

以下オススメリンク。

デスクトップアプリ R用のIDE: RStudio
RのIDEツール。主な機能、使い方は次の通り。

1) RStudioのコード補完や、関数へのジャンプ、組み込み関数のソース閲覧機能
2) Apacheとの連携によるレポート表示、Emailでのレポート配信
3) Gitによるソースコード管理のサポート機能

knitr
RStudio上でReproducible Researchを実現するためのパッケージ。Markdown, Rファイル => HTMLへの変換を行う事ができる。また、後述のRPubsへのアップロード機能も備える。

RPubs
RStudio.comが提供するサービス。再現可能な研究(Reproducible Research)を実現するために、Rの実行手順などをWebでシェアできるサービス。

RStudio Server
RStudioのWebアプリケーション版。

Pythonによる機械学習

Pythonで使えるデータ分析に有効なライブラリ群についての紹介。

データ分析と機械学習

データ分析において、機械学習は単なる手段の一つ。

* ある説明変数に対して目的変数がどのくらい影響を与えているかを知ることができる
* 複数の変数から重要な変数を絞り込む事ができる

Python環境の準備

もしMachのPython環境を作る場合は、拙著[Python開発環境構築 徹底解説pyenv, macがオススメです。

Python のライブラリ

利用するライブラリのインストール。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
# 効率的な数値計算のためのライブラリ
pip install numpy

# 最適化、補間、積分、統計、画像処理、クラスタリング解析、信号処理、空間的解析などができる
pip install scipy

# matplotlibのインストール用の設定
ln -s /usr/local/opt/freetype/include/freetype2 /usr/local/include/freetype

# グラフ描画ライブラリ
pip install matplotlib

# ファイルを作成してパラメータを記述
vim ~/.matplotlib/matplotlibrc
# ↓ 以下を追加
backend : TkAgg

# 機械学習ライブラリ: 単回帰分析〜SVMなどまで対応
# 教師あり学習、教師なし学習、モデル選択のためのクロスバリデーション
pip install scikit-learn

# データ構造・データ解析のライブラリ:Rのdata.frameなどがある
pip install pandas

# PythonからR言語を呼び出すためのライブラリ
pip install rpy2

# pythonを対話的に実行するためのライブラリ
pip install ipython

あやめのデータを知る

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
wget https://dataminingproject.googlecode.com/svn-history/r44/DataMiningApp/datasets/Iris/iris.csv

# IRIS.csvのあるフォルダでipythonを起動
ipython

# iris.csvを読み込み
import pandas as pd
iris = pd.read_csv("iris.csv")
iris.info()

# データ加工
setosa = iris[iris["Species"] == "setosa"]
versicolor = iris[iris["Species"] == "versicolor"]
virginica = iris[iris["Species"] == "virginica"]

# 基本的な統計量 ###########################
# 合計
virginica.sum()
# Sepal Length                                                329.4
# Sepal Width                                                 148.7
# Petal Length                                                277.6
# Petal Width                                                 101.3

# 平均
virginica.mean()
# Sepal Length    6.588
# Sepal Width     2.974
# Petal Length    5.552
# Petal Width     2.026

# 中央値
virginica.median()
# Sepal Length    6.50
# Sepal Width     3.00
# Petal Length    5.55
# Petal Width     2.00

# 最小値
verginica.min()
# Sepal Length          4.9
# Sepal Width           2.2
# Petal Length          4.5
# Petal Width           1.4

# 最大値
virginica.max()
# Sepal Length          7.9
# Sepal Width           3.8
# Petal Length          6.9
# Petal Width           2.5

# 相関係数
virginica.corr()
#               Sepal Length  Sepal Width  Petal Length  Petal Width
# Sepal Length      1.000000     0.457228      0.864225     0.281108
# Sepal Width       0.457228     1.000000      0.401045     0.537728
# Petal Length      0.864225     0.401045      1.000000     0.322108
# Petal Width       0.281108     0.537728      0.322108     1.000000

# 分散
virginica.var()
# Sepal Length    0.404343
# Sepal Width     0.104004
# Petal Length    0.304588
# Petal Width     0.075433

# 標準偏差
virginica.std()
# Sepal Length    0.635880
# Sepal Width     0.322497
# Petal Length    0.551895
# Petal Width     0.274650

# 共分散
virginica.cov()
#               Sepal Length  Sepal Width  Petal Length  Petal Width
# Sepal Length      0.404343     0.093763      0.303290     0.049094
# Sepal Width       0.093763     0.104004      0.071380     0.047629
# Petal Length      0.303290     0.071380      0.304588     0.048824
# Petal Width       0.049094     0.047629      0.048824     0.075433

# ピポッテーブルの作成
import numpy as np
pd.pivot_table(iris, rows=["Species"], aggfunc=np.mean)
#             Petal Length  Petal Width  Sepal Length  Sepal Width
# Species
# setosa             1.464        0.244         5.006        3.418
# versicolor         4.260        1.326         5.936        2.770
# virginica          5.552        2.026         6.588        2.974

irisでヒストグラムの実行

1
2
3
4
5
import matplotlib.pyplot as plt
plt.hist(iris["Sepal Length"])
plt.xlabel("Sepal Length")
plt.ylabl("Freq")
plt.show()

スクリーンショット 2014-03-29 9.37.35

箱ひげ図

1
2
3
4
5
6
7
data = [setosa["Petal Length"], versiclor["Petal Length"], virginica["Petal Length"]]
plt.boxplot(data)
plt.xlabel("Species")
plt.ylabel("Petal Length")
ax = plt.gca()
plt.setp(ax, xticklabels=["setosa", "versicolor", "virginica"])
plt.show()

スクリーンショット 2014-03-29 9.55.08

散布図

1
2
3
4
plt.scatter(versicolor["Sepal Width"], versicolor["Petal Width"])
plt.xlabel("Sepal Width")
plt.ylabel("Petal Width")
plt.show()

スクリーンショット 2014-03-29 10.01.24

相関係数

1
2
3
4
5
# versicolorのSepal WidthとPetal Widthの相関係数
import numpy as np
corr = np.corrcoef(versicolor["Sepal Width"], versicolor["Petal Width"])
print(corr[0,1])
#=> 0.663998720024

単回帰分析

1
2
3
4
5
6
7
8
9
10
11
12
from sklearn import linear_model
LinerRegr = linear_model.LinearRegression()
X = setosa[["Sepal Length"]]
Y = setosa[["Sepal Width"]]
LinerRegr.fit(X, Y)
plt.scatter(X, Y, color="black")
px = np.arange(X.min(), X.max(), .01)[:,np.newaxis]
py = LinerRegr.predict(px)
plt.plot(px, py, color="blue", linewidth=3)
plt.xlabel("Sepal Length")
plt.ylabel("Sepal Width")
plt.show()

スクリーンショット 2014-03-29 11.43.26

1
2
3
4
5
6
7
# 回帰変数
print(LinerRegr.coef_)
#=> [[ 0.80723367]]

# 切片
print(LinerRegr.inercept_)
#=> [-0.62301173]

Rの機能を使って重回帰分析 => rpy2

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
import rpy2.robjects as robjects
robjects.globalenv["SepalLength"] = robjects.FloatVector(setosa["Sepal Length"])
robjects.globalenv["SepalWidth"] = robjects.FloatVector(setosa["Sepal Width"])
robjects.globalenv["PetalLength"] = robjects.FloatVector(setosa["Petal Length"])
robjects.globalenv["PetalWidth"] = robjects.FloatVector(setosa["Petal Width"])
r = robjects.r

# Rの重回帰分析
rlm = r.lm("SepalWidth~SepalLength+PetalLength+PetalWidth")
print((r.summary(rlm)))
# Residuals:
#      Min       1Q   Median       3Q      Max
# -0.75025 -0.17983 -0.00411  0.15933  0.57704
#
# Coefficients:
#             Estimate Std. Error t value Pr(>|t|)
# (Intercept) -0.48721    0.56900  -0.856    0.396
# SepalLength  0.79304    0.11162   7.105 6.35e-09 ***
# PetalLength -0.09678    0.22875  -0.423    0.674
# PetalWidth   0.31530    0.37186   0.848    0.401
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#
# Residual standard error: 0.2594 on 46 degrees of freedom
# Multiple R-squared:  0.5649,  Adjusted R-squared:  0.5366
# F-statistic: 19.91 on 3 and 46 DF,  p-value: 2.042e-08

ダミー変数の利用

1
2
3
4
5
# 数値データ以外も分析に活用する
# Spiciesをダミー変数(品種の0/1のカラム)に変換
import pandas as pd
dummies = pd.get_dummies(iris["Species"])
iris = pd.concat([iris, dummies], axis=1)

回帰係数と切片を求める

1
2
3
4
5
6
7
8
9
10
LinearRegr = linear_model.LinearRegression()
X = iris[["setosa", "versicolor"]]
Y = iris[["Sepal Length"]]
LinearRegr.fit(X, Y)

# 回帰係数
print(LinearRegr.coef_) #=> [[-1.582 -0.652]]

# 切片
print(LinearRegr.intercept_) #=> [ 6.588]

 データマイニングの流れ

1) ゴールを決める
2) 変数の種類・表している事象を知る
3) 1変量の解析を行う。平均値を算出し、ヒストグラム・箱ひげ図を描写
4) 2変量の解析を行う。散布図などのグラフから理解を深める
5) 説明変数を減らして、解釈しやすいモデルに当てはめる
6) 推定された結果を(2)-(4)に振り返ってモデル検証を行う

ロジスティック回帰モデル

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
import sklearn
use_data = np.logical_or(iris["Species"] == "setosa", iris["Species"] == "virginica")
setosa_virginica = iris[use_data]
X = setosa_virginica[["Sepal Length", "Sepal Width"]]
Y = setosa_virginica[["setosa"]]
LogRegr = sklearn.linear_model.LogisticRegression(C=1.0)
LogRegr.fit(X,Y)

# 偏回帰係数
print(LogRegr.coef_)
#=> [[-2.15056433  3.53948104]]

# 切片
print(LogRegr.intercept_)
#=> [ 0.94024115]

なんか間違っている気がするので、後で直す予定です。まずはちょっと触れるところから。

モデル結果をグラフにして確認

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
xMin = X["Sepal Length"].min()
xMax = X["Sepal Length"].max()
yMin = X["Sepal Width"].min()
yMax = X["Sepal Width"].max()
xx,yy = np.meshgrid(np.arange(xMin, xMax, 0.01), np.arange(yMin, yMax, 0.01))
Z = LogRegr.predict(np.c_[xx.ravel(), yy.ravel()])
Z = Z.reshape(xx.shape)
plt.figure()
plt.xlim(xx.min(), xx.max())
plt.ylim(yy.min(), yy.max())
plt.xlabel("Sepal Length")
plt.ylabel("Sepal Width")
plt.pcolormesh(xx, yy, Z, cmap=plt.cm.Paired)
plt.scatter(X["Sepal Length"], X["Sepal Width"], c=np.array(Y), cmap=plt.cm.Paired)
plt.show()

スクリーンショット 2014-03-29 17.36.16

決定木

まずはhomebrewでGraphvizとpythonで使うためのライブラリpydotをインストール。

1
2
3
4
5
6
7
brew update
brew install gts
brew install graphviz

pip uninstall pyparsing
pip install -Iv https://pypi.python.org/packages/source/p/pyparsing/pyparsing-1.5.7.tar.gz#md5=9be0fcdcc595199c646ab317c1d9a709
pip install pydot

続いてIPythonのコンソール。

1
2
3
4
5
from sklearn import tree
X = iris[["Sepal Length", "Sepal Width", "Petal Length", "Petal Width"]]
Y = iris[["Species"]]
treeClf = tree.DecisionTreeClassifier(max_depth=2)
treeClf.fit(X, Y)

決定木

同じくグラフがなんか違うので、後で見直す予定。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
## 決定木の構築
from sklearn import tree
X = iris[["Sepal Length", "Sepal Width", "Petal Length", "Petal Width"]]
Y = iris[["Species"]]
treeClf = tree.DecisionTreeClassifier(max_depth=2)
treeClf.fit(X, Y)

## 決定木の可視化
import StringIO
import pydot
dot_data = StringIO.StringIO()
tree.export_graphviz(treeClf, out_file=dot_data, feature_names=["Sepal Length", "Sepal Width", "Petal Length", "Petal Width"])
graph = pydot.graph_from_dot_data(dot_data.getvalue())
graph.write_pdf("iris_decision_tree.pdf")

スクリーンショット 2014-03-29 18.43.37

k-means法

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
# k-meansの実行
from sklearn import cluster
X = iris[["Sepal Length", "Sepal Width"]]
kmeansCls = cluster.KMeans(n_clusters=3)
kmeansCls.fit(X)
print(kmeansCls.predict(X))
# [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
#  0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
#  1 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 2 2 2 2 1 2 2 2 2
#  2 2 1 1 2 2 2 2 1 2 1 2 1 2 2 1 1 2 2 2 2 2 1 2 2 2 2 1 2 2 2 1 2 2 2 1 2
#  2 1]

# k-means法の可視化
import numpy as np
import matplotlib.pyplot as plt
def category2int(x):
    category = {"setosa": 0, "versicolor": 1, "virginica": 2}
    return category[x]

f = lambda x: category2int(x)
Y = iris["Species"].map(f)
xMin = X["Sepal Length"].min()
xMax = X["Sepal Length"].max()
yMin = X["Sepal Width"].min()
yMax = X["Sepal Width"].max()
xx,yy = np.meshgrid(np.arange(xMin, xMax, 0.01), np.arange(yMin, yMax, 0.01))
Z = kmeansCls.predict(np.c_[xx.ravel(), yy.ravel()])
Z = Z.reshape(xx.shape)
plt.figure()
plt.xlim(xx.min(), xx.max())
plt.ylim(yy.min(), yy.max())
plt.xlabel("Sepal Length")
plt.ylabel("Sepal Width")
plt.pcolormesh(xx, yy, Z, cmap=plt.cm.Reds)
plt.scatter(X["Sepal Length"], X["Sepal Width"], c=np.array(Y), cmap=plt.cm.Blues)
plt.show()

スクリーンショット 2014-03-29 18.40.37

データマイニング 10のアルゴリズム

C4.5

* 決定木を構築する計算アルゴリズム
* 結果のイメージが直感的て理解しやすい
* 分類の精度は高くない

k-meansアルゴリズム

クラスタリングでは次のステップを反復して計算。

1) データにクラスタを割り振る(中心と一番近いクラスタを割り振る)
2) クラスタ事の平均値の計算

サポートベクターマシン

* 少なめのサンプルで説明変数が多い場合に適した手法
* SVMは特徴量が膨大でも高速に計算できる
* サンプル数が増えると計算時間が急に増加する

アプリオリアルゴリズム

* 大量のトランザクションの中から価値あるつながりを見つけるために使われる手法
* 強力な結果が得られる上、実装が比較的容易

EMアルゴリズム

Expectation-maximization algorithm(期待値最大化法)。数式展開では計算出来ないような複雑なモデルのパラメータ推定に用いる。以下はその手順。

1) パラメータの初期値を決める
2) Exceptationステップ: 未観測データの期待値を求める
3) Maximizationステップ: 未観測データの期待値を求める

アダブースト

* いくつかの学習木を組み合わせることで、強力な予測性能を得る
* 実装が容易で解釈がしやすい
* ディープラーニングが出るまでは最も有力な予測モデルと言われてきた

ナイーブベイズ

* クラスを予測するための手法で構築が容易
* テキスト分類やパターン認識で広く利用されている

Rによるマーケティング分析

R言語で学ぶマーケティング分析 競争ポジショニング戦略
この記事の作者さんの発表。エッセンスがうまくまとまっています。

主成分分析の実行

1
2
3
4
5
6
7
8
9
setwd("~/Dropbox/study-R/")
# データの読み込み
sp.user.data <- read.csv("data/sp_user_research_data.csv")

# 主成分分析 => ユーザー間の類似関係の把握
sp.user.pca <- prcomp(sp.user.data[,-1], scale =T)

# バイプロットの表示
biplot(sp.user.pca)

スクリーンショット 2014-03-30 13.33.16

k-meansによるクラスタリングの実行

1
2
3
4
5
6
7
8
9
10
11
sp.user.km <- kmeans(sp.user.data[,-1], 4)

# 主成分分析の結果にクラスタ情報を付与
sp.user.pca.df <- data.frame(sp.user.pca$x)
sp.user.pca.df$id <- sp.user.data$id
sp.user.pca.df$cluster <- as.factor(sp.user.km$cluster)

# 描画
install.packages("ggplot2")
library("ggplot2")
ggplot(sp.user.pca.df, aes(x = PC1, y = PC2, label = id, col = cluster)) + geom_text() + theme_bw(16)

スクリーンショット 2014-03-30 13.40.50

レーダーチャートの作成

1
2
3
4
5
6
7
8
9
10
11
12
install.packages("fmsb")
library("fmsb")

# レーダーチャート用にデータを整形
df <- data.frame(scale(sp.user.km$centers))
dfmax <- apply(df, 2, max) + 1
dfmin <- apply(df, 2, min) - 1
df <- rbind(dfmax, dfmin, df)

# レーダーチャートの描画
radarchart(df, seg = 5, plty = 1, pcol = rainbow(4))
legend("topright", legend = 1:4, col = rainbow(4), lty = 1)

スクリーンショット 2014-03-30 13.54.09

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
## MDS(Mult-dimentional scaling: 多次元尺度構成法) - 知覚マップ
install.packages("MASS")
library("MASS")

# 先行データの読み込み
setwd("~/Dropbox/study-R/data")
target.data <- read.csv("target_preference_data.csv", header = T)

# 非計算MDSの実行
service.dist <- dist(t(target.data[, -1]))
service.map <- isoMDS(service.dist)

# 描画用データの整形
service.map.df <- data.frame(scale(service.map$points))
service.map.df$service_name <- names(target.data[, -1])

# 描画
ggplot(service.map.df, aes(x = X1, y = X2, label = service_name)) + geom_text() + theme_bw(16)

スクリーンショット 2014-03-30 14.07.50

選好ベクトルの推定

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
# 選好ベクトルの推定
user.preference.data <-
  do.call(rbind,
          lapply(1:nrow(target.data),
                 function(i){
                   preference.data <- data.frame(
                     p=as.numeric(target.data[i, -1]),
                     X1=service.map.df$X1,
                     X2=service.map.df$X2)
                   fit <- lm(p~., data=preference.data)
                   b <- 2/sqrt(fit$coef["X1"]^2+fit$coef["X2"]^2)
                   data.frame(X1=b*fit$coef["X1"],
                              X2=b*fit$coef["X2"],
                              service_name=i)
                 }))

# 選好ベクトルの描画
ggplot(service.map.df, aes(x=X1, y=X2, label=service_name)) +
  geom_text() +
  theme_bw(16) +
  xlim(-2,2) +
  ylim(-2,2) +
  geom_point(data=user.preference.data, aes(x=X1, y=X2))

スクリーンショット 2014-03-30 14.27.52

仮想データの読み込みと散布図の描画

1
2
3
4
5
6
7
8
9
10
11
12
13
14
# 利用する仮想データの読み込み
library(ggplot2)
library(scales)

# GRPと売上データの読み込み
grp.data <- read.csv("http://image.gihyo.co.jp/assets/files/book/2013/978-4-7741-5896-9/download/grp.csv", header = T)

# 散布図の描画
ggplot(grp.data, aes(x=grp, y=amount)) +
  geom_point() +
  scale_y_continuous(label=comma, limits=c(0,360000)) +
  ylab("Sales") +
  xlab("GRP") +
  theme_bw(16)

スクリーンショット 2014-03-30 20.50.20

線形型モデルの構築

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
# 可視化はgeom_smooth関数ですぐにできる
ggplot(grp.data, aes(x=grp, y=amount)) +
  geom_point() +
  scale_y_continuous(label=comma, limits=c(0, 360000)) +
  ylab("Sales") +
  xlab("GRP") +
  geom_smooth(method = "lm") +
  theme_bw(16)

# モデル構築
fit <- lm(amount~grp, data=grp.data)

# モデル概要の表示
summary(fit)
# Call:
#   lm(formula = amount ~ grp, data = grp.data)
# 
# Residuals:
#   Min     1Q Median     3Q    Max 
# -31521  -8275    321   8701  36962 
# 
# Coefficients:
#   Estimate Std. Error t value Pr(>|t|)    
# (Intercept) 229812.53    2561.73   89.71   <2e-16 ***
#   grp            455.54      15.93   28.60   <2e-16 ***
#   ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# Residual standard error: 13070 on 198 degrees of freedom
# Multiple R-squared:  0.8051,  Adjusted R-squared:  0.8042 
# F-statistic: 818.1 on 1 and 198 DF,  p-value: < 2.2e-16

スクリーンショット 2014-03-30 20.58.40

xx減型モデルの構築

読めなかった。日本語苦手っすw

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
fit <- lm(log(amount)~log(grp), data=grp.data)

# モデル概要
summary(fit)
# Call:
#   lm(formula = log(amount) ~ log(grp), data = grp.data)
# 
# Residuals:
#   Min        1Q    Median        3Q       Max 
# -0.121782 -0.023589 -0.003298  0.025672  0.116504 
# 
# Coefficients:
#   Estimate Std. Error t value Pr(>|t|)    
# (Intercept) 11.546963   0.032379  356.62   <2e-16 ***
#   log(grp)     0.213933   0.006551   32.66   <2e-16 ***
#   ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# Residual standard error: 0.04064 on 198 degrees of freedom
# Multiple R-squared:  0.8434,  Adjusted R-squared:  0.8426 
# F-statistic:  1067 on 1 and 198 DF,  p-value: < 2.2e-16

# 予測結果を描画
fit.data <- data.frame(grp=grp.data$grp, amount=exp(fit$fitted.values))
ggplot(grp.data, aes(x=grp, y=amount)) +
  geom_point() +
  geom_line(data=fit.data, aes(x=grp, y=amount)) +
  scale_y_continuous(label=comma, limits=c(0,360000)) +
  ylab("Sales") +
  xlab("GRP") +
  theme_bw(16)

スクリーンショット 2014-03-30 21.19.14

AとBが同じ場合のA/Bテストのシミュレーション

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
library("plyr")

# A/Bの心のコンバージョン率を設定
A1.CVR <- 0.09
B1.CVR <- 0.09

# サンプル数
n <- 10000
set.seed(2)

# シミュレーション
AB1 <- data.frame(Pattern=c(rep("A", n), rep("B", n)), CV=c(rbinom(n, 1, A1.CVR), rbinom(n, 1, B1.CVR)))

# コンバージョン率の算出
ddply(AB1, .(Pattern), summarize, CVR=mean(CV))
# Pattern    CVR
# 1       A 0.0942
# 2       B 0.0912

# カイ二乗検定
chisq.test(table(AB1))
# Pearson's Chi-squared test with Yates' continuity correction
# 
# data:  table(AB1)
# X-squared = 0.5, df = 1, p-value = 0.4795

Bのほうが高い場合のA/Bテストのシミュレーション

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
A2.CVR <- 0.095
B2.CVR <- 0.097
n <- 10000
set.seed(2)

# シミュレーション
AB2 <- data.frame(Pattern=c(rep("A", n), rep("B", n)), CV=c(rbinom(n, 1, A2.CVR), rbinom(n, 1, B2.CVR)))

# コンバージョン率の算出
ddply(AB2, .(Pattern), summarize, CVR = mean(CV))
# Pattern   CVR
# 1       A 0.100
# 2       B 0.097

# カイ二乗検定
chisq.test(table(AB2))
# Pearson's Chi-squared test with Yates' continuity correction
# 
# data:  table(AB2)
# X-squared = 0.4735, df = 1, p-value = 0.4914

50万件でのシミュレーション

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
A3.CVR <- 0.095
B3.CVR <- 0.097
n <- 500000
set.seed(2)

# シミュレーション
AB3 <- data.frame(Pattern=c(rep("A", n), rep("B", n)),
                  CV=c(rbinom(n, 1, A3.CVR), rbinom(n, 1, B3.CVR)))

# コンバージョン率の算出
ddply(AB3, .(Pattern), summarize, CVR=mean(CV))
# Pattern      CVR
# 1       A 0.095192
# 2       B 0.097040

# カイ二乗検定の実行
chisq.test(table(AB3))
# Pearson's Chi-squared test with Yates' continuity correction
# 
# data:  table(AB3)
# X-squared = 9.8061, df = 1, p-value = 0.001739

直交表の作成

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
install.packages("conjoint")
library(conjoint)

# 構成要素
experiment <- expand.grid(
  imgA = c("ImageA1", "ImageA2"),
  imgB = c("ImageB1", "ImageB2"),
  txtA = c("TextA1", "TextA2"),
  txtB = c("TextB1", "TextB2"))

# 交流表の作成
design.ort <- caFactorialDesign(data= experiment, type="orthogonal")
design.ort
# imgA    imgB   txtA   txtB
# 2  ImageA2 ImageB1 TextA1 TextB1
# 3  ImageA1 ImageB2 TextA1 TextB1
# 5  ImageA1 ImageB1 TextA2 TextB1
# 8  ImageA2 ImageB2 TextA2 TextB1
# 9  ImageA1 ImageB1 TextA1 TextB2
# 12 ImageA2 ImageB2 TextA1 TextB2
# 14 ImageA2 ImageB1 TextA2 TextB2
# 15 ImageA1 ImageB2 TextA2 TextB2

ロジスティック回帰モデルの構築

1
2
3
4
5
6
7
8
9
10
11
# 仮想のテスト結果データの読み込み(csがないので以下実行していません)
web.test.data <- read.csv("web_test_sample.csv")

# 上位6件だけ表示
head(web.test.data)

# ロジスティック回帰モデルの構築
fit <- step(glm(cv~., data = web.test.data[,-1], family=binomial))

# モデル概要
summary(fit)

Fluentd 入門

ログ回収をストリーミングで行ってくれるツール。ツールについての詳しい説明は、創業者の古橋さんのインタビュー記事『グリー技術者が聞いた、fluentdの新機能とTreasure Data古橋氏の野心』がわかりやすくてオススメです。

ログ解析の従来の課題

1) ログデータの転送に長い時間が必要、帯域の局所的な消費
2) データが解析可能になるまでのタイムラグ
3) ネットワーク障害などによる不安定さ

Fluentdの利点

1) JSON形式で構造化されている
2) 豊富なプラグイン(プラガブル)
3) 導入が容易
4) 耐障害性、バッファ、レジューム機能、可用性

アーキテクチャ

1) input: さまざまなソースからデータを読み取りFluentd用の形式に変換
2) Engine: 入ってきたデータをどのBuffe/Outpuプラグインに渡すかをルーティング
3) Bufferプラグイン: Outputプラグインの前にデータを一時保持する仕組み
4) Outputプラグイン: Bufferプラグインから渡されたデータを書き出す

ログ構造

tag: ログの種類。ルーティングで利用
time: ログが記録されたUNIX時刻
record: JSON形式のログ内容

活用方法

1) ログ収集: ログを集めてS3やHadoopに渡す
2) 統計 + 可視性: ログの統計情報を取得し、グラフ表示系のツールと連携する
3) 検知/通知: しきい値ベースや、データマイニングの異常値検出の通知など

システム構成

forwarder: ログ発生場所に立ち上げてログ回収 => aggregatorに送信
aggregator: 複数workerのためのロードバランシング、複数台構成による冗長化
worker: ログのパースや書き換えなど負荷の高い処理を行う
serializer: workerの処理結果のログを集約して別システムに書き出す
watcher: ログからメトリックスを収集、検知や通知を行う

事例紹介

1) ログの回収: DSPの広告効果のレポーティングのための収集の仕組み
2) nginx ログの可視化: Real Time Bidding(RTB)のツール

あとがき

中身濃すぎる!超面白かったです^^

Speical Thanks

python - Sudo pip install matplotlib fails to find freetype headers. [OS X Mavericks / 10.9] - Stack Overflow

python - pydot and graphviz error: Couldn’t import dot_parser, loading of dot files will not be possible - Stack Overflow

おすすめの書籍