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

廿TT

譬如水怙牛過窓櫺 頭角四蹄都過了 因甚麼尾巴過不得

[RStan]同時確率に基づく検索キーワードのクラスタリング(改訂版)

[RStan]同時確率に基づく検索キーワードのクラスタリング - 廿TT の改訂版です。

上エントリはそもそもやろうとしてることに無理があったと思い、ある程度クラスタを人力で与えることにしました。

モデル

SEO 的には本当は語順にも意味があるんだろうけれども、ここでは気にしないことにする。

ある単語がクエリに登場したかのみを見る。

単語はぜんぶで M 種類あるとする。

検索クエリは L 種類あるとする。

単語はそれぞれ独立に、クエリに登場する確率 p[m] を持っているとする。

たとえば M=3 で考えると、あるクエリ「A B C」が出て来る確率は p[1] p[2] p[3]、クエリ「A B」が出て来る確率は p[1] p[2] (1-p[3]) になる。

インプレッションの回数やクリック数も考慮したほうがいいのかもしれないけれども、今回は無視。

この分布が混合比率 \rm \phi_k で K 個混ざっているとする。

クエリの背後には検索者のモチベーション(トピック)が K 種類あって、それによって出やすい単語、出にくい単語が変わるという解釈です。

w[l,m] をクエリ l に単語 m が登場したとき 1、さもなくば 0 の値をとる変数とする。

いくつかのワードにはトピックのラベルが与えられているとして、対数尤度は、

\rm log(\phi[label[l]])

と、

\rm \sum_{m=1}^{M} \{w[l,m] \log p[label[l],m]+\\~~ (1-w[l,m]) \log(1-p[label[l],m])\})

の和である。

トピックが与えられなかったクエリについては、k について足し合わせて、周辺分布にして対数尤度を書くと、

 \rm \log  \sum_{k=1}^{K} \exp(\log \phi_k + \\~~\sum_{j=1}^{M} \{w[l,m] \log p[k,m]+\\~~~~~~ (1-w[l,m]) \log(1-p[k,m])\})

となる。

ここで l はクエリの数だけ動く( \rm l=1,\ldots,L)。

さて、確率関数を \rm f_k、観測値を y とまとめて書くと, 混合分布を考えたとき、観測 y が得られたという条件のもとで、その観測がトピック k から得られたという確率は、ベイズの定理より、

 \rm \phi_k f_k (y;p)/\left( \sum_{k=1}^{K} \phi_k f_k (y;p)\right)

となる。

属する確率が最大になるトピックをそのクエリのトピックとすることで、クラスタリングができる。

Stan のコード

w[i,j] と 1-w[i,j] を直接 Stan の上で計算すると for 文で回さないといけないので、w[i,j] を要素とする行列 W1 と 1-w[i,j] を要素とする行列 W2 をデータとして渡すことする。

generated quantities ブロックがベイズの定理でクラスタリングするところです。

//multi_mix.stan
data{
  int<lower=1> K; #num of topics
  int<lower=1> L; #num of queries
  int<lower=1> M; #num of words
  int<lower=0,upper=M> J; #num of labeled words
  matrix<lower=0,upper=1>[L,M-J] W1;
  matrix<lower=0,upper=1>[L,M-J] W2;
  int<lower=0,upper=K> label[L];
}
parameters{
  simplex[K] phi;
  vector<lower=0,upper=1>[M-J] p[K];
}
model{
  for(l in 1:L){
    if(label[l]!=0){
      target += log(phi[label[l]]);
      target += W1[l,] * log(p[label[l]]) + W2 * log1m(p[label[l]]);
    }else{
      vector[K] gamma;
      gamma = log(phi);
      for(k in 1:K){
       gamma[k] = gamma[k] + W1[l,] * log(p[k]) + W2[l,] * log1m(p[k]);
      }
      target += log_sum_exp(gamma);
    }
  }
}
generated quantities{
  matrix[L,K] labels_mat;
  vector[K] tmp;
  vector[K] gamma;
  for(l in 1:L){
  if(label[l]!=0){
    for(k in 1:K){
      tmp = rep_vector(0,K);
      tmp[label[l]] =1;
      labels_mat[l,] =tmp';
      }
    }else{
    gamma = log(phi);
    for(k in 1:K){
       gamma[k] = gamma[k] + W1[l,] * log(p[k]) + W2[l,] * log1m(p[k]);
    }
      labels_mat[l,] = softmax(gamma)';
    }
  }
}

結果

ページ http://abrahamcow.hatenablog.com/entry/2014/12/19/221230 が表示されたクエリのうち、クリック数が 1 以上あるクエリを分析対象とした。

ユニークなクエリの数 L は 85、単語の数 M は 38 だった。

クラスタの数 K は 2 とした。

勘と経験にもとづき、"エクセル" または "excel" を含むクエリは、必ずクラスタ 1 に分類することにした。

4 chain 回して、トレースプロットの目視と Rhat 指標で収束を確認した。

f:id:abrahamcow:20170501021840p:plain

2000回のイテレーションが 1 chain あたり 1 分くらいで終わった。

クラスタごとに p[m] の推定値を見る。下図は事後期待値の棒グラフと 95% ベイズ信頼区間のエラーバーです。

f:id:abrahamcow:20170501021945p:plain

"エクセル" または "excel" を含むまないクラスタ 2 は、「とは」や「意味」の出現確率が高い。

また、"エクセル" または "excel" と親和性が高いクラスタ 1 では、「グラフ」や「反比例」の出現確率が高いことがわかる。

クラスタ 1 はエクセルの近似曲線機能について知りたい層、クラスタ 2 はそもそも近似曲線とはなにかに興味がある層と解釈できる。

検索アナリティクスの指標をクラスタごとに見るために下の箱ひげ図を書いた。

f:id:abrahamcow:20170501022731p:plain

大きな差はつかなかったがクリック、インプレッション、CTR(クリックスルーレート)、ポジション(ポジションは数字が小さいほうが上)いずれもクラスタ 1 が若干高めで、やっぱりエクセルって人気なんだなとわかる。

R のコード

処理の流れはこんな感じ。

  1. searchConsoleR でデータを読み込む
  2. RMeCab で形態素解析
  3. 下処理
  4. Stan で推定
  5. ggplot で図示、解釈
#あとで使うパッケージをまとめてロードしておく
library(rstan)
library(RMeCab)
library(dplyr)
library(tidyr)
library(cowplot)
library(searchConsoleR)
library(mixtools)
#######
#サーチコンソールでは 90 日分のデータが保持されるようなので、90 日分とってくる
scr_auth()
sc_websites <- list_websites()
start <- as.character(Sys.Date() - 93)
end <- as.character(Sys.Date()-3)
scdata <- search_analytics(sc_websites[1,1], 
                           startDate = start, endDate = end,
                           dimensions = c("page","query"),
                           rowLimit = 20000)
#head(scdata)
#######
#分析したいランディングページを決める
url1 <- "http://abrahamcow.hatenablog.com/entry/2014/12/19/221230"
LP1 <-dplyr::filter(scdata,page==url1,clicks>0)
#形態素にバラす
query <-sapply(LP1$query,FUN=function(x)unlist(RMeCabC(x)))
#######
#下処理
UW <-unique(unlist(query))
M<-length(UW)
BoW <- t(sapply(query, function(x){as.integer(UW %in% x)}))
L=length(query)
M=length(UW)
flag <- (UW=="エクセル")|(UW=="excel")
W1 = BoW[,!flag]
W2 = 1-W1
q_labels <- as.integer(BoW[,which(UW == "エクセル")]>0|BoW[,which(UW == "excel")]>0)
J = sum(flag)
#######
#Stan の出番
multi_mix <- stan_model("multi_mix.stan")
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
K=2 #トピック数を決める
dat4stan <-list(K=K,J=J,M=M,L=L,W1=W1,W2=W2,label=q_labels)
fit1 <- sampling(multi_mix,dat4stan,iter=2000, chain=4)
traceplot(fit1,pars="phi")
all(summary(fit1)$summary[,"Rhat"]<1.1,na.rm = TRUE)
#######
#結果の取り出し
phat<-matrix(get_posterior_mean(fit1,"p"),K,M-J,byrow = TRUE)
get_posterior_mean(fit1,"phi")

df_mean <-as_data_frame(t(phat)) %>% 
  mutate(id=1:(M-J),words=UW[!flag]) %>% 
  gather(class,p,-id,-words)

df_upper<-as_data_frame(t(matrix(summary(fit1,pars="p")$summary[,"97.5%"],K,M-J,byrow = TRUE))) %>% 
  mutate(id=1:(M-J)) %>% 
  gather(class,upper,-id)

df_lower<-as_data_frame(t(matrix(summary(fit1,pars="p")$summary[,"2.5%"],K,M-J,byrow = TRUE))) %>% 
  mutate(id=1:(M-J)) %>% 
  gather(class,lower,-id)

df_p <- left_join(df_mean,df_lower,by=c("id","class")) %>% 
  left_join(df_upper,by=c("id","class"))

#棒グラフ
ggplot(df_p,aes(x=reorder(words,p),y=p))+
  geom_col(alpha=0.5)+
  geom_errorbar(aes(ymin=lower,ymax=upper),size=1,width=0.5,colour="salmon")+
  coord_flip()+
  facet_wrap(~class)+
  ylim(c(0,1))+xlab("")+
  theme_cowplot(font_family = "HiraKakuProN-W3",font_size = 12)

lab <-matrix(get_posterior_mean(fit1,"labels_mat"),L,K,byrow = TRUE)
LP1_g <- mutate(LP1,class=apply(lab,1,which.max)) %>% 
  select(-page) %>% 
  gather(metrics,value,-query,-class) %>% 
  mutate(class=factor(class))

ggplot(LP1_g,aes(x=class,y=value))+
  geom_boxplot()+
  facet_wrap(~metrics,scales = "free_y")+
  theme(legend.position = "none")

参考文献

StanとRでベイズ統計モデリング (Wonderful R)

StanとRでベイズ統計モデリング (Wonderful R)

実践 ベイズモデリング -解析技法と認知モデル-

実践 ベイズモデリング -解析技法と認知モデル-

関連エントリ

abrahamcow.hatenablog.com