R – ページ 5 – かものはしの分析ブログ

RstanでCVRの前後比較をするためのコード

目的

データサイエンス界隈の方がP値での意思決定に警鐘を鳴らしている昨今、施策実施に関するCVRの前後比較をχ2乗検定のP値を用いるのではなく、ベイズ統計学によるアプローチにチャレンジしてみたいと思いました。『基礎からのベイズ統計学』の8章で取り上げられていた比率データに対してのベイズ統計学的アプローチをもとに、stanを用いて事後分布から意思決定をするための進め方を紹介します。

進め方

・データの整形
・stanコード作成
・rstanでの引数の指定
・rでの可視化

データの特徴

Webマーケティング界隈では大変に多用するデータだと思いますが、実験を行ったユーザーに対しての開封・非開封、これまで通りのユーザーの開封・非開封の自然数からなるデータです。

stanコード

stanコードは
・dataブロック
・parametersブロック
・transformed parametersブロック（今回は不使用）
・modelブロック
・generated quantitiesブロック
からなります。
今回は自然数のデータであることから、ディリクレ分布を事前分布に設定するために、parametersブロックにおいてsimplexを指定しています。（教科書の比率データのものをそのまま使っています。）
modelは二項データしか出てこないので、二項分布を用いています。generated quantitiesブロックでは各々の比率、比率の差、比率の差が0を超える確率・0.01を超える確率、リスク比、リスク比が1を超える確率、オッズ比などを出力するようにしています。

data{
	int<lower=0> N[2]; #自然数からなるベクトルNの指定
	int n[2,2]; #整数からなる行列nの指定
}

parameters{
	simplex[2] p[2]; #ディリクレ分布を事前分布に設定したもとでの確率pを指定
}

model{
	for(i in 1:2){
		for(j in 1:2){
			n[i,j] ~ binomial(N[j], p[j][i]); #二項分布
		}
	}
}
generated quantities{
	real p11;
	real p10;
	real p01;
	real p00;
	real d;
	real delta_over;
	real delta_over_onep;
	real RR;
	real RRover;
	real OR;
	p11 <- p[1][1];
	p10 <- p[1][2];
	p01 <- p[2][1];
	p00 <- p[2][2];
	d <- p11 - p01; #比率の差
	delta_over <- if_else(d > 0,1,0);
	delta_over_onep <- if_else(d > 0.01,1,0);
	RR <- p11/p01; #リスク比
	RRover <- if_else(RR > 2,1,0);
	OR <- (p11/p10) / (p01/p00); #オッズ比
}

data{

int<lower=0> N[2]; #自然数からなるベクトルNの指定

int n[2,2]; #整数からなる行列nの指定

}

parameters{

simplex[2] p[2]; #ディリクレ分布を事前分布に設定したもとでの確率pを指定

}

model{

for(i in 1:2){

for(j in 1:2){

n[i,j] ~ binomial(N[j], p[j][i]); #二項分布

}

generated quantities{

real p11;

real p10;

real p01;

real p00;

real d;

real delta_over;

real delta_over_onep;

real RR;

real RRover;

real OR;

p11 <- p[1][1];

p10 <- p[1][2];

p01 <- p[2][1];

p00 <- p[2][2];

d <- p11 - p01; #比率の差

delta_over <- if_else(d > 0,1,0);

delta_over_onep <- if_else(d > 0.01,1,0);

RR <- p11/p01; #リスク比

RRover <- if_else(RR > 2,1,0);

OR <- (p11/p10) / (p01/p00); #オッズ比

}

rコード

以下は、stanをrで実行し、ggplot2などで可視化するためのコードが記されています。

library(rstan)
library (reshape)
library (ggplot2)

scr <- "model871.stan" #stanコード名

#カウントデータを用意します。
N <- c(1123, 1200)
n <- structure(.Data = c(106, 1017, 46, 1154), .Dim = c(2, 2))

#stanで用いるデータの型に変換
data <-list(N=N, n=n)

#パラメータの設定
par<-c("p","d","delta_over","delta_over_onep","RR","RRover","OR")
#pは各々の比率
#dは比率の差
#delta_overは比率の差が0を超える確率
#delta_over_onepは比率の差が0.1を超える確率
#RRはリスク比
#RRoverはRRが1を超える確率
#ORはオッズ比

war<-1000               #バーンインの期間を指定しています。
ite<-11000              #試行回数をしていしています。
see<-12345              #乱数の種
dig<-3                  #有効数字
cha<-1                  #連鎖構成数

#stanの実行
fit <- stan(file = scr, data = data, warm=war, iter=ite, seed=see,
	      pars=par,chains=cha)

#結果の出力
print(fit,pars=par,digits_summary=dig)

#事後分布の可視化
post   <- extract (fit, permuted = F)
m.post <- melt (post)
graph  <- ggplot (m.post, aes(x = value))
graph  <- graph + geom_density () + facet_grid(. ~ parameters, scales = "free") + theme_bw()
plot (graph)

library(rstan)

library (reshape)

library (ggplot2)

scr <- "model871.stan" #stanコード名

#カウントデータを用意します。

N <- c(1123, 1200)

n <- structure(.Data = c(106, 1017, 46, 1154), .Dim = c(2, 2))

#stanで用いるデータの型に変換

data <-list(N=N, n=n)

#パラメータの設定

par<-c("p","d","delta_over","delta_over_onep","RR","RRover","OR")

#pは各々の比率

#dは比率の差

#delta_overは比率の差が0を超える確率

#delta_over_onepは比率の差が0.1を超える確率

#RRはリスク比

#RRoverはRRが1を超える確率

#ORはオッズ比

war<-1000 #バーンインの期間を指定しています。

ite<-11000 #試行回数をしていしています。

see<-12345 #乱数の種

dig<-3 #有効数字

cha<-1 #連鎖構成数

#stanの実行

fit <- stan(file = scr, data = data, warm=war, iter=ite, seed=see,

pars=par,chains=cha)

#結果の出力

print(fit,pars=par,digits_summary=dig)

#事後分布の可視化

post <- extract (fit, permuted = F)

m.post <- melt (post)

graph <- ggplot (m.post, aes(x = value))

graph <- graph + geom_density () + facet_grid(. ~ parameters, scales = "free") + theme_bw()

plot (graph)

推定結果&可視化

今回の例では、実験を行ったユーザーのCVRの差が0以上の確率（delta_over）が1.0なので、ほぼ確実に差があると言えそうです。0.01以上差がある確率も1.0なので1%以上は差があると言えそうです。リスク比（RR）に関しては2.47と実験しない場合と比べて2.47倍程度CVを高めています。オッズ比（OR）は2.63とあるので、実験によるCV増大効果が2.63倍あると考えることができます。χ2乗検定では、二つの集団が独立かどうかを検定していますが、ベイズ統計学に従えば、「1%を超える確率」を算出することが容易なので、ディレクターなどに説明する際は圧倒的に理解を得られそうな気がします。

> print(fit,pars=par,digits_summary=dig)
Inference for Stan model: model871.
1 chains, each with iter=11000; warmup=1000; thin=1; 
post-warmup draws per chain=10000, total post-warmup draws=10000.

                 mean se_mean    sd  2.5%   25%   50%   75% 97.5% n_eff Rhat
p[1,1]          0.095   0.000 0.006 0.083 0.091 0.095 0.099 0.107  5645    1
p[1,2]          0.905   0.000 0.006 0.893 0.901 0.905 0.909 0.917  5645    1
p[2,1]          0.039   0.000 0.004 0.031 0.036 0.039 0.041 0.047  5488    1
p[2,2]          0.961   0.000 0.004 0.953 0.959 0.961 0.964 0.969  5488    1
d               0.056   0.000 0.007 0.042 0.051 0.056 0.061 0.071  5662    1
delta_over      1.000   0.000 0.000 1.000 1.000 1.000 1.000 1.000 10000  NaN
delta_over_onep 1.000   0.000 0.000 1.000 1.000 1.000 1.000 1.000 10000  NaN
RR              2.477   0.004 0.303 1.933 2.267 2.453 2.662 3.136  5454    1
RRover          0.955   0.003 0.207 0.000 1.000 1.000 1.000 1.000  5749    1
OR              2.633   0.005 0.341 2.026 2.396 2.606 2.842 3.377  5463    1

Samples were drawn using NUTS(diag_e) at Sun Mar 13 23:19:36 2016.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).

> print(fit,pars=par,digits_summary=dig)

Inference for Stan model: model871.

1 chains, each with iter=11000; warmup=1000; thin=1;

post-warmup draws per chain=10000, total post-warmup draws=10000.

mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat

p[1,1] 0.095 0.000 0.006 0.083 0.091 0.095 0.099 0.107 5645 1

p[1,2] 0.905 0.000 0.006 0.893 0.901 0.905 0.909 0.917 5645 1

p[2,1] 0.039 0.000 0.004 0.031 0.036 0.039 0.041 0.047 5488 1

p[2,2] 0.961 0.000 0.004 0.953 0.959 0.961 0.964 0.969 5488 1

d 0.056 0.000 0.007 0.042 0.051 0.056 0.061 0.071 5662 1

delta_over 1.000 0.000 0.000 1.000 1.000 1.000 1.000 1.000 10000 NaN

delta_over_onep 1.000 0.000 0.000 1.000 1.000 1.000 1.000 1.000 10000 NaN

RR 2.477 0.004 0.303 1.933 2.267 2.453 2.662 3.136 5454 1

RRover 0.955 0.003 0.207 0.000 1.000 1.000 1.000 1.000 5749 1

OR 2.633 0.005 0.341 2.026 2.396 2.606 2.842 3.377 5463 1

Samples were drawn using NUTS(diag_e) at Sun Mar 13 23:19:36 2016.

For each parameter, n_eff is a crude measure of effective sample size,

and Rhat is the potential scale reduction factor on split chains (at

convergence, Rhat=1).

参考文献

基礎からのベイズ統計学: ハミルトニアンモンテカルロ法による実践的入門
 rstanでちょこちょこ

RのChannelAttributionパッケージでマルコフ連鎖モデルのアトリビューション分析にチャレンジ

【目次】
・アトリビューションモデルとは
・アトリビューションモデルにおける評価の種類
・均等配分モデルと数理モデル
・ChannelAttributionパッケージについて
・Rでの実行
・実行結果
・参考文献

アトリビューションモデルとは

コンバージョンに至った広告の貢献度を、コンバージョンした経路のタッチポイントにどのように割り振るかを決めるルール。直接コンバージョンに至らなかったが、コンバージョンに繋がる働きをした広告の評価ができるモデルです。

アトリビューションモデルにおける評価の種類

Googleアナリティクスに載っていたアトリビューションモデルの種類です。一例なので、これが全てだとは言えないです。

終点アトリビューションモデル

　コンバージョンに繋がった最後の経路を評価

最後の間接クリックアトリビューションモデル

　コンバージョンに至る経路の直前の経路を評価

起点アトリビューションモデル

　最初に関わった経路を評価

線形アトリビューションモデル

　コンバージョンに至るまでに関わった経路を均等に評価

減衰アトリビューションモデル

　コンバージョンした際に、時間的に最も近いものを評価

接点ベースアトリビューションモデル

　起点と終点それぞれに40%ずつ評価、残りの20%は中間点に対して均等に評価

均等配分モデルと数理モデル

均等配分モデル

　コンバージョンに至るまでに関わった経路を均等に評価。数学に明るくなくても実行でき、集計が比較的簡単かつ解釈から施策までが容易とされる。一方で、恣意的な判断が入り込みやすかったり、各広告が独立であるという過程があるため、任意の広告のIMPが変化したときに他の広告への影響が加味できないなど欠点も多い。他の評価方法との比較のためのベンチマークとして用いられる。

数理モデル

マルコフ連鎖モデル
　ユーザーの2点間の行動のみに基づいて分析を行う。
　例えば、広告の経路が10通りある場合、広告1を踏んだユーザーは確率p11で広告1を次に触れる、確率p12で広告2を次に触れるというような確率で広告の接触確率が決まるモデルとなる。確率から接触したページの価値を計算することで、評価が可能となる。複雑怪奇なユーザー行動を2点間のデータに落とし込むことで、分析する際の情報量を制限することができる。

ChannelAttributionパッケージについて

Package ‘ChannelAttribution’
※Rのバージョンが古いとインストールできません。3.2以上にしましょう。

今回はRのChannelAttributionパッケージを用いて、マルコフ連鎖モデルを実践してみます。
このChannelAttributionパッケージにはデータセットがついているため、その用意されたデータを用います。

実際、このようなデータを元に分析を進めます。

内訳としては以下の通りです。
・コンバージョンに至ったパス
・コンバージョン数
・コンバージョン価値(価値に関してはサイトによって異なってくると思います。)
・コンバージョンに至らなかったパスの数(計算上使っていないようです。除外しても結果が変わらなかったので。)

Rでの実行

#ライブラリとデータを読み込む
library(ChannelAttribution)
library(reshape2)
library(ggplot2)
data(PathData)

#ヒューリスティックモデルを推定
H=heuristic_models(Data,'path','total_conversions',var_value='total_conversion_value')

#マルコフ連鎖モデルを推定
M=markov_model(Data, 'path', 'total_conversions', var_value='total_conversion_value')

#合計コンバージョンをプロット
R=merge(H,M,by='channel_name')
R1=R[,(colnames(R)%in%c("channel_name","first_touch_conversions","last_touch_conversions"
,"linear_touch_conversions","total_conversion"))]
colnames(R1)=c('channel_name','first_touch','last_touch','linear_touch','markov_model')
R1=melt(R1,id="channel_name")
ggplot(R1, aes(channel_name, value, fill = variable)) +
 geom_bar(stat="identity", position = "dodge") +
 ggtitle("TOTAL CONVERSIONS")+
 theme(axis.title.x = element_text(vjust=-2))+
 theme(axis.title.y = element_text(vjust=+2))+
 theme(title = element_text(vjust=2))+
 theme(text = element_text(size=16)) +
 theme(plot.title=element_text(size=20)) +
 ylab("")

#合計価値をプロット
R2=R[,(colnames(R)%in%c("channel_name","first_touch_value","last_touch_value","linear_tou
ch_value","total_conversion_value"))]
colnames(R2)=c('channel_name','first_touch','last_touch','linear_touch','markov_model')
R2=melt(R2,id="channel_name")
ggplot(R2, aes(channel_name, value, fill = variable)) +
 geom_bar(stat="identity", position = "dodge") +
 ggtitle("TOTAL VALUE")+
 theme(axis.title.x = element_text(vjust=-2))+
 theme(axis.title.y = element_text(vjust=+2))+
 theme(title = element_text(vjust=2))+
 theme(text = element_text(size=16)) +
 theme(plot.title=element_text(size=20)) +
 ylab("")

#ライブラリとデータを読み込む

library(ChannelAttribution)

library(reshape2)

library(ggplot2)

data(PathData)

#ヒューリスティックモデルを推定

H=heuristic_models(Data,'path','total_conversions',var_value='total_conversion_value')

#マルコフ連鎖モデルを推定

M=markov_model(Data, 'path', 'total_conversions', var_value='total_conversion_value')

#合計コンバージョンをプロット

R=merge(H,M,by='channel_name')

R1=R[,(colnames(R)%in%c("channel_name","first_touch_conversions","last_touch_conversions"

,"linear_touch_conversions","total_conversion"))]

colnames(R1)=c('channel_name','first_touch','last_touch','linear_touch','markov_model')

R1=melt(R1,id="channel_name")

ggplot(R1, aes(channel_name, value, fill = variable)) +

geom_bar(stat="identity", position = "dodge") +

ggtitle("TOTAL CONVERSIONS")+

theme(axis.title.x = element_text(vjust=-2))+

theme(axis.title.y = element_text(vjust=+2))+

theme(title = element_text(vjust=2))+

theme(text = element_text(size=16)) +

theme(plot.title=element_text(size=20)) +

ylab("")

#合計価値をプロット

R2=R[,(colnames(R)%in%c("channel_name","first_touch_value","last_touch_value","linear_tou

ch_value","total_conversion_value"))]

colnames(R2)=c('channel_name','first_touch','last_touch','linear_touch','markov_model')

R2=melt(R2,id="channel_name")

ggplot(R2, aes(channel_name, value, fill = variable)) +

geom_bar(stat="identity", position = "dodge") +

ggtitle("TOTAL VALUE")+

theme(axis.title.x = element_text(vjust=-2))+

theme(axis.title.y = element_text(vjust=+2))+

theme(title = element_text(vjust=2))+

theme(text = element_text(size=16)) +

theme(plot.title=element_text(size=20)) +

ylab("")

実行結果

ぱっと見ですが、alphaやtheta以外では均等配分モデルとマルコフ連鎖モデルに大きな差はなさそうです。データにもよるとは思いますが。自社のデータで実行して、違いがあるか見てみたいですね。

こちらはコンバージョン数で評価したもの

こちらはコンバージョン価値で評価したもの

参考文献

Markov model for the online multichannel attribution problem

アトリビューションモデルの概要

マルコフ連鎖モデルによるアトリビューション分析（インターネット広告への応用）

Attribution models—数理モデリング入門

マルコフ連鎖を説明してみる。

LDA（潜在的ディリクレ配分法）まとめ手法の概要と試行まで

【目次】
・トピックモデルとは
・トピックモデルの歴史
・トピックモデルでできること
・トピックモデルを理解するために必要な知識
・トピックモデルの手法について
・トピックモデルの実行方法について（R言語）
・トピックモデルの評価方法について
・Correlated Topic Models (CTM)について
・PAM:Pachinko Allocation Modelについて
・Relational Topic Models(RTM)について
・参考文献

トピックモデルとは

・一つの文書に複数の潜在的なトピック（話題・分野・文体・著者など）、文書ごとのトピック分布、トピックごとの単語分布を仮定して、テキストからトピックを推定するモデル。文書に限らず、様々な離散データに隠された潜在的なトピックを推定するベイジアンモデル。幅広いドメインでの離散データで有効とされている。

トピックモデルの歴史

1998年：pLSA(probabilistic Latent Semantic Analysis)
2003年：LDA(Latent Dirichlet Allocation)
2004年〜：拡張モデル
2007年〜：大規模データのための高速化

トピックモデルでできること

・人を介することなく、大量の文書集合から話題になっているトピックを抽出できる。
・文書データだけでなく、画像処理、推薦システム、ソーシャルネットワーク解析、バイオインフォマティクス、音楽情報処理などにも応用されている。
・確率過程を用いて、ノイズを取り除き、本質的な情報を抽出できる。

トピックモデルを理解するために必要な知識

・確率
　-確率分布
　　-ベルヌーイ分布
　　-カテゴリ分布
　　-ベータ分布・ガンマ分布
　　-ディリクレ分布
・ラグランジュ未定乗数法
・ユニグラム
　-BOW(Bag of words)
　-ユニグラムモデル
・混合ユニグラムモデル
・混合モデル
・EMアルゴリズム
・最尤推定
・ベイズの定理
　-ベイズ推定
　 -ベイズ予測分布
　 -ハイパーパラメータ推定
　 -モデル選択
　 -変分ベイズ推定
　 -ギブスサンプリング

トピックモデルの手法について

推定方法としては以下の三つが提案されている。
・最尤推定
・変分ベイズ推定
・ギブスサンプリング

ギブスサンプリングによる方法

一部のパラメータを積分消去して、トピック集合の事後分布を推定
↓
文書集合とトピック集合の同時分布を導出
↓
単語ごとにトピックをサンプリングする
↓
サンプリングされたトピックからトピック分布と単語分布を推定
↓
周辺同時尤度を最大化させるハイパーパラメータαとβを推定する

※LDAのギブスサンプリングはLDAが共役性に基づいたモデリングであるため効率的な計算が可能とされる。

ハイパーパラメータを適切にデータから学習すれば、MAP推定・変分ベイズ推定・ギブスサンプリングの性能の差は大きく出ないという研究結果があるらしい。（なお、MAP推定だと、クロスバリデーションを行い、尤度を最大化させるハイパーパラメータを見つけなければならない。）

トピックモデルの実行方法について（R言語）

以下のパッケージで実行可能だが、新しい手法などにはまだ対応していない。
lda（CRANのPDF）
topicmodels（CRANのPDF）

以下はtopicmodelsパッケージの実行コードであるが、BOW形式のデータがあれば、実行はすぐにできる。
ただし、パープレキシティなどを計算しながら、ハイパーパラメータのチューニングをする必要がある。

library(topicmodels)
k <- 10 #トピック数
LDA_estimate <- LDA(bagofwords, k, method="Gibbs",control=list(alpha=alpha,verbose=1,
                                                                            iter=10000,burnin=1000,delta=delta))

library(topicmodels)

k <- 10 #トピック数

LDA_estimate <- LDA(bagofwords, k, method="Gibbs",control=list(alpha=alpha,verbose=1,

iter=10000,burnin=1000,delta=delta))

時間があれば、ソースコードを見て自分で書けるようにしたい。
ちなみに、HDP-LDAはPythonのgensimに用意されているようです。(gensimへのリンク)

トピックモデルの評価方法について

パープレキシティ（Perplexity）

-確率モデルの性能を評価する尺度として、テストデータを用いて計算する。
-負の対数尤度から計算できる。
-低いパープレキシティはテストデータを高い精度で予測できる良い確率モデルであることを示す。
-全ての語彙が一葉の確率で出現するモデルのパープレキシティは語彙数Vになり、性能が悪いと考えることができる。
-このパープレキシティが小さくなるようなトピック数やパラメータの設定が必要。

Correlated Topic Models (CTM)について

トピックモデルは「各トピック k の間には相関がない」という前提のもとに成り立つ手法ですが、本当にトピック間に相関はないのでしょうか。「本当は相関のあるトピック」を無理やり「相関のないトピック」に分割している可能性が高いので、相関を加味したモデルが必要です。そこで、トピックの生成割合を決める際に、トピック間の相関を持つように多次元の正規分布を用います。ただし、その代わりに共役的な分布を持たないため、従来のギブスサンプリングではサンプリングが非効率になることから変分ベイズを用いる必要があります。

PAM:Pachinko Allocation Modelについて

CTMは各トピックが同じレベルにあるため、トピック間の階層構造を表すことができません。
自動車にとっての、セダンやトラック、あるいはお酒にとってのワインやビールのようにトピックに関しても階層構造があることは想像に難くありません。トピック間の関係・相関を一般的に表現するモデルで、トピック間の階層構造を基本として、パチンコ玉が落ちるように単語を生成します。パチンコ玉が落ちるというのは、向きが一方行の有向かつ非巡回ノードが連想されます。分布に関しても共役なので、ギブスサンプリングでサンプリングが可能です。
2016/04/24の段階ではまだGensimでモデルの開発はなされていないようです。
Pachinko Allocation Model
ちなみに、論文はこちらです。
Pachinko Allocation: DAG-Structured Mixture Models of Topic Correlations

Relational Topic Models(RTM)について

文書の中身だけでなく、文書間のリンクの生成過程も同時に確率モデル化する手法。論文や特許データに対して活用できる。過去の購買行動に応じた、顧客のセグメント解析や商品のレコメンデーションに活用できるかもしれない。

参考文献

『トピックモデル (機械学習プロフェッショナルシリーズ)』
『トピックモデルによる統計的潜在意味解析 (自然言語処理シリーズ)』
Tokyo Webmining 46th 『トピックモデルことはじめ』
machine_learning_python/topic.md

統計数理研究所 H24年度公開講座「確率的トピックモデル」サポートページ

Rのleafletパッケージを使って蒙古タンメン中本をプロットしてみた

Rのleafletパッケージで蒙古タンメン中本の店舗をマッピングしました。
中本の全店舗地図は未だ観察できていないので、中本信者の方には嬉しい情報かもしれません。
元となるデータは中本のホームページから収集し、GoogleMapのAPIで住所から緯度・経度を求めました。

Rに読み込ませるデータはこんな感じになりました。

コードはこんな感じです。如何に簡単かがおわかりいただけると思います。

library(leaflet)
library(htmlwidgets)
setwd("任意")

nakamoto <- read.csv("nakamoto.csv",header=T)

m <- leaflet(nakamoto) %>%
  addTiles() %>%
  addMarkers(lng=~longitude,lat=~latitude) %>%
  addPopups(lng=~longitude,lat=~latitude,popup=~store_name)
m

saveWidget(m, file="m.html")

library(leaflet)

library(htmlwidgets)

setwd("任意")

nakamoto <- read.csv("nakamoto.csv",header=T)

m <- leaflet(nakamoto) %>%

addTiles() %>%

addMarkers(lng=~longitude,lat=~latitude) %>%

addPopups(lng=~longitude,lat=~latitude,popup=~store_name)

saveWidget(m, file="m.html")

参考情報

embed R leaflet map in wordpress

Japan.R 2015に参加して得た情報まとめ

Japan.Rで私にとって使えそうなアイデアや手法だけを厳選して載せたいと思います。（結局、長くなった気がしますが。。）
全てを見たい場合は、この方のブログがよいと思います。
Japan.R 2015に参加してきた

Spark Streaming

Spark Streamingストリーミングを使って、挙動に応じて自動的にメールやSlackに送るとかでサイトの監視に使えるそうな。

jupyter

RでいうところのRStudioのような分析しやすいPython用のGUI
　

rangerパッケージ

RFなどができる。C++で実装されているので速いらしい。
最近のRのランダムフォレストパッケージ -ranger/Rborist-

RcppRollパッケージ

移動平均とかの計算が楽になるらしい。
インターネットの予約のステップの伸び率をモデリングし、ある一定速度だとしたら表示を出し分けるなどの施策に活かしているそうな。

多腕バンディットアルゴリズム

サイト内のバナーなどの出し分けの最適化に使っているらしい。Rでは重くてできないそうな。

leafletパッケージ

これを使えば地理データの可視化が簡単にできる。
→デフォルトの地図がOpenStreetMapなので、サイトに載せる際も好都合。
→「ぐりぐり動かせる・レイヤー構造・多機能・Shinyでアプリケーションも作れる」
→Githubの方がバージョン上
→データテーブルで読み込んで緯度経度を指定する。
　サークルや線の太さも指定でき、ポップアップも立てれる。
　地点間でのラインも引ける。addPolylines
　プロットでエクスポートでWebページにすれば図が作れる。
leafletでMapbox（とかいろいろな地図）を使う

purrrパッケージ

→複数のモデルをリストにして、推定結果を一気に出すことができる。（都度、モデルを作る必要がない。）
{purrr} による非テーブルデータの処理
 purrr 0.1.0を使ってみる(1) map

ggfortifyパッケージ

purrrで実行した分析結果のプロットなどが簡単にできるらしい。
Concepts and Basics of ggfortify

Shiny

分析結果を可視化したりできるツールとしてはカジュアルにできていいかもしれない。
無料版では複数のアクセスに弱いので社内用のツールとしてはいい。

rApacheパッケージ

rとアパッチでwebアプリケーションを作れる。

deployR

APIを作るためのフレームワーク

Gradient Boosting Decision Tree

モデルの性能が良くなるように特徴量の選択ができる手法らしいです。これは進捗が良くなる、イケメンな手法ですね。
今度深掘ってみたいです。
Gradient Boosting Decision Treeでの特徴選択 in R

SparkRパッケージ

Rからスパークの処理を実行できる。
SparkRをEC2上で動かして分散処理してみる

Parquet（パーケ）

ビッグデータで一般的になってきているフォーマットのことらしい。
10GBのデータが300MBになった。→これはすごい！！
DMMゲームのログ解析~Parquetについて~

TinySegmenter

MeCabとかを入れてごにょごにょすることが許されていない方向け？
TinySegmenter:Javascriptだけで書かれたコンパクトな分かち書きソフトウェア

knitrパッケージ

Rの分析結果をhtmlでレポーティングができる。
R でレポートを作成するのに knitr が超便利

R Presentation

Studioでスライド作成できる。
Rとウェブの融合（５）ーslidify、R Presentationー

rvestパッケージ

RでWebスクレイピングができてしまう。（私はPythonでやっていますが。）
rvest でログインしてスクレイピング #rstatsj

pforeachパッケージ

デフォルトで並列計算するようにしたパッケージ
R で超簡単に並列処理を書けるパッケージ pforeach を作った

RstanでGoogleトレンドの推定

rstanを用いて、Googleトレンドデータの予測モデルを推定してみます。

ほとんど岩波データサイエンスのものですが、Googleトレンドのデータを月ごとの季節性を加味した状態空間モデルを用いて予測してみました。

今回の分析では、
・modelのstanコード(stan)
・Rでstanを動かすためのコード(R)
（・可視化のためのコード(R)）←必須ではない
を用意します。

データですが、GoogleTrendのサイトで任意のキーワードで検索して、
その時系列データをCSVでダウンロードすれば手に入ります。（ちょっと見つけにくい）

データの形式はシンプルで、
先頭にY
とおいて後はトレンドの値を行ごとに置いていけばいけます。

要はN行1列データをテキストファイルに保存すればOKです。（1行目はY）

まずstanのコードですが、岩波データサイエンスのサンプルコードの季節を4から12に変えています。（たったこれだけ）
Googleトレンドのデータは月単位でも結構値がふれることがあるので、月ごとに応じた潜在的な変数が必要だと思いました。

data {
   int<lower=1> T;
   int<lower=1> T_next;
   real Y[T];
}

parameters {
   real mu[T];
   real s[T];
   real<lower=0> s_mu;
   real<lower=0> s_s;
   real<lower=0> s_r;
}

model {
   for(t in 2:T)
      mu[t] ~ normal(mu[t-1], s_mu);
   for(t in 12:T)
      s[t] ~ normal(-s[t-1]-s[t-2]-s[t-3]-s[t-4]-s[t-5]-s[t-6]-s[t-7]-s[t-8]-s[t-9]-s[t-10]-s[t-11], s_s);
   for(t in 1:T)
      Y[t] ~ normal(mu[t]+s[t], s_r);
}

generated quantities {
   real mu_all[T+T_next];
   real s_all[T+T_next];
   real y_next[T_next];

   for (t in 1:T){
      mu_all[t] <- mu[t];
      s_all[t] <- s[t];
   }
   for (t in (T+1):(T+T_next)){
      mu_all[t] <- normal_rng(mu_all[t-1], s_mu);
      s_all[t] <- normal_rng(-s_all[t-1]-s_all[t-2]-s_all[t-3]-s_all[t-4]-s_all[t-5]-s_all[t-6]-s_all[t-7]-s_all[t-8]-s_all[t-9]-s_all[t-10]-s_all[t-11], s_s);
   }
   for (t in 1:T_next)
      y_next[t] <- normal_rng(mu_all[T+t]+s_all[T+t], s_r);
}

data {

int<lower=1> T;

int<lower=1> T_next;

real Y[T];

}

parameters {

real mu[T];

real s[T];

real<lower=0> s_mu;

real<lower=0> s_s;

real<lower=0> s_r;

}

model {

for(t in 2:T)

mu[t] ~ normal(mu[t-1], s_mu);

for(t in 12:T)

s[t] ~ normal(-s[t-1]-s[t-2]-s[t-3]-s[t-4]-s[t-5]-s[t-6]-s[t-7]-s[t-8]-s[t-9]-s[t-10]-s[t-11], s_s);

for(t in 1:T)

Y[t] ~ normal(mu[t]+s[t], s_r);

}

generated quantities {

real mu_all[T+T_next];

real s_all[T+T_next];

real y_next[T_next];

for (t in 1:T){

mu_all[t] <- mu[t];

s_all[t] <- s[t];

}

for (t in (T+1):(T+T_next)){

mu_all[t] <- normal_rng(mu_all[t-1], s_mu);

s_all[t] <- normal_rng(-s_all[t-1]-s_all[t-2]-s_all[t-3]-s_all[t-4]-s_all[t-5]-s_all[t-6]-s_all[t-7]-s_all[t-8]-s_all[t-9]-s_all[t-10]-s_all[t-11], s_s);

}

for (t in 1:T_next)

y_next[t] <- normal_rng(mu_all[T+t]+s_all[T+t], s_r);

}

Rでstanを動かすためのコードですが、ここはサンプルコードとほぼ一緒です。

library(rstan)
setwd("任意")

d <- read.csv('data-trend4.txt', header=TRUE)
T <- nrow(d)
T_next <- 8
data <- list(T=T, T_next=T_next, Y=d$Y)

stanmodel <- stan_model(file='model_custom.stan')
fit <- sampling(
   stanmodel, data=data, pars=c('mu_all','s_all','y_next','s_mu','s_s','s_r'),
   iter=10200, warmup=200, thin=10, chains=3,
   seed=123
)

library(rstan)

setwd("任意")

d <- read.csv('data-trend4.txt', header=TRUE)

T <- nrow(d)

T_next <- 8

data <- list(T=T, T_next=T_next, Y=d$Y)

stanmodel <- stan_model(file='model_custom.stan')

fit <- sampling(

stanmodel, data=data, pars=c('mu_all','s_all','y_next','s_mu','s_s','s_r'),

iter=10200, warmup=200, thin=10, chains=3,

seed=123

)

可視化のためのコードについてもサンプルコードとほぼ一緒です。

library(ggplot2)

# after estimation

d_obs <- data.frame(X=1:T, Y=d$Y)
p <- ggplot()
p <- p + theme_bw() + theme(text=element_text(size=18))
p <- p + geom_line(data=d_obs, aes(x=X, y=Y), color='black', alpha=0.8, size=2)
p <- p + labs(x='Time [month]', y='trend')
p <- p + coord_cartesian(xlim=c(0.9, 152.1))
ggsave(file='fig2-top-left.png', plot=p, dpi=300, width=6, height=4)


makeDataFrameQuantile <- function(x, y_smp){
  qua <- apply(y_smp, 2, quantile, prob=c(0.1, 0.25, 0.5, 0.75, 0.9))
  d_est <- data.frame(X=x, t(qua))
  colnames(d_est) <- c('X', 'p10', 'p25', 'p50', 'p75', 'p90')
  return(d_est)
}

plotTimecourse <- function(file, d_est, d_obs){
  p <- ggplot()
  p <- p + theme_bw() + theme(text=element_text(size=18))
  p <- p + geom_vline(xintercept=T, linetype='dashed')
  p <- p + geom_ribbon(data=d_est, aes(x=X, ymin=p10, ymax=p90), fill='black', alpha=0.25)
  p <- p + geom_ribbon(data=d_est, aes(x=X, ymin=p25, ymax=p75), fill='black', alpha=0.5)
  p <- p + geom_line(data=d_est, aes(x=X, y=p10), color='black', size=0.2)
  p <- p + geom_line(data=d_est, aes(x=X, y=p90), color='black', size=0.2)
  p <- p + geom_line(data=d_est, aes(x=X, y=p25), color='black', size=0.2)
  p <- p + geom_line(data=d_est, aes(x=X, y=p75), color='black', size=0.2)
  p <- p + geom_line(data=d_est, aes(x=X, y=p50), color='black', size=0.4)
  if (!is.null(d_obs)){
    p <- p + geom_line(data=d_obs, aes(x=X, y=Y), color='black', size=2, alpha=0.9)
  }
  p <- p + labs(x='month', y='trend')
  p <- p + coord_cartesian(xlim=c(0.9, 152.1))
  ggsave(file=file, plot=p, dpi=300, width=6, height=4)
}

la <- rstan::extract(fit)
d_est <- makeDataFrameQuantile(x=1:(T+T_next), y_smp=la$mu_all)
plotTimecourse(file='fig2-bottom-left.png', d_est=d_est, d_obs=d_obs)

d_est <- makeDataFrameQuantile(x=1:(T+T_next), y_smp=la$s_all)
plotTimecourse(file='fig2-bottom-right.png', d_est=d_est, d_obs=NULL)

d_est <- makeDataFrameQuantile(x=(T+1):(T+T_next), y_smp=la$y_next)
d_est <- rbind(data.frame(X=T, p10=d$Y[T], p25=d$Y[T], p50=d$Y[T], p75=d$Y[T], p90=d$Y[T]), d_est)
plotTimecourse(file='fig2-top-right.png', d_est=d_est, d_obs=d_obs)

library(ggplot2)

# after estimation

d_obs <- data.frame(X=1:T, Y=d$Y)

p <- ggplot()

p <- p + theme_bw() + theme(text=element_text(size=18))

p <- p + geom_line(data=d_obs, aes(x=X, y=Y), color='black', alpha=0.8, size=2)

p <- p + labs(x='Time [month]', y='trend')

p <- p + coord_cartesian(xlim=c(0.9, 152.1))

ggsave(file='fig2-top-left.png', plot=p, dpi=300, width=6, height=4)

makeDataFrameQuantile <- function(x, y_smp){

qua <- apply(y_smp, 2, quantile, prob=c(0.1, 0.25, 0.5, 0.75, 0.9))

d_est <- data.frame(X=x, t(qua))

colnames(d_est) <- c('X', 'p10', 'p25', 'p50', 'p75', 'p90')

return(d_est)

}

plotTimecourse <- function(file, d_est, d_obs){

p <- ggplot()

p <- p + theme_bw() + theme(text=element_text(size=18))

p <- p + geom_vline(xintercept=T, linetype='dashed')

p <- p + geom_ribbon(data=d_est, aes(x=X, ymin=p10, ymax=p90), fill='black', alpha=0.25)

p <- p + geom_ribbon(data=d_est, aes(x=X, ymin=p25, ymax=p75), fill='black', alpha=0.5)

p <- p + geom_line(data=d_est, aes(x=X, y=p10), color='black', size=0.2)

p <- p + geom_line(data=d_est, aes(x=X, y=p90), color='black', size=0.2)

p <- p + geom_line(data=d_est, aes(x=X, y=p25), color='black', size=0.2)

p <- p + geom_line(data=d_est, aes(x=X, y=p75), color='black', size=0.2)

p <- p + geom_line(data=d_est, aes(x=X, y=p50), color='black', size=0.4)

if (!is.null(d_obs)){

p <- p + geom_line(data=d_obs, aes(x=X, y=Y), color='black', size=2, alpha=0.9)

}

p <- p + labs(x='month', y='trend')

p <- p + coord_cartesian(xlim=c(0.9, 152.1))

ggsave(file=file, plot=p, dpi=300, width=6, height=4)

}

la <- rstan::extract(fit)

d_est <- makeDataFrameQuantile(x=1:(T+T_next), y_smp=la$mu_all)

plotTimecourse(file='fig2-bottom-left.png', d_est=d_est, d_obs=d_obs)

d_est <- makeDataFrameQuantile(x=1:(T+T_next), y_smp=la$s_all)

plotTimecourse(file='fig2-bottom-right.png', d_est=d_est, d_obs=NULL)

d_est <- makeDataFrameQuantile(x=(T+1):(T+T_next), y_smp=la$y_next)

d_est <- rbind(data.frame(X=T, p10=d$Y[T], p25=d$Y[T], p50=d$Y[T], p75=d$Y[T], p90=d$Y[T]), d_est)

plotTimecourse(file='fig2-top-right.png', d_est=d_est, d_obs=d_obs)

以上を実行した結果、以下のような図が出てきます。

こちらは実際の時系列データのプロットです。

8期先までの予測です。

8期先までの予測範囲です。信頼区間90%までの範囲となっています。

推定した潜在的な季節性のデータをプロットしています。

ついでに、4月までのデータを用いて、5~9月の予測を行い、その比較を行っています。

5月が大きく外れましたが、その後はある程度当てれているように見えます。
5月も当てれるようなモデルを作りたいものですね。

参考文献

<br />

Sparkの仕組み・導入メリット・インストール方法など

【目次】
・Sparkについて
・Sparkの歴史
・Sparkの仕組み
・データ分析におけるSpark導入のメリット
・Sparkの導入方法
・参考文献

Sparkについて

OSS(Open Source Software)のインメモリ処理の分散並列基盤
※インメモリ(メインメモリ上で処理するため、ディスクストレージを介するよりも高速かつ安定しているらしい。メモリに乗り切らないようなケースでも、ディスクを利用するなどしてアプリケーションが問題なく動作するように作られている。)
ビッグデータ活用を支えるプラットフォームとして注目されている。
比較の際はHadoopで言う所の、MapReduceが引き合いに出されます。MapReduceがバッチ処理に特化した仕組みになっているのに対して、Sparkはバッチ処理だけでなく、対話的な処理や繰り返し処理・ストリーミング処理といった多様なデータ処理のパターンに対応できる仕組みになっていることから、Sparkに注目が集まっています。

分散処理エンジンとして、Apache Spark Coreがあり、以下の様な機能が用意されているようです。

Spark SQL・・・SQL処理用
Spark Streaming・・・ストリーミング処理
MLib・・・機械学習処理用
　基本API
　・統計情報の取得
　・相関関係の算出
　・層化抽出
　・仮説検定
　・ランダムデータ生成
　機械学習アルゴリズム
　・分類回帰（SVM、ロジット、線形回帰、ナイーブベイズ、決定木、ランダムフォレスト、isotonic回帰）
　・協調フィルタリング（ALS(Alternating Least Squares)）
　・クラスタリング（K平均法、混合ガウシアン、PIC、LDA、ストリーミングK平均法）
　・次元削減（SVD、PCA）
　・特徴抽出&変換（TF-IDF、Word2Vec、etc…）
　・頻出パターンマイニング（FP-Groth、アソシエーションルール、PrefixSpan）
GraphX・・・グラフの処理
※機械学習の機能があるのは嬉しいですね。スクラッチで作るのはペインフルなので。

データソースとしては、HDFS(Hadoop Distributed File System)だけでなく、Hive、HBase、PodgreSQL、MySQL、CSVファイルなども対応している。

Sparkの歴史

2009年にUC BarkeleyのMatei Zaharia氏がScalaを用いて開発

Google Trendで検索トレンドを見たところ、2014年頃から今に至るまで伸び続けています。

Sparkの仕組み

RDD(Resilient Distributed Dataset)
不変で並列実行可能な分割されたコレクションのことを指します。
これによりユーザーが並列分散処理を意識しなくとも、それができる環境が提供されます。
こちらのRDDを
map,filter,union,flatMapなどのTransformations
で生み出し、RDDのデータを
reduce,countなどのActions
を実行することで動かします。
RDDを任意の時点でメモリやディスクに保存することもできるため、処理の再実行のみならず、同じデータを複数回利用するなど、複雑なデータ処理に関しても効率的に実行できます。
SparkはScalaという言語で書かれていますが、Python・Java・Rなどに対してAPIが提供されているので、幅広いユーザーが扱うことが可能です。しかしながら、痒いところに手を届かせるためには、Scalaの勉強は避けられないのかもしれません。実際のところ、大体のデータマイニング系エンジニアの知り合いはScalaの勉強などをしていた印象があります。

データ分析におけるSpark導入のメリット

・並列分散の手法を意識せずに大規模データの処理が可能になる。
・Rなどの苦手とする反復処理を高速化させることができる。
・データ取得、前処理、分析処理、レポーティングを一つのプラットフォームで実行できる。
・コマンドラインによるインタラクティブな試行錯誤がしやすいので、データ分析のタスクに向いている。(PythonやRのように)

※レコメンドエンジンとしてSparkを用いるケースも観察されているようです。
　Yahoo!台湾

Sparkの導入方法

ここでの導入はあくまでもプログラミングを学ぶためであって、ローカル環境での実行環境の構築ではSparkのメリットを感じることはできません。
※SparkはJava仮想マシン上で動作するため、事前にOracle JDKをインストールしておきましょう。
Java SE Development Kit 8 Downloads

■Windowsの場合
Apache Sparkで始めるお手軽リアルタイムウインドウ集計
こちらで手順が詳しく書かれています。

■Macの場合
Macの方が導入はシンプルです。
Download Spark
からダウンロードします。（SparkRはSparkのver1.4以降からサポートされているそうです。）
HOW TO INSTALL APACHE SPARK ON MAC OS X YOSEMITE
Apache Spark 1.6.0 setup on Mac OS X Yosemite
そこで、

build/mvn -Pyarn -Phadoop-2.4 -Dhadoop.version=2.4.0 -DskipTests clean package

を実行すればコンパイルできるようです。（10分ほどかかる。私は30分もかかったが。。）

そのあと、パスを通します。
PATH=/usr/local/share/spark/bin:$PATH; export PATH

それから、pythonでsparkを動かしてみます。
pyspark

を実行して、試しにReadmeのテキスト行数をカウントさせるコマンドを入力します。これは、textFileメソッドを用いて、Readmeファイルを読み込み、linesという名前のRDDを生成していることを意味します。
lines = sc.textFile(“/usr/local/share/spark/README.md”)

そして、このコマンドで値が返されたら、うまく動いていることが確認できるでしょう。
print lines.count()

加えて、ワードカウントを行うコマンドを入力します。（スペース区切りで分解し、mapメソッドで（単語, 1）という組み合わせに変換したあと、最終的にreduceByKeyで同一の単語を足し上げるという処理をしています。）

counts = lines.flatMap(lambda line: line.split(” “)) \
.map(lambda word: (word, 1)) \
.reduceByKey(lambda a, b: a + b)

そして、これで結果を取得します。
counts.collect()