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

廿TT

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

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

このモデルは推定するたびに結果が変わることが判明しました。申しわけありません。以下の記述はなかったことにしてください。

改訂版を書きました → [RStan]同時確率に基づく検索キーワードのクラスタリング 2 - 廿TT

はじめに

検索キーワードのグルーピングがしたい。

StanとRでベイズ統計モデリング (Wonderful R)実践 ベイズモデリング -解析技法と認知モデル- では Stan を使った LDA (Latent Dirichlet Allocation) のわかりやすい解説が書かれている。

LDA は文章が複数の潜在的なトピックから確率的に生成されると仮定したモデルである。

文章に複数のトピックがあると仮定するのは自然だし魅力的ではあるけれども、ぼくはとりあえず検索キーワードのグルーピングがしたい。

検索クエリは短いので、トピックは一種類持っていれば十分だろうと思った。

そこで LDA よりも素朴なモデルを考えた。

モデル

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

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

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

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

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

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

w[i,j] をクエリ i に単語 j が登場したとき 1、さもなくば 0 の値をとる変数とすると、p[j] についての対数尤度は、

 \rm \sum_{i=1}^{L} \sum_{j=1}^{M} w[i,j] \log p[j]+ \\~~(1-w[i,j]) \log(1-p[j])

となる。

さらに、各クエリはインプレッションの回数 IMP[i] だけ登場したとして対数尤度を書くと、

\rm \sum_{i=1}^{L} \sum_{j=1}^{M} {\rm IMP[i]}\{w[i,j] \log p[j]+\\~~ (1-w[i,j]) \log(1-p[j])\}

となる。

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

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

トピックは観測されないので、k について足し合わせて、周辺分布にして対数尤度を書くと、

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

となる。

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

 \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 をデータとして渡すことする。

最後の softmax がベイズの定理でクラスタリングするところです。

//NB.stan
data {
  int<lower=1> K; #num of topics
  int<lower=1> L; #num of queries
  int<lower=1> M; #num of words 
  matrix<lower=0,upper=1>[L,M] W1; # words
  matrix<lower=0,upper=1>[L,M] W2; # words
  vector<lower=0>[L] IMP;
}
parameters {
  simplex[K] phi; #mixing ratio
  vector<lower=0,upper=1>[M] p[K]; #word dist
}
transformed parameters{
  matrix[L,K] lik;
  for(i in 1:L){
    for(j in 1:K){
      lik[i,j] = log(phi[j]) + IMP[i]*(W1[i,]*log(p[j])+W2[i,]*log1m(p[j]));
    }
  }  
}
model {
  for(i in 1:L){
  target += log_sum_exp(lik[i,]);
  }
}
generated quantities{
  matrix[L,K] label;
  for(i in 1:L){
   label[i,] = softmax(lik[i,]')'; 
  }
}

結果

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

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

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

ラベルスイッチングの問題があるので、1 chain しか回してないけれども、トレースプロットの目視と Rhat 指標で一応収束を確認した。

f:id:abrahamcow:20170429024745p:plain

10000回のシミュレーションがだいたい 10 秒くらいで終わった。

混合比率はだいたい 7 対 3 くらいだった。

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

f:id:abrahamcow:20170429024748p:plain

クラスタ 1 と 2 は、明らかに「とは」を含むか含まないかで別れている。

下のほうを見ると、クラスタ 1 は「エクセル」「excel」というワードに出現確率があるのに、クラスタ 2 はほぼない。

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

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

f:id:abrahamcow:20170429024750p:plain

クリック、インプレッション、CTR(クリックスルーレート)いずれもクラスタ 1 が高く、やっぱりエクセルって人気なんだなとわかる。

ただ、ポジションはクラスタ 2 のほうが高く……いやまちがえた、ポジションは数字が小さいほうが上だ。

クラスタ 2 いいとこなしです。

エクセル関連エントリが狙い目です。

R のコード

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

  1. searchConsoleR でデータを読み込む
  2. RMeCab で形態素解析
  3. 下処理
  4. Stan で推定
  5. ggplot で図示、解釈
#あとで使うパッケージをまとめてロードしておく
library(rstan)
library(RMeCab)
library(dplyr)
library(tidyr)
library(cowplot)
library(searchConsoleR)
#######
#サーチコンソールでは 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)
#######
#分析したいランディングページを決める
url1 <- "http://abrahamcow.hatenablog.com/entry/2014/12/19/221230"
LP1 <-dplyr::filter(scdata,page==url1,impressions>=100)
#形態素にバラす
query <-sapply(LP1$query,FUN=function(x)unlist(RMeCabC(x)))
UW <-unique(unlist(query))
M<-length(UW)
BoW1 <- t(sapply(query, function(x){as.integer(UW %in% x)}))
BoW2 <- t(sapply(query, function(x){as.integer(!(UW %in% x))}))
L=nrow(BoW1)
K=2 #トピック数を決める
#
NB <- stan_model("NB.stan")
#######
#Stan の出番
rstan_options(auto_write = TRUE)
dat4stan <-list(M=M,L=L,W1=BoW1,W2=BoW2,IMP=LP1$impressions,K=K)
fit <- sampling(NB,dat4stan,iter=10000, chain=1, seed=49)

traceplot(fit,c("phi","p"))
all(summary(fit)$summary[,"Rhat"]<1.1,na.rm = TRUE)

#結果の取り出し
phat<-matrix(get_posterior_mean(fit,"p"),K,M,byrow = TRUE)
lab <-matrix(get_posterior_mean(fit,"label"),L,K,byrow = TRUE)
get_posterior_mean(fit,"phi")

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

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

df_lower<-as_data_frame(t(matrix(summary(fit,pars="p")$summary[,"2.5%"],K,M,byrow = TRUE))) %>% 
  mutate(id=1:M) %>% 
  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)

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

参考文献

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

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

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

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

関連エントリ

abrahamcow.hatenablog.com