Tokyo.R #76の後日譚々

Tokyo.R #76 @DeNA にてLTを行ってきました。そのスライドがこちらになります。

 

speakerdeck.com

 

発表での伝わりやすさを考慮して、厳密さに関して少し欠落しているところやLDAなのにディリクレ分布の説明が出てこないところ等、至らぬところは目を瞑って、参考文献を読んでいただけると幸いです。

 

コードはこちらの記事にまとめました。

moratoriamuo.hatenablog.com

 

当初のタイトルとしては、『トピックモデル ~江戸時代のレシピを添えて~』というタイトルでした。
人文学オープンデータ共同利用センター

データセット | 人文学オープンデータ共同利用センター
江戸料理レシピデータセットがあり、
これを現代風のレシピにしたものがクックパッドにあるのを発見しました。

cookpad.com


これにトピックモデルを適用してみたが、あまり卵料理が多かったのと、単純に量が少なかった(70件ほど)せいか、あまり面白い結果が出なかったので、急遽路線を現代のレシピに変更して、今回の発表になりました。


それ以前のLT案として、

・『占いは統計学?占い自動生成人工知能 uranAIを作ろう』
データが集められなかったため断念しました。あと、作ってしまうと何かしらまずいことが起きるかもしれないという恐れもありました。
ただ、やってみたい気持ちもあるので、もしも占いのテキストデータを持っている人がいれば、お話ししたいところです。占いではなく、おみくじでも可能です。

・『お酒の失敗談をテキストマイニングする』
少しいじっていたのですが、発表の何かしらの規約に引っかかる可能性を危惧して断念しました。

 

というものがありましたが、これらはボツ案としてお蔵入りです。

 

スライド作成については、Rmarkdownでスライド作ってみたかったのですが、時間の関係上、習得が間に合わず断念し、パワポで作りました。スライドのデザインとしては、クッキング番組風に、「このように作ったものがこちらになります」 みたいなイメージで作りました。色合いもそんな感じで、オレンジ主体になってます。

 

登壇する機会を与えてくださったTokyo.R運営の皆さま、参加者や関係者の皆さま、ありがとうございました。またネタを考えて登壇したいと思っているので、引き続きよろしくお願いいたします。

 

 

トピックモデルで1週間の献立をレコメンドする衝動遊戯:第1ゲーム

 

【ABSTRACT】

料理において、最も面倒くさいことの1つが献立決めである。献立決めのレコメンドエンジンを作れば、この料理への心理的障壁を乗り越えることができる。料理の献立決めのレコメンドエンジンの要件として、同じようなものが続けてレコメンドされないことが求められる。なぜならば、同じものばかりでは飽きてしまうし、栄養の偏りも出てしまうためだ。したがって、この記事では、クックパッドのレシピに、トピックモデルを適用して、レシピをトピックごとにクラスタリングし、クラスタ別で1週間分のレシピを提案するレコメンドエンジンを作成する。

【OBJECTIVE】

料理において、最も面倒くさいことの1つが献立決めである。献立を考える手間を省くため、レコメンドエンジンを作成したい。最も簡単なものとして、レシピ集をランダムに開いて決定する方法が思い浮かぶが、完全にランダムに決定してしまうと、同じものばかり出てしまう可能性が生じる。レシピ集の料理の種類の割合によって出てくるものに偏りが出るのだ。同じものばかりでは飽きてしまうし、栄養のバランスも崩れてしまう。これを避けるために、レシピをクラスタリングして、1週間分の献立を、クラスタが被らないようにしてランダムにレコメンドする方法が考えられる。以降、このようなレコメンドエンジンを作成する。

【METHODS】

f:id:moratoriamuo271:20190305004151p:plain

このような手順でRで作成する。

【RESULTS】

材料確保(スクレイピングする)

使用するライブラリと乱数の種。

library(tidyverse)
library(stringr)
library(rvest)
library(RMeCab)
set.seed(20190302)

クックパッドのカテゴリのリストが
https://cookpad.com/category/list にある。「今日のご飯・おかず」、「お菓子」、「パン」、「離乳食」、「その他」のカテゴリから、「今日のご飯・おかず」を選択する。13214ページあって、レシピ141355品が掲載されている(この件数は日々変動しているので注意)。

14万件全てを取ってくるのは大変なので、50ページ分だけランダムに取ってくることにする。

url <- "https://cookpad.com/category/177?page="
allrepnum <- 1:13214
repset <- sample(allrepnum, 50, replace = FALSE)

ページからレシピの番号(https://cookpad.com/recipe/の後に続く数字)を入手する関数を作成

fget_recipenumlist <- function(url,repset){
  rec <- vector("list", length(repset))
  for(i in 1:length(repset)){
    urlfull <- str_c(url, repset[i] %>% as.character())
    rec[[i]] <- read_html(urlfull) %>% html_nodes(css = "a") %>% 
      html_attr("href") %>% str_subset("/recipe/(\\d\\d\\d\\d)") 
  }
  return(rec)
}

関数を適用する。重複して取られてしまうケースがあるので重複削除。

reclist <- fget_recipenumlist(url, repset)
reclist <- reclist %>% unlist() %>% unique()    #重複を消す

入手したレシピの番号からレシピの中身を入手する関数を作成。ごり押し気味。

fget_cookpads <- function(reclist){
  rec <- vector("list", length(reclist))
  for(i in 1:length(reclist)){
    urlfull <- str_c("https://cookpad.com",reclist[i])
    temp_text  <- read_html(urlfull) %>% html_nodes(css = "p") %>%html_text()
    loc_start <- temp_text %>% str_which("\nメンバー検索") + 1
    loc_end <- temp_text %>% str_which("(\\d\\d)/(\\d\\d)/(\\d\\d)") %>% as.list()
    if(length(loc_end)==0){
      loc_end <- temp_text %>% str_which("お困りの方はこちら") %>% as.list()
    }
    loc_end <- loc_end[[1]] - 1
    rec[[i]] <- temp_text[loc_start:loc_end]
    if(i%%100==0)print(i)    #経過状況を見るため100個おきにプリント
    Sys.sleep(1)
  }
  return(rec)
}

関数を適用する。改行を消すのと、同じレシピは一つの文書にまとめる。

cookpads <- fget_cookpads(reclist)
#改行の\nを消して各レシピを一つの文書にまとめる
cookpads <- purrr::map(cookpads, str_remove_all,"\n") %>% purrr::map(str_c, collapse="")

入手したレシピを各レシピごとで保存する。加えて、レシピを一つの文書にまとめて保存する。

#レシピの番号
reclistnum <- str_extract(reclist, "(\\d){4,}")

#レシピの保存(各レシピ)
for(i in 1:length(cookpads)){
  filename <- str_c("docs_nowrecipe/",reclistnum[i] %>% as.character()) %>% str_c(".txt")
  write(cookpads[[i]],filename)
}

#全レシピ連結で保存
allcookpads <- cookpads %>% unlist() %>% str_c(collapse = "")
write(allcookpads, "allnowdocs.txt")

これにより、500件のレシピが保存できた。

材料確認(ワードクラウドで遊んでみる)

library(tidyverse)
library(stringr)
library(RMeCab)
library(tm)
library(topicmodels)

どんな単語があるかを確認。自立語(動詞や形容詞)で見るケースと名詞で見るケース。

all_freq <- RMeCabFreq("allnowdocs.txt")
## file = allnowdocs.txt 
## length = 4110
all_freq %>% filter(Info2=="自立") %>% arrange(-Freq) %>% head(20)
##          Term  Info1 Info2 Freq
## 1        する   動詞  自立 1366
## 2      入れる   動詞  自立  897
## 3        切る   動詞  自立  412
## 4      加える   動詞  自立  332
## 5      炒める   動詞  自立  325
## 6      混ぜる   動詞  自立  285
## 7        なる   動詞  自立  193
## 8      かける   動詞  自立  174
## 9  出来上がる   動詞  自立  171
## 10       焼く   動詞  自立  167
## 11     茹でる   動詞  自立  153
## 12     食べる   動詞  自立  138
## 13       熱す   動詞  自立  117
## 14       煮る   動詞  自立  116
## 15       軽い 形容詞  自立  113
## 16     大きい 形容詞  自立  107
## 17       作る   動詞  自立   98
## 18       取る   動詞  自立   87
## 19       盛る   動詞  自立   81
## 20     つける   動詞  自立   80
all_freq %>% filter(Info1=="名詞"&Info2!="数"&Info2!="接尾") %>% arrange(-Freq) %>% head(20)
##          Term Info1    Info2 Freq
## 1  フライパン  名詞     一般  307
## 2          塩  名詞     一般  272
## 3          火  名詞     一般  254
## 4          水  名詞     一般  241
## 5          油  名詞     一般  182
## 6          ♪  名詞 サ変接続  178
## 7          鍋  名詞     一般  163
## 8          肉  名詞     一般  160
## 9          卵  名詞     一般  136
## 10          +  名詞 サ変接続  131
## 11          U  名詞     一般  130
## 12     玉ねぎ  名詞     一般  129
## 13       よう  名詞   非自立  129
## 14         皮  名詞     一般  121
## 15     レシピ  名詞     一般  117
## 16         味  名詞     一般  117
## 17       調味  名詞 サ変接続  114
## 18       材料  名詞     一般  113
## 19         皿  名詞     一般  112
## 20       好み  名詞     一般  106

名詞と自立語でワードクラウドを作ってみる。(警告が現れるが、楽しく可視化したいだけなので無視している。)

set.seed(271)
corpus <- docDF("alldocs.txt", type = 1)
## file_name =  ./alldocs.txt opened
## number of extracted terms = 1166
## now making a data frame. wait a while!
corpus <- corpus[81:nrow(corpus),]    #記号とか余計なものを消した
options(warn=-1)    #ワードクラウドの警告が邪魔なので消した
corpus1 <- corpus %>% filter(POS1 == "名詞"&POS2=="一般"| POS2 == "自立") 
wordcloud::wordcloud(corpus1$TERM, corpus1$alldocs.txt, min.freq= 5, scale=c(15,0.5),
                     family ="JP1", random.color=TRUE, color = rainbow(5))

options(warn=0)

名詞でワードクラウドを作ってみる。

options(warn=-1)
corpus2 <- corpus %>% filter(POS1 == "名詞"&POS2=="一般")
wordcloud::wordcloud(corpus2$TERM, corpus2$alldocs.txt, min.freq= 5, scale=c(15,0.5),
                     family ="JP1", random.color=TRUE, color = rainbow(5))

options(warn=0)

調理(トピックモデル)

トピックモデルの1つのLDA(Latent Dirichlet Allocation)を用いてクラスタリングを行う。

名詞と動詞に絞って、文書ターム行列を作成。記号等も削除している。削除の仕方はごり押し気味。

dir_recipe <- "docs_nowrecipe"
DTM <- docDF(dir_recipe, type = 1)
DTM <- DTM[512:4106,]

LDAを適用するために文書ターム行列の形に整形する。

DTM_nouns_verbs <- DTM %>% filter( (POS1 == "名詞"&POS2 == "一般") | POS2 == "動詞" )%>% 
  select(-c(POS1,POS2)) %>% 
  data.frame(row.names = 1) %>% t() %>% 
  as.DocumentTermMatrix(weighting = weightTf)

後でperplexity計算するため8:2でtrainとtestに分割する。

#500*1720の行列
DTM_nouns_verbs_train <- DTM_nouns_verbs[1:400,]
DTM_nouns_verbs_test <- DTM_nouns_verbs[401:500,]

モデル選択(トピック数の決定)の方法として

  • パープレキシティ perplexityによる評価 負の対数尤度から計算される値。testデータを使用。低い方が良いモデル。

  • {ldatuning}パッケージの使用 4つの論文で提案された指標で評価する。

  • 変分下限でモデル評価 変分ベイズを使っているならば評価に使える。

  • ディリクレ過程でトピック数もモデルに組み込む 階層ディリクレ過程を用いるとトピック数の推定が可能。勉強不足のため、これについてはまだよく分かっていない。pythonのgenismというものを使うとできるかもしれない。以下のブログが参考になるかもしれない。 https://qiita.com/yamano357/items/90863cf15326ebd19231

以下では、{ldatuning}パッケージで、だいたいのあたりをつけて、その中からperplexityで評価して、トピック数の選択を行う。
(トピック数のベストな設定方法は結局のところ分かっていないので、もう少し調べていきたいところ。混合ユニグラムモデルまでは過去にシミュレーションデータが作れているので、LDAのシミュレーションデータを作って各指標をチェックする記事をそのうち書きたいと思っている。)

ここでは、{ldatuning}のパッケージでトピック数を2から97のトピック数を5個刻みにして評価する。

result <- ldatuning::FindTopicsNumber(
  DTM_nouns_verbs_train,
  topics = seq(from = 2, to = 97, by = 5),
  metrics = c("Griffiths2004", "CaoJuan2009", "Arun2010", "Deveaud2014"),
  method = "Gibbs",
  control = list(seed = 777),
  mc.cores = 4L,
  verbose = TRUE
)
## fit models... done.
## calculate metrics:
##   Griffiths2004... done.
##   CaoJuan2009... done.
##   Arun2010... done.
##   Deveaud2014... done.
ldatuning::FindTopicsNumber_plot(values = result)

これにより、12~27あたりが良いだろうと見当がつくので、12,17,22,27のトピック数でperplexityを評価する。

result_tp12 <- LDA(DTM_nouns_verbs_train,12,method = "Gibbs")
result_tp17 <- LDA(DTM_nouns_verbs_train,17,method = "Gibbs")
result_tp22 <- LDA(DTM_nouns_verbs_train,22,method = "Gibbs")
result_tp27 <- LDA(DTM_nouns_verbs_train,27,method = "Gibbs")
perplexity(result_tp12, DTM_nouns_verbs_test)
## [1] 488.4454
perplexity(result_tp17, DTM_nouns_verbs_test)
## [1] 474.552
perplexity(result_tp22, DTM_nouns_verbs_test)
## [1] 471.9647
perplexity(result_tp27, DTM_nouns_verbs_test)
## [1] 473.9953

以上のldatuningとperpleixtyからトピック数22に決定。

各トピックの頻出単語を確認

terms(result_tp22,5)
##      Topic 1      Topic 2 Topic 3    Topic 4    Topic 5  Topic 6   
## [1,] "水"         "皮"    "チーズ"   "塩"       "レンジ" "火"      
## [2,] "量"         "大根"  "牛乳"     "コショウ" "耐熱"   "ベーコン"
## [3,] "ご飯"       "写真"  "ソース"   "中火"     "皿"     "好み"    
## [4,] "アク"       "作者"  "バター"   "上"       "ラップ" "とろみ"  
## [5,] "マヨネーズ" "味"    "コンソメ" "皿"       "容器"   "小麦粉"  
##      Topic 7 Topic 8 Topic 9 Topic 10 Topic 11 Topic 12     Topic 13  
## [1,] "部分"  "肉"    "ネギ"  "玉ねぎ" "材料"   "鍋"         "油"      
## [2,] "砂"    "鶏"    "生姜"  "パスタ" "器"     "火"         "こしょう"
## [3,] "肝"    "むね"  "煮汁"  "トマト" "鶏肉"   "蓋"         "火"      
## [4,] "醤油"  "塩"    "鯛"    "オイル" "油"     "弱火"       "豚"      
## [5,] "黒"    "衣"    "酒"    "ソース" "一口"   "じゃがいも" "強火"    
##      Topic 14     Topic 15   Topic 16     Topic 17   Topic 18  
## [1,] "フライパン" "塩"       "フライパン" "豆腐"     "キャベツ"
## [2,] "卵"         "きゅうり" "ボウル"     "水気"     "豚肉"    
## [3,] "大さじ"     "ボール"   "サラダ油"   "皿"       "ごま油"  
## [4,] "油"         "胡椒"     "両面"       "水分"     "好み"    
## [5,] "白"         "玉ねぎ"   "玉ねぎ"     "キッチン" "火"      
##      Topic 19       Topic 20 Topic 21 Topic 22
## [1,] "レシピ"       "冷蔵庫" "野菜"   "白菜"  
## [2,] "ドレッシング" "酢"     "水"     "ご飯"  
## [3,] "ねぎ"         "大葉"   "味"     "だし"  
## [4,] "中華"         "好み"   "最後"   "幅"    
## [5,] "サラダ"       "水気"   "鍋"     "葉"

解釈しやすそうなものを4つ取り上げた。

terms(result_tp22,5)[,c(3,8,15,18)] %>% as.data.frame()
##    Topic 3 Topic 8 Topic 15 Topic 18
## 1   チーズ      肉       塩 キャベツ
## 2     牛乳      鶏 きゅうり     豚肉
## 3   ソース    むね   ボール   ごま油
## 4   バター      塩     胡椒     好み
## 5 コンソメ      衣   玉ねぎ       火
  • Topic 3 グラタン, シチュー, スープ
  • Topic 8 唐揚げなど鶏肉料理
  • Topic 15 サラダ
  • Topic 18 野菜炒め

といったところだろうか。

文書ごとのトピック分類を行う。推定されたトピック分布の最大の値をとるトピックをその文書のトピックとして分類することにする。

result_tp22_cdoc <- result_tp22 %>% tidytext::tidy(matrix = "gamma") %>% as.data.frame()
result_tp22_cdoc$document <- result_tp22_cdoc$document %>% str_remove_all("X")
doctop <- result_tp22_cdoc %>% group_by(document) %>% filter(gamma==max(gamma)) %>% ungroup() %>% select(-gamma)
head(doctop)
## # A tibble: 6 x 2
##   document    topic
##   <chr>       <int>
## 1 109612.txt      1
## 2 1289543.txt     1
## 3 1307357.txt     1
## 4 139698.txt      1
## 5 1426704.txt     1
## 6 1457975.txt     1

400になると思ったら、もっとあった。重複した文書があるのは、おそらくgammaの最大値が一緒のトピックがあったのだろう。ここでは、重複も含めてしまうことにする。

これを保存する。

write_csv(doctop, "doctop.csv")
doctop <- read_csv("doctop.csv")

ここまでやったものを保存しておけば、いちいち上の面倒くさい計算を回さず、レコメンドしたいときに読み込んで使える。

完成(レコメンド)

#一週間分の料理(7食分)をレコメンド
frecommend7url <- function(doctop){
  seventopics <- sample(x = unique(doctop$topic), size = 7, replace = FALSE)
  rec7 <- character(7)
  for(i in 1:7){
    doctopfil <- doctop %>% filter(topic==seventopics[i])
    rec7[i] <- sample(x = doctopfil$document, size = 1)
  }
  rec7 <- rec7 %>% str_remove(".txt") %>% str_c("https://cookpad.com/recipe/", .)
  return(rec7)
}

frecommend7url(doctop)
## [1] "https://cookpad.com/recipe/3169287"
## [2] "https://cookpad.com/recipe/1370545"
## [3] "https://cookpad.com/recipe/2733198"
## [4] "https://cookpad.com/recipe/2866455"
## [5] "https://cookpad.com/recipe/2721802"
## [6] "https://cookpad.com/recipe/2430458"
## [7] "https://cookpad.com/recipe/1120158"

トピックの異なるレシピのurlが7つ提案されるレコメンドエンジンが完成された。提案されたレシピを実際に見てみると、トピックが違うレシピが提案されていることが確認される。

【CONCLUSIONS】

LDAを用いてクックパッドのレシピをトピックで分類して、トピックが異なるレシピを7つ提案するレコメンドエンジンが作成された。提案されたレシピを実際に見てみると、トピックが違うレシピが提案されていることが確認された。
このレコメンドエンジンは、とりあえず、トピック別にレシピが提案されるが、レコメンドエンジンとしての性能があまり良いものとは言えないだろう。要件として、冷蔵庫の中身とか手持ちの材料の賞味期限とか季節の食べ物とか必要な分量とかの情報を組み込んだレコメンドが望ましい。そんな性能のよいレコメンドエンジンが作れているなら、クックパッドにとっくに売り込んでいるという話である。

【REFERENCES】

トピックモデル(混合ユニグラムモデル・EMアルゴリズム)の学習記録

 

【学習動機】

前回、ユニグラムモデルを学習した。

moratoriamuo.hatenablog.com

これからは、混合ユニグラムモデルを学習する。混合ユニグラムモデルをEMアルゴリズムで推定する。

【学習内容】

混合ユニグラムモデル(mixture of unigram models)は各文書がそれぞれ1つのトピックをもっていてトピックごとに異なった単語分布をもち、そこから単語が生成されるモデル。文書のトピックはトピック分布から生成され、割り当てられる。Rで、このモデルの人工データを生成し、EMアルゴリズムでパラメータを最尤推定する。

人工データ生成

ライブラリの読み込み。

library(tidyverse)
library(ggplot2)
library(MCMCpack)    #rdirichletでディリクレ乱数を生成する
## Warning: package 'MCMCpack' was built under R version 3.5.2
set.seed(271)

トピックとトピックごとに頻出する語彙を設定する。トピックとして、しりとりとRとアジカンを選択した。しりとりを選んでしまったが、分布から独立に単語が生成されることを考えると例としてふさわしくないことに後で気付いたがそのままにする。語彙の種類が多すぎると大変なので、この程度に絞った。

topics <- c("shiritori","r","akg")
allword <- c("りんご","ゴリラ","ラッパ",
             "宇宙","神","パイプ",
             "スタンダード","センスレス","リライト","ループ" )
nvoc <- length(allword)    #語彙の総数

各トピックに対応する単語分布の設定。トピックを想定しながら手動で設定した。ここは、パラメータを上手く決めてディリクレ分布から乱数発生させてもよい。各トピックの単語分布を決めたら、リストにしてまとめておく。

#トピック: しりとり
df_word_dist_shiritori <- data.frame(
  words = allword,
  prob = c(0.25,0.2,0.2,
           rep(0.05,7))
)

#トピック: R
df_word_dist_r <- data.frame(
  words = allword,
  prob = c(0.01,0.01,0.02,
           0.3,0.2,0.31,
           0.01,0.01,0.015,0.115)
)

#トピック: AKG
df_word_dist_akg <- data.frame(
  words = allword,
  prob = c(rep(0.01,3),
           0.06,0.1,0.01,
           0.2,0.2,0.2,0.2)
)

word_dists <- list(df_word_dist_shiritori, df_word_dist_r, df_word_dist_akg)
names(word_dists) <- c("shiritori","r","akg")

文書の数と各文書の単語数(決めるのが面倒なので一様乱数で生成)を決定。

numD <- 100
vnumW <- sample(140:271,numD,replace = TRUE)
numW <- sum(vnumW)    #総単語数

トピック分布をディリクレ分布によって生成し、トピック分布から文書へのトピックの割り当てを生成。

theta <- rdirichlet(1,c(1,2,7))    #alpha=c(1,2,7)
vtopic <- sample(topics, numD, replace = TRUE, prob = theta)

あとで推定されたパラメータと比較するために真のパラメータ(thetaとphi)を保存しておく。

truepara <- list(
  truetheta = theta %>% as.vector(),
  truephi = word_dists
)

以上で指定してものを使えば混合ユニグラムモデルからのデータを生成する関数が作れる。引数は文書の数と各文書の単語数のベクトル、単語分布のデータフレームのリストをとる。

sim_mixunigram <- function(numD, vnumW, word_dists){
  list_word <- vector("list", numD)
  for(d in 1:numD){
    list_word[[d]] <- sample(word_dists[[vtopic[d]]]$words, vnumW[d], replace=TRUE,
                             prob = word_dists[[vtopic[d]]]$prob)
  }
  return(list_word)
}

人口データを発生させてみる。前回同様、頻度テーブルにして文書1の頻度テーブルを見てみる。

DW <- sim_mixunigram(numD=numD, vnumW = vnumW, word_dists = word_dists)    #DocumentWordsList
DWtab <- DW %>% purrr::map(factor,levels = allword) %>% purrr::map(table)
DWtab[[1]]
## 
##       りんご       ゴリラ       ラッパ         宇宙           神 
##            0            2            2           53           29 
##       パイプ スタンダード   センスレス     リライト       ループ 
##           51            1            1            1           33

1つ目の文書のトピックはRだったに違いない。

後の計算をしやすいように、文書ターム行列(Document Term Matrix)に変換しておく。

DTM <- matrix(DWtab %>% unlist() %>% as.numeric(), ncol = nvoc, byrow=T)    #DocumentTermMatrix
colnames(DTM) <- allword

EMアルゴリズム最尤推定

トピックの数は決めておかないといけないので、トピックは3つと天下り的に仮定しておく。

numT <- 3    ##ofTopics

対数尤度を計算する関数を定義しておく。引数には、文書ターム行列とトピック分布のパラメータthetaと各トピックの単語分布のデータフレームのリストphiをとる。

loglikelihood_mug <- function(DTM, theta, phi){
  ##mixture of unigramの対数尤度を計算する
  oneL <- numeric(nrow(DTM))    #文書ごとの尤度ベクトル
  pr <- numeric(nrow(DTM))    #計算用のproductで計算したもの
  
  for(k in 1:length(theta)){
    for(d in 1:nrow(DTM)){
      pr[d] <- prod(purrr::map2_dbl(phi[[k]]$prob, DTM[d,], `^`))    #p(wd|phik)
    }
    oneL <- oneL + theta[k]*pr
  }
  LL <- oneL %>% log() %>% sum()
  return(LL)
}

岩田『トピックモデル』のEMアルゴリズム擬似コードを見ながら、EMアルゴリズム関数を作成。引数には、文書ターム行列とトピック数、全ての語彙、繰り返し回数(対数尤度を見ている限り、想像以上に局所最適への収束が速いので、デフォで10としている)をとる。返り値として、推定されたthetaとphi、対数尤度の履歴のリストを返す。

fEM <- function(DTM, numT, allword, numRep=10){

  #対数尤度履歴
  LLhistory <- numeric(numRep)
  
  #基準化関数
  fScale <- function(x){
    y <- x/sum(x)
    return(y)
  }
  
  #初期値: 全て同じ値で入れてしまうと上手く動かなかった
  theta <- runif(3) %>% fScale()
  phi <- vector("list",numT)
  for(i in 1:numT){
    phi[[i]] <- data.frame(
      words = allword,
      prob = runif(10) %>% fScale()
    )
  }
  
  for(nrep in 1:numRep){
    #次ステップのパラメータを0に初期化
    theta_new <- rep(0,numT); phi_new <- vector("list",numT);
    for(k in 1:numT){
      phi_new[[k]] <- data.frame(
        words = allword,
        prob = rep(0,length(allword))
      )
    }

    for(d in 1:nrow(DTM)){
      Qd <- numeric(numT)
      for(k in 1:numT){
        Qd[k] <- prod(purrr::map2_dbl(phi[[k]]$prob, DTM[d,], `^`))*theta[k]
      }

      Qd <- fScale(Qd)
      for(k in 1:numT){
        theta_new[k] <- theta_new[k] + Qd[k]
        phi_new[[k]]$prob <- phi_new[[k]]$prob + Qd[k]*as.vector(DTM[d,])
      }
    }
    
    theta <- fScale(theta_new)
    for(k in 1:numT){
      phi[[k]]$prob <- phi_new[[k]]$prob %>% fScale()
    }
    
    LLhistory[nrep] <- loglikelihood_mug(DTM, theta, phi)
  }
  
  
  ret <- list(
    EMtheta = theta,
    EMphi = phi,
    LLhistory = LLhistory
  )
}

実際に推定してみて、対数尤度の動きを見てみる。

result <- fEM(DTM = DTM, numT = 3, allword = allword)
result$LLhistory
##  [1] -37824.38 -37803.02 -37803.02 -37803.02 -37803.02 -37803.02 -37803.02
##  [8] -37803.02 -37803.02 -37803.02

二回目には収束している感じだろうか。

推定されたphiをみていく。

result$EMphi[[1]]
##           words       prob
## 1        りんご 0.01059916
## 2        ゴリラ 0.01078621
## 3        ラッパ 0.00885342
## 4          宇宙 0.06004115
## 5            神 0.09626535
## 6        パイプ 0.01028742
## 7  スタンダード 0.19963838
## 8    センスレス 0.20163352
## 9      リライト 0.19970073
## 10       ループ 0.20219465

どうやらアジカンのトピックの単語分布を推定しているようなので、真の分布と比較してみると

truepara[["truephi"]][["akg"]]
##           words prob
## 1        りんご 0.01
## 2        ゴリラ 0.01
## 3        ラッパ 0.01
## 4          宇宙 0.06
## 5            神 0.10
## 6        パイプ 0.01
## 7  スタンダード 0.20
## 8    センスレス 0.20
## 9      リライト 0.20
## 10       ループ 0.20

次は

result$EMphi[[2]]
##           words        prob
## 1        りんご 0.011339662
## 2        ゴリラ 0.008175105
## 3        ラッパ 0.021360759
## 4          宇宙 0.294831224
## 5            神 0.209651899
## 6        パイプ 0.300105485
## 7  スタンダード 0.010812236
## 8    センスレス 0.010284810
## 9      リライト 0.017405063
## 10       ループ 0.116033755

Rのトピックの単語分布を推定しているようなので、真の分布と比較してみると

truepara[["truephi"]][["r"]]
##           words  prob
## 1        りんご 0.010
## 2        ゴリラ 0.010
## 3        ラッパ 0.020
## 4          宇宙 0.300
## 5            神 0.200
## 6        パイプ 0.310
## 7  スタンダード 0.010
## 8    センスレス 0.010
## 9      リライト 0.015
## 10       ループ 0.115

最後に

result$EMphi[[3]]
##           words       prob
## 1        りんご 0.24165554
## 2        ゴリラ 0.19225634
## 3        ラッパ 0.19759680
## 4          宇宙 0.05740988
## 5            神 0.04405874
## 6        パイプ 0.04806409
## 7  スタンダード 0.04806409
## 8    センスレス 0.06809079
## 9      リライト 0.04138852
## 10       ループ 0.06141522

しりとりの真の分布と比較する。

truepara[["truephi"]][["shiritori"]]
##           words prob
## 1        りんご 0.25
## 2        ゴリラ 0.20
## 3        ラッパ 0.20
## 4          宇宙 0.05
## 5            神 0.05
## 6        パイプ 0.05
## 7  スタンダード 0.05
## 8    センスレス 0.05
## 9      リライト 0.05
## 10       ループ 0.05

推定されたtheta(ここではakg,r,shiritoriの順)と真のtheta(shiritori,r,akgの順)を比較してみる。

result$EMtheta %>% round(digits = 3)
## [1] 0.78 0.18 0.04
truepara[["truetheta"]]%>% round(digits = 3)
## [1] 0.031 0.217 0.752

そこそこ近くまでは来ているみたい。

もう少し精度を上げようと単語数や文書数を増やしてみようとしたのだが、fEM関数で小数の積をとりまくることになることがおそらくの原因で、NaNが発生してしまった。こういうときはどうやって対処するのか分からなかったので、数値計算・数値解析のお話しが気になった。

【学習予定】

次は混合ユニグラムモデルを変分ベイズで推定していく。

トピックモデル(ユニグラムモデル)の学習記録

 

【学習動機】

以前、テキストマイニングをRで使って遊んでみたとき, トピックモデルが出てきたが、理論を全く知らないのも悲しいので、岩田『トピックモデル』

www.amazon.co.jp

で学ぶことにした。まずはユニグラムモデルについて学ぶ。

【学習内容】

ユニグラムモデル(unigram model)はBOW(bag of words)表現された文書集合の最も簡単な生成モデル。このモデルの理論を学んだので、Rで人工データを作成し、最尤推定とMAP推定でパラメータの推定を行う。

ライブラリの読み込みと表示桁数の指定。

library(tidyverse)
library(ggplot2)
options(digits=3)
set.seed(536329)

単語分布の設定をする。語彙の種類とそれに付与される確率を決定。

df_word_dist <- data.frame(
  words = c("あじ","いくら","うに","えんがわ","かれい","こはだ","さわら","さんま","しゃこ","すずき"),
  prob = c(0.2,0.1,0.075,0.075,0.05,0.075,0.1,0.2,0.075,0.05)
)
numV <- nrow(df_word_dist)    #単語分布のボキャブラリーの種類(number of Vocabulary)
allword <- df_word_dist$words %>% as.vector()   #全単語(あとでファクタ型にしたときの並べ替えで使う)

文書数と各文書の単語数を指定する。

numD <- 12    #文書の数(number of Documents)
vnumW <- c(20,70,10,80,20,80,10,80,20,80,40,60)    #各文書に含まれる単語数(number of words)のベクトル
numW <- sum(vnumW)    #単語の総数

シミュレーションデータ作成関数が以下のようになる。引数には文書の数と各文書の単語数ベクトル、単語分布のデータフレームをとる。

sim_unigram <- function(numD, vnumW, df_word_dist){
  list_word <- vector("list", numD)
  for(d in 1:numD){
    list_word[[d]] <- sample(df_word_dist$words, vnumW[d], replace=TRUE, prob = df_word_dist$prob)
  }
  return(list_word)
}

これにより人工データを作成し、生成された文書1の単語を見てみる。

DW <- sim_unigram(numD=numD, vnumW = vnumW, df_word_dist = df_word_dist)
DW[[1]]
##  [1] あじ     こはだ   さんま   すずき   えんがわ さんま   うに    
##  [8] さわら   さんま   かれい   すずき   しゃこ   さんま   あじ    
## [15] さわら   あじ     さんま   すずき   あじ     さんま  
## 10 Levels: あじ いくら うに えんがわ かれい こはだ さわら ... すずき

各語彙の頻度テーブルを作成し、生成された文書1のテーブルを見てみる。

DWtab <- DW %>% purrr::map(factor,levels = allword) %>% purrr::map(table)
DWtab[[1]]
## 
##     あじ   いくら     うに えんがわ   かれい   こはだ   さわら   さんま 
##        4        0        1        1        1        1        2        6 
##   しゃこ   すずき 
##        1        3

ユニグラムモデルでは文書によらず単語分布が一緒なので、語彙の頻度テーブルを統一してしまって問題ない。

allWtab <- DW %>% unlist() %>% factor(levels = allword) %>% table()

最尤推定(phiv = Nv/N)を行う。

MLphi <- allWtab / numW
MLphi
## .
##     あじ   いくら     うに えんがわ   かれい   こはだ   さわら   さんま 
##   0.2035   0.1035   0.0842   0.0930   0.0526   0.0579   0.0860   0.1930 
##   しゃこ   すずき 
##   0.0807   0.0456

事前分布にディリクレ分布を設定し、パラメータを(beta,beta,beta,…,beta)として、MAP推定(phiv = (Nv+beta-1)/(N+(beta-1)V) beta:paramofDirichlet)を行う。

beta <- 2    #観測がない単語の推定値を負にしないためにはbetaが1以上とする必要がある
MAPphi <- (allWtab + beta - 1) / (numW + (beta - 1)*numV)
MAPphi
## .
##     あじ   いくら     うに えんがわ   かれい   こはだ   さわら   さんま 
##   0.2017   0.1034   0.0845   0.0931   0.0534   0.0586   0.0862   0.1914 
##   しゃこ   すずき 
##   0.0810   0.0466

ちなみに、ベイズ推定では事前分布にディリクレ分布として、パラメータを(beta,beta,beta,…,beta)とすると、事後分布はDirichlet(phi|beta+N1,beta+N2,…,beta+Nv)が得られるが、ここでは求めない。

MAP推定と最尤推定(ML)と真のパラメータを棒グラフにして並べる。赤がMAP推定値、緑が最尤推定値、青が真の値を示している。

df_est <- data.frame(
    voc = df_word_dist$words,
    phi = df_word_dist$prob,
    MLphi = MLphi %>% as.vector(),
    MAPphi = MAPphi %>% as.vector()
)
df_est_wide <- gather(df_est, key = "estimator",value,-voc)
df_est_wide$voc <- df_est_wide$voc %>% factor(levels = df_word_dist$words)

ggplot(df_est_wide, aes(x=estimator, y=value, fill=estimator)) +
  geom_bar(stat = "identity", position = "dodge") +
  facet_wrap(~voc)

サンプルサイズを10%にすると当然ながら推定精度が落ちる。

set.seed(536329)
numD <- 12    #文書の数(number of Documents)
vnumW <- c(20,70,10,80,20,80,10,80,20,80,40,60)*0.1   #各文書に含まれる単語数(number of words)のベクトル

【学習予定】

次は混合ユニグラムモデルを学ぶ。

moratoriamuo.hatenablog.com

2018年度のRの学習記録

【学習動機】

R初級者からR中級者くらいになりたいと思って、Rをこの1年間メインに使っていた。Tokyo.RやJapan.Rに参加してみたりもした。色々と読んでみたので、それらの本について、簡単に記録しておく。

 

【学習記録】

<おすすめ>

・『RユーザのためのRStudio[実践]入門−tidyverseによるモダンな分析フローの世界−』

いわゆる宇宙本。tidyverseを何も知らない人はまずこれを読むべき。

 

・『Rプログラミング本格入門: 達人データサイエンティストへの道 』

いわゆるパンダ本。一回見たことあるけど、よく分かってなかったとことか気になってたとことか、説明されてて、頭の中整理できて良かった。

 

・『Rではじめるデータサイエンス』

宇宙本の次に読んでみると良いかもしれない。map()やnest()など、分析の効率性を高めるものが知れる。

 

・『Rによるハイパフォーマンスコンピューティング』

並列化の仕方が学べる。

 

・『Rによるスクレイピング入門』

スクレイピングの仕方が学べる。

 

・『欠測データ処理: Rによる単一代入法と多重代入法 (統計学One Point)』

欠測データの代入の仕方が学べる。

 

・『時系列分析と状態空間モデルの基礎: RとStanで学ぶ理論と実装』 

いわゆる隼本。Rの時系列周りの本では、これがかなり良かった。

 

・『R言語徹底解説』

辞書的に使ってる。 {rcpp}パッケージについてのところを読んだりした。

 

<その他読んだもの>

 ・『Rによるやさしいテキストマイニング: 機械学習編』

テキストマイニングで遊べた。

 ・『再現可能性のすゝめ (Wonderful R 3)』

再現可能性ソウルが手に入る。

・『Rグラフィックスクックブック ―ggplot2によるグラフ作成のレシピ集』

ggplot2のグラフの作り方が学べる。練習で使った。読んだのが古い版だからか分からないが、記載されているコードのggplot2の仕様などが古かったりするので注意が必要。

・『一般化線形モデル (Rで学ぶデータサイエンス 10)』

一般化線形モデルのRでのやり方が分かる。

・『Rで学ぶベイズ統計学入門』

{LearnBayes}パッケージにデータセットが豊富。和訳の問題かもしれないが、密度や分布、周辺尤度あたりの言葉の使い方が気になった。

 

<読んでるもの・読みたいもの>

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

・『自然科学研究のためのR入門―再現可能なレポート執筆実践―(Wonderful R 4)』

・『Rパッケージ開発入門 ―テスト、文書化、コード共有の手法を学ぶ』

・『前処理大全[データ分析のためのSQL/R/Python実践テクニック]』

・『Seamless R and C++ Integration with Rcpp (Use R!)』


【学習予定】

上の読んでるもの・読みたいものを消化したいのと、pythonも扱えるようにしていきたい。あと、アルゴリズムとデータ構造周りやコードの書き方も学んでいきたいところ。

 

 

 

何かしらデータベースやデータセットを使いたいときの覚書

【用途】

何かしらデータ分析してみたいけど、よくあるデータセット(irisとかmtcarsとかflightsとかNMISTとか)は見飽きたし、自分がビビッと来て興味が湧いてくるデータをいじってみたいなんてときのためのメモ。

【内容】

・人文学オープンデータ共同利用センター データセット 日本古典籍データセット, 江戸料理レシピデータセット, 日本古典籍くずし字データセット, KMNISTデータセット, 武鑑全集, 顔貌コレクション, 近代雑誌データセット, Geoshapeリポジトリ http://codh.rois.ac.jp/dataset/

・J-TEXT 日本文学電子図書館 http://www.j-texts.com/index.html

・浄土宗全書検索システム http://www.jozensearch.jp/pc/

国立国語研究所 https://www.ninjal.ac.jp/database/

国際日本文化研究センター http://db.nichibun.ac.jp/ja/

東京大学史料編纂所 http://wwwap.hi.u-tokyo.ac.jp/ships/db.html

国立歴史民俗博物館 『データベースれきはく』 https://www.rekihaku.ac.jp/doc/t-db-index.html

・データベース『日本列島の旧石器時代遺跡』 http://palaeolithic.jp/data/index.htm

・身装画像データベース<近代日本の身装文化> http://htq.minpaku.ac.jp/databases/mcd/shinsou.html

・怪異・妖怪伝承データベース http://www.nichibun.ac.jp/youkaidb/

・怪異・妖怪画像データベース http://www.nichibun.ac.jp/YoukaiGazouMenu/

まんが日本昔ばなし〜データベース〜 http://nihon.syoukoukai.com/

・少年犯罪データベース http://kangaeru.s59.xrea.com/

名探偵コナン 全事件レポート編纂室 https://websunday.net/conandb/top.html

・だじゃれステーション https://dajare.jp/

・楽曲動画印象データセット配布サイト https://nkmr.io/mood/index.html

国立感染症研究所 IDWR速報データ https://www.niid.go.jp/niid/ja/data.html

・NDBオープンデータ https://www.mhlw.go.jp/stf/seisakunitsuite/bunya/0000177182.html

・日本蛋白質構造データバンク https://pdbj.org/

・国際ヒトエピゲノムコンソーシアム 日本チームのデータベース http://crest-ihec.jp/database/index.html

・CEBESゲノムデータベース http://www.nies.go.jp/genome/genome_list.html

・組織学カラースライドデータベース http://db.kobegakuin.ac.jp/kaibo/his_pp/Japanese.html

・侵入生物データベース https://www.nies.go.jp/biodiversity/invasive/

・日本産アリ類画像データベース http://ant.miyakyo-u.ac.jp/J/index.html

・イネ品種・特性データベース https://ineweb.narcc.affrc.go.jp/

・橋梁年鑑データベース http://www.jasbc.or.jp/kyoryodb/index.cgi

・大島てる 事故物件 http://www.oshimaland.co.jp/

・ラーメンデータベース https://ramendb.supleks.jp/

【記憶の検索キーワード】

データベース, データセット

グラフ構成の原則の学習記録

【学習動機】

"The Elements of Graphing Data"の第2章では, グラフを構成する際の原則が, グラフの例を挙げながら紹介されている.

https://www.amazon.co.jp/Elements-Graphing-Data-William-Cleveland/dp/0963488414

そこの原則を表す文をメモしておく.

【学習内容】

用語の定義

Fig 2.1とFig 2.2 を引用. これらの用語を知っておくと, 可視化ライブラリの関数の引数の意味が分かってきてよい気がする.

f:id:moratoriamuo271:20190118125451p:plain

f:id:moratoriamuo271:20190118125516p:plain

Clear Vision

  • Make the data stand out. Avoid superfluity. (データを目立たせ, 過剰は避けろ.)

  • Use visually prominent graphical elements to show the data. (視覚的に目立ったグラフの要素を使え.)

  • Use a pair of scale lines for each variable. Make the data rectangle slightly smaller than the scale-line rectangle. Tick marks should point outward. (上下左右に目盛線を.データの箱は目盛線の箱より僅かに小さくすべき. 目盛は外側に.)

  • Do not clutter the interior of the scale-line rectangle. (目盛線の箱の内側を散らかすな.)

  • Do not overdo the number of tick marks. (目盛を使いすぎるな.3~10個で十分.)

  • Use a reference line when there is an important value that must be seen across the entire graph, but do not let the line interfere with the data. (グラフ全体を横切らなければならない重要な値があるなら基準線を使え.ただし,データの邪魔はするな.)

  • Do not allow data labels in the interior of the scale-line rectangle to interfere with the quantative data or to clutter the graph. (量的データを邪魔したり,グラフを散らかすようなデータラベルを目盛の箱の内側に置くな.)

  • Avoid putting notes and keys inside the scale-line rectangle. Put a key outside, and put notes in the caption or in the next. (注釈や凡例を目盛線の箱の内側に置くのは避けよ. 凡例は外側に. 注釈は説明文か文中に.)

  • Overlapping plotting symbols must be visually distinguishable. (重なり合うプロットは見分けがつくようにすべし.)

  • Superposed data sets must be readily visually assembled. (重ね合わせるデータセットは視覚的に分かりやすく整理すべし.)

  • Visual clarity must be preserved under reduction and reproduction. (複製や縮小を経ても視覚的明確さが損なわれないようにすべし.)

Clear Understanding

  • Put major conclusions into graphical form. Make captions comprehensive and informative. (主要な結論をグラフの形式に落とし込むこと.説明文は包括的で有益に.)

  • Error bars should be clearly explained.(エラーバーは明確に説明されるべし.)

  • When logarithms of a variable are graphed, the scale label should correspond to the tick mark labels. (対数グラフを書くときは, スケールラベルを目盛ラベルと対応させるべし.)

  • Proofread graphs. (グラフは校正せよ.)

  • Strive for clarity. (今までの原則, 全部やれ.)

【学習予定】

色調が与える印象などについても学んでいきたいところ.