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

廿TT

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

[searchConsoleR]テキストマイニングことはじめ:検索キーワードの視覚化(ワードクラウド、ワードカウント、共起ネットワーク)

R SEO

はじめに

フリーソフトだけでテキストマイニングしたい。

R と MeCab を使います。

R のパッケージもいっぱい使います。

library("googleAuthR")
library("searchConsoleR")
library("geomnet")
library("wordcloud")
library("RColorBrewer")
library("dplyr")
library("RMeCab")
library("tidyr")
library("ggrepel")
library("cowplot")

データ取得

Google サーチコンソールからデータを取ってきます。

gar_auth()
sc_websites <- list_websites() #ウェブサイトのリストを取得
dat_sa <- search_analytics(siteURL = sc_websites[1,1], #分析対象のURLを指定
                          startDate = "2017-01-01",
                          endDate = "2017-01-31", #90日前のデータまでが取得可能
                          dimensions = c("query"),
                          searchType = "web", rowLimit = 5000)

ワードクラウド

下処理はこんなです。

subfunc <- function(x){
tmp <-  unlist(RMeCabC(x))
dplyr::as_data_frame(matrix(tmp[names(tmp)!="助詞"],nrow=1))
}

queries <-sapply(dat_sa$query[dat_sa$clicks>0],subfunc)

queries_df <- bind_rows(queries) %>% 
  mutate(clicks=dat_sa$clicks[dat_sa$clicks>0],
         combn=unname(sapply(queries,length)),
         rownumber=1:n())%>% 
  gather(position,query,-clicks,-combn,-rownumber) %>% 
  dplyr::filter(!is.na(query)) %>% 
  arrange(rownumber)

まず RMeCabC 関数で, 形態素解析して検索クエリをバラバラにしています。

次にクリック数ベースで語句の頻度をカウントします。

助詞はカウントから除きました。

上記の処理をすると、こんなデータフレームが得られます。

> head(queries_df)
# A tibble: 6 × 5
  clicks combn rownumber position    query
   <dbl> <int>     <int>    <chr>    <chr>
1    172     2         1       V1     花札
2    172     2         1       V2   ルール
3     99     3         2       V1     花札
4     99     3         2       V2       花
5     99     3         2       V3   合わせ
6     98     2         3       V1 ポアソン
  • clicks: 一個のクエリに対するクリック数
  • combn: 一個のクエリ内の語の数
  • rownumber: 各クエリに降った通し番号
  • position: 検索後の何番目に対応するか(V1 が「花札」、V2 が「ルール」の場合、元のクエリは「花札 ルール」)
  • query: 検索クエリに現れた語句

これを集計して語句と頻度だけのテーブルを作りましょう。

頻度は重複を認めてカウントすることにします。

例えば、「花札」というクエリが 10 クリック、「花札 ルール」というクエリが 5 クリックあったとしたら、「花札」は頻度は 15、「ルール」の頻度は 5 です。

wordcount_df <- queries_df %>% group_by(query) %>%
  summarise(freq=sum(clicks)) %>% 
  arrange(desc(freq))

集計結果はこんなふうです。

> head(wordcount_df)
# A tibble: 6 × 2
   query  freq
   <chr> <dbl>
1   花札  1162
2486
3 合わせ   478
4   近似   476
5 ルール   455
6391

wordcloud パッケージを使うとかんたんにワードクラウドが作成できます。

#プロットは pdf に書き出すことにする
pdf("wordcloud.pdf",family="Japan1GothicBBB") 
wordcloud(wordcount_df$query,wordcount_df$freq,
          colors = brewer.pal(8,"Dark2"),
          use.r.layout=FALSE)
dev.off()

頻度が大きいほど文字が大きくなります。

今回、色は brewer.pal(8,"Dark2") にしました。

RColorBrewer | Rのカラーパレットの使い方 などを参考に各位好きな色を使ってください。

f:id:abrahamcow:20170219155532p:plain

ちなみに wordcloud2 というパッケージを使うとインタラクティブなワードクラウドが作成できます。が、ここでは解説しません。

ワードクラウドはぱっとみかっこいいけど、頻度を正確に把握するにはふつうの棒グラフのほうがいいでしょう。

ワードカウント

ふつうの棒グラフです。

頻度が 100 以上の語句だけプロットします。

p_wordcount <-ggplot(wordcount_df[wordcount_df$freq>=100,])+
  geom_bar(aes(x=reorder(query,freq),y=freq),stat = "identity")+
  coord_flip()+
  theme_cowplot()+
  xlab("")+ylab("")

ggsave("wordcount.pdf", p_wordcount, family="Japan1GothicBBB")

f:id:abrahamcow:20170219160424p:plain

共起ネットワーク

単語どうしがおなじクエリ内に存在することを「共起」と呼ぶことにします。

共起した語の関係をみることで、キーワードをグルーピングしたり、検索者のニーズをおおまかに把握したりできます。

図示には geomnet パッケージを使います。

集計です。

queries_single <- dplyr::filter(queries_df,combn == 1)
queries_comb <- dplyr::filter(queries_df,combn > 1)

unique_queries <- unique(queries_comb$query)
ulen <- length(unique_queries)
wordforce <- vector("list",length=ulen)
for(i in 1:ulen){
  num = with(queries_comb,rownumber[query == unique_queries[i]])
  to  = with(queries_comb,query[rownumber %in% num])
  wordforce[[i]] <- data_frame(from=unique_queries[i],to=to)
}

wordforce_df <-bind_rows(wordforce) %>% 
  dplyr::filter(from != to) %>% 
  group_by(from,to) %>% 
  ungroup()

wordforce_df_out <- data_frame(from=queries_single$query,to=queries_single$query) %>% 
  bind_rows(wordforce_df) %>% 
  mutate(query=from) %>% 
  left_join(wordcount_df,by="query") %>% 
  arrange(desc(freq)) %>% 
  select(-query)

最終的に得られるデータフレームはこんなふうです。

> head(wordforce_df_out)
# A tibble: 6 × 3
   from     to  freq
  <chr>  <chr> <dbl>
1  花札 ルール  1162
2  花札     花  1162
3  花札 合わせ  1162
4  花札      4  1162
5  花札     人  1162
6  花札      5  1162

from と to は共起関係を示しています。

from が「花札」、to が「ルール」の場合、「花札→ルール」というネットワークのパスが引かれます。

図示の都合上 、from と to を作成していますが、今回共起には矢印の「向き」がないと考えています。

そのため「花札→ルール」というパスが存在すれば「ルール→花札」というパスも必ずあります。

freq は from の語の頻度です。

頻度が 200 を超える単語を含むネットワークのみを図示します。

wordforce_df_top <- wordforce_df_out %>% 
  dplyr::filter(freq>=200 | to %in% from[freq > 200])

p_wordforce<-ggplot(wordforce_df_top,aes(from_id = from, to_id = to)) +
  geom_net(layout.alg = "fruchtermanreingold",
           aes(size=freq),
           directed = FALSE,
           colour="lightpink",
           labelcolour = 'black',
           ecolour = "grey70",
           labelon = TRUE,repel = TRUE,
           labelgeom = "text") +
  theme_net()+
  theme(legend.position="none")+
  xlim(c(-0.1, 1.1))

ggsave("wordforce.pdf", p_wordforce, family="Japan1GothicBBB")

f:id:abrahamcow:20170219162804p:plain

丸の大きさが語句の頻度に対応します。

点と点の距離に意味はありません。見やすくなる位置に適当に並べてくれています。

左にトランプや花札関連のキーワード、右にラップ関連のキーワード、下のほうに統計・データ処理関係のワードが見て取れます。

これで当ブログのユーザー層がだいたい把握できますね。

思い切って全キーワードのネットワークも書いてみます。

p_wordforce_all <- ggplot(wordforce_df_out,aes(from_id = from, to_id = to)) +
  geom_net(layout.alg = "fruchtermanreingold",
           aes(size=freq),
           fontsize=1,
           directed = FALSE,
           colour="lightpink",
           labelcolour = 'black',
           ecolour = "grey70",
           labelon = TRUE,repel = TRUE,
           labelgeom = "text") +
  theme_net()+
  theme(legend.position="none")+
  xlim(c(-0.1, 1.1))

ggsave("wordforce_all.pdf", p_wordforce_all, family="Japan1GothicBBB")

はい。

f:id:abrahamcow:20170219163323p:plain


エンジニアのための データ可視化[実践]入門 ~D3.jsによるWebの可視化 (Software Design plus)

エンジニアのための データ可視化[実践]入門 ~D3.jsによるWebの可視化 (Software Design plus)

未知の変化点があるモデルでは AIC が使えない

R 時系列

モデル

時系列データ x_t (t=1,\ldots, n) があるとします.
このデータが, 変化点( \tau)以前では平均 \mu _1, 標準偏差 1 の正規分布に従い, 変化点から後には平均 \mu_2, 標準偏差 1 の正規分布に従うと考えます.
標準偏差は既知とします.

 \displaystyle x_t =\begin{cases} \mu_1 + \varepsilon _t & t \le \tau \\ \mu_2 + \varepsilon _t & t > \tau \end{cases}

ここで \varepsilon _t は標準正規分布に従う確率変数です.

変化点 \tau最尤推定するには, 対数尤度関数に \tau が与えられたときの \mu _1, \mu _2最尤推定量(標本平均)を代入して尤度が最大になる点を探してやればよさそうです.

最大化すべき対数尤度関数は以下です.

 \displaystyle l(\tau)= \sum_{i=1}^{\tau} (\log(\phi (x_i | \hat \mu_1))+ \sum_{j=\tau+1}^{n} (\log(\phi (x_j | \hat \mu_2)).

ここで \phi標準偏差 1 の正規分布の密度関数,
 \displaystyle \hat \mu_1 = \frac{1}{\tau}\sum_{i =1}^{\tau} x_i,  \displaystyle \hat \mu_2 = \frac{1}{n-\tau}\sum_{i = \tau+1}^{n} x_i
です.

R で推定

乱数で適当なデータを作って, \tau, \mu _1, \mu _2 を推定してみます.

\tau=50, \mu _1=-1, \mu _2=1, n=100と設定しました.

set.seed(1)
x=c(rnorm(50,-1),rnorm(50,1)) #データの生成
ll1_f <- function(tau,n,x){ #尤度関数の定義
  sum(dnorm(x[1:tau],mean(x[1:tau]),log=TRUE))+
    sum(dnorm(x[(tau+1):n],mean(x[(tau+1):n]),log=TRUE))
}
l1v <- sapply(1:99,ll1_f,n=100,x=x) #尤度が最大になる点を総当りで探す
cp <-which.max(l1v)
mu1hat = mean(x[1:cp])
mu2hat = mean(x[(cp+1):100]) 
plot(x,type="b")
segments(x0=0,y0=mu1hat,x1=cp,y1=mu1hat,col="red",lwd=2)
segments(x0=cp+1,y0=mu2hat,x1=100,y1=mu2hat,col="red",lwd=2)

f:id:abrahamcow:20170214010307p:plain

無事, 変化点と \mu _1, \mu _2 が推定されました.

次に AIC で変化点ありのモデルと変化点なしのモデルを比較してみます.

変化点なしのモデルは,

 \displaystyle x_t = \mu_0 + \varepsilon _t

です.

変化点なしのモデルのパラメータは  \mu_0 の 1 個, 変化点ありのモデルのパラメータは \tau, \mu _1, \mu _2 の 3 個です.

ll1 <-max(l1v) #変化点ありのモデルの対数尤度
ll0 <- sum(dnorm(x,mean(x),log=TRUE)) #変化点なしのモデルの対数尤度
-2*ll0+2 #AIC
-2*ll1+2*3

変化点なしのモデルの AIC は367.345, 変化点ありのモデルの AIC は 269.65 となり, 変化点ありのモデルのほうが AIC が小さい=良いモデルとして選択されました。

ここからが本題です。

シミュレーション

あえて変化点なしのモデルで乱数を生成して, 変化点ありのモデルとなしのモデルで AIC を比較してみます.

真のモデルは変化点なしです.

どちらの AIC が小さくなるでしょうか.

simAIC <- function(i,n){
  x=rnorm(n)
  l1v <- sapply(1:(n-1),ll1_f,n=n,x=x)
  cp <-which.max(l1v)
  ll1 <-max(l1v)
  ll0 <- sum(dnorm(x,mean(x),log=TRUE))
 c(AIC0 = -2*ll0 + 2,
   AIC1 = -2*ll1 + 2*3)
}

library(parallel)
n10 <-simplify2array(mclapply(1:10000,simAIC,n=10,mc.cores=detectCores()-1))
n100 <-simplify2array(mclapply(1:10000,simAIC,n=100,mc.cores=detectCores()-1))
n1000 <-simplify2array(mclapply(1:10000,simAIC,n=1000,mc.cores=detectCores()-1))

mean(n10["AIC0",] < n10["AIC1",])
mean(n100["AIC0",] < n100["AIC1",])
mean(n1000["AIC0",] < n1000["AIC1",])

サンプルサイズ n を 10, 100, 1000 と変化させて, それぞれ 1 万回のシミュレーションをおこないました.

変化点なしのモデルが選択された割合は以下のようになりました.

 n 割合
10 76%
100 47%
1000 28%

サンプルサイズが増えても, 正しいモデルが選択されない, というかサンプルサイズが増えるほど, 変化点ありのモデルが選ばれやすくなることがわかります.

なんでこうなるのか, ではどうしたらいいのかはそのうち書きます.