Stanコードを書くとはてブ数は減る?
こんな話がある。
記事内にStanコード書くと、はてブ数が減るって知ってます？
— ホクソエム (@berobero11) 2015, 5月 15
そもそも はてブ受け気にしてないでしょ、、とも思うが、この方のブログのコードにはいつも大変お世話になっているので、今後コードが貼られないようになるととても困る。そこで、事実をデータから確かめることにした。
記事一覧の取得
まずは、この方のブログからコンテンツ一覧とURLを取得する。ご丁寧に記事の目次を作ってくださっているので、それを {httr} にパースさせて指定のタグだけを読み込めばよい。
library(httr)
archives <- httr::GET('http://heartruptcy.blog.fc2.com/archives.html')
contents <- httr::content(archives, as = 'parsed')
library(XML)
# entry 一覧中のリンクタグだけを取得
xpath = "//div[@class='entry']/..//a"
berobero_df <- data.frame(title = sapply(XML::getNodeSet(contents, xpath), XML::xmlValue),
url = XML::xpathSApply(contents, xpath, XML::xmlGetAttr, "href"),
stringsAsFactors=FALSE)
library(magrittr)
library(dplyr)
# entry の URL のみにフィルタ
berobero_df %<>%
filter(grepl('entry', url))
library(xtable)
print(xtable::xtable(head(berobero_df)), type = "html", include.rownames = FALSE)| title | url |
|---|---|
| [Stan] kivantiumさんのブログ記事「アニメキャラのバストサイズとPixiv R-18タグ率の関係」の追加解析 | http://heartruptcy.blog.fc2.com/blog-entry-175.html |
| [コンピュータ将棋] 高校生の時に作った傑作詰将棋 | http://heartruptcy.blog.fc2.com/blog-entry-174.html |
| [JAGS,Stan] ノンパラベイズ（ディリクレ過程）の実装 | http://heartruptcy.blog.fc2.com/blog-entry-173.html |
| [Stan] 生存時間分析 - ハザード関数に時間相関の制約を入れる | http://heartruptcy.blog.fc2.com/blog-entry-172.html |
| [Stan] 不等間隔の状態空間モデル | http://heartruptcy.blog.fc2.com/blog-entry-171.html |
| [JAGS] 2次元のCAR model | http://heartruptcy.blog.fc2.com/blog-entry-170.html |
各記事内にプログラムが含まれるかどうかは、記事中の特定のタグから判断できそうだ。Stan のプログラムかどうか？は直接判断できない (ものによってはファイル名が付いているものもあるのだが、、) 。まあ Stan の記事に Stan 以外のプログラムが貼られることはないだろうから気にしないことにする。
library(stringr)
get_body_text <- function(url, xpath = "//div[@class='entry']") {
f <- function(x) {
# x は factor 型になる
archives <- httr::GET(as.character(x))
parsed <- httr::content(archives, as = 'parsed')
body <- sapply(XML::getNodeSet(parsed, xpath), XML::xmlValue)
return(stringr::str_replace_all(body, '\n|\r|<br />', ''))
}
as.character(lapply(url, f))
}
berobero_df %<>%
filter(grepl('entry', url)) %>%
mutate(body_all = get_body_text(url))
# プログラムを含むと思われるタグ
program_tag <- '<!\\[CDATA\\[.*\\]\\]>'
berobero_df %<>%
mutate(has_program = stringr::str_detect(body_all, program_tag),
body = stringr::str_replace_all(body_all, program_tag, ''))
# 上記タグが入っていないがプログラム添付のある記事 (他にもあるかも)
berobero_df[berobero_df$url == 'http://heartruptcy.blog.fc2.com/blog-entry-166.html', 'has_program'] <- TRUE
print(xtable::xtable(head(berobero_df[c('title', 'url', 'has_program')])),
type = "html", include.rownames = FALSE)| title | url | has_program |
|---|---|---|
| [Stan] kivantiumさんのブログ記事「アニメキャラのバストサイズとPixiv R-18タグ率の関係」の追加解析 | http://heartruptcy.blog.fc2.com/blog-entry-175.html | TRUE |
| [コンピュータ将棋] 高校生の時に作った傑作詰将棋 | http://heartruptcy.blog.fc2.com/blog-entry-174.html | FALSE |
| [JAGS,Stan] ノンパラベイズ（ディリクレ過程）の実装 | http://heartruptcy.blog.fc2.com/blog-entry-173.html | TRUE |
| [Stan] 生存時間分析 - ハザード関数に時間相関の制約を入れる | http://heartruptcy.blog.fc2.com/blog-entry-172.html | TRUE |
| [Stan] 不等間隔の状態空間モデル | http://heartruptcy.blog.fc2.com/blog-entry-171.html | TRUE |
| [JAGS] 2次元のCAR model | http://heartruptcy.blog.fc2.com/blog-entry-170.html | TRUE |
はてブ数の取得
こちらのパッケージを利用させていただく。
# library(devtools)
# install_github('dichika/jaguchi')
# install_github('dichika/hatenab')
library(jaguchi)
hateb <- function(x) {
jaguchi('hatenab', pkgcheck=TRUE, x)[['count']]
}
# 複数URLではソート順が変わるため、1 URL ずつ取得する
berobero_df['hatenab'] <- unlist(lapply(berobero_df[['url']], hateb))
print(xtable::xtable(head(berobero_df[c('title', 'has_program', 'hatenab')])),
type = "html", include.rownames = FALSE)| title | has_program | hatenab |
|---|---|---|
| [Stan] kivantiumさんのブログ記事「アニメキャラのバストサイズとPixiv R-18タグ率の関係」の追加解析 | TRUE | 34 |
| [コンピュータ将棋] 高校生の時に作った傑作詰将棋 | FALSE | 0 |
| [JAGS,Stan] ノンパラベイズ（ディリクレ過程）の実装 | TRUE | 1 |
| [Stan] 生存時間分析 - ハザード関数に時間相関の制約を入れる | TRUE | 2 |
| [Stan] 不等間隔の状態空間モデル | TRUE | 4 |
| [JAGS] 2次元のCAR model | TRUE | 1 |
確認
まずはプログラムを含む記事かどうかで分布をみる。プログラムがあるからブクマが減る、とは言えなさそうだ。一部 プログラムなしでブクマ数が多い記事があるため、これらが印象に影響を与えているのかもしれない。
library(ggplot2)
ggplot(berobero_df) + geom_violin(aes(x = has_program, y = hatenab))
Stan を含む記事だけを抽出して、同様の比較を行う。ほとんどの記事 (20記事中19記事) にコードが貼ってあるためカテゴリ内での比較には意味ない。上と比べると、Stan についてプログラムがある記事のほうが、分布でみてブクマされている傾向があるように見える。
# タイトルに "Stan" を含むもののみ
berobero_stan <- berobero_df %>%
filter(grepl('Stan', title)) %>%
filter(!grepl('Standard', title))
ggplot(berobero_stan) + geom_violin(aes(x = has_program, y = hatenab))
結論
気のせいでは。
おまけ
また、Stan の記事の中でのブクマ数の差が印象の差に繋がっているのかな？と思い、記事内のトピックでさらにカテゴリわけし、比較ができないかを検討する。
辞書は neologd を使っている。neologd の設定はこちらから。
http://d.hatena.ne.jp/dichika/20150328/p1
library(RMeCab)
parsed <- RMeCab::docMatrixDF(berobero_df[, 'body'], pos = c("名詞"))## to make data frameparsed <- parsed[!str_detect(rownames(parsed), '.*[\\.\\*\\+\\(\\)=;%\\#（）].*'), ]
parsed <- parsed[!str_detect(rownames(parsed), '.*[0-9]+.*'), ]
# 1回しか出現しない Term を除外
parsed <- parsed[!(rowSums(parsed) <= 1), ]
# ストップワードを選定すると膨大になるので、
# パースされた単語の中から特徴になりそうな単語を適当に選ぶ
not_stopwords <- c("AIC", "AR", "AWS", "Amazon", "BIC", "BONANZA", "BUGS",
"Bootstrap", "CAR", "CMDSTAN", "CPU", "CUDA", "CmdStan", "Dirichlet", "EMアルゴリズム",
"Fisher", "GLM", "GLMM", "GPS", "GPSShogi", "GPSfish", "GPS将棋", "GPU", "Gaussian",
"GeoBUGS", "Gibbs", "Gurobi", "HMC", "HMM", "Hamiltonian", "Hierarchical", "IJulia",
"JAGS", "Julia", "LDA", "Lasso","MCMC", "NUTS", "OpenBUGS", "PCA", "PRML", "Poisson",
"Ponanza", "PyMC", "Python", "R", "RSRuby", "RStan", "Rhat", "Rtools", "Ruby", "SEM",
"STAN", "SVM", "Sampler", "Sampling", "Stan", "Theano", "Unigram", "WAIC", "WBIC",
"WinBUGS", "bernoulli", "binomial", "cholesky", "poisson", "rstan",
"カルマンフィルタ", "カーネル", "ガウス分布", "ガウス過程",
"ガンマ関数", "クラスタリング", "グラフ", "グラフィカルモデル", "グリーン関数",
"ゲーム理論", "コレスキー分解", "コンピュータ将棋", "コーシー分布", "サンプリング", "シミュレーション",
"スクリプト", "スクレイピング", "スケーリング", "ディリクレ", "ディリクレ分布",
"データセット", "ノンパラベイズ", "ハイパーパラメータ",
"ハザード", "バイオ", "バイオインフォ", "パターン認識", "パラメータ",
"ヒストグラム", "フィッティング", "フィルタ", "フィルター",
"ベイジアン", "ベイジアンネットワーク", "ベイズ", "ベイズ推定",
"モデリング", "モンテカルロ", "モンテカルロ法", "ヤコビアン", "ランダム",
"ロジスティック回帰", "一様分布", "主成分分析", "事前", "事後", "事後確率", "二項分布",
"信頼区間", "個体差", "傾き", "共分散構造分析", "共役", "分子生物学", "分散共分散行列",
"分配関数", "切片", "制約", "勝ち", "勝敗", "勝負", "同時分布", "名人",
"周期", "回帰", "因果", "因果関係", "固有値", "多価関数", "多項分布", "対局", "対数", "対数正規分布",
"将棋", "将棋電王戦", "常微分方程式", "微分方程式", "投了", "持ち時間", "持ち駒", "持将棋",
"数値計算", "数学", "時系列", "時間発展", "最尤推定", "最尤法", "最適化", "有向グラフ", "有意",
"構造", "標本", "標準偏差", "機械学習", "次元", "正規分布", "汎関数", "混合", "混合比",
"状態空間", "生データ", "生存", "相関係数", "確率", "確率変数", "確率過程",
"統計", "統計学", "統計的", "自己相関", "自由エネルギー", "英語", "評価関数", "認識",
"近似", "近傍", "逆関数", "過程", "量子力学", "階層", "集合", "電王戦", "非線形", "項目応答理論")
parsed <- parsed[not_stopwords, ]
# 1回しか出現しない Term を除外
parsed <- parsed[!(rowSums(parsed) <= 1), ]
tfidf <- function(data) {
tf <- data
idf <- log(nrow(data) / (colSums(data >= 1)))
result <- data
for(word in names(idf)){
result[, word] <- tf[, word] * idf[word]
}
result
}
tf <- tfidf(parsed)
tf[is.na(tf)] <- 0
important <- names(sort(apply(tf, 1, mean), decreasing = TRUE))[1:150]
parsed <- parsed[important, ]
# Term を含まない Document を除外
parsed <- parsed[, !(colSums(parsed) == 0)]
library(slam)
stm <- slam::as.simple_triplet_matrix(t(parsed))
library(tm)
dtm <- tm::as.DocumentTermMatrix(stm, weighting = tm::weightTf)トピック数を変えながらいくつか試してみたが、あまりうまく分離できているようには見えない、、。だいたいがモデリングの話をされているので、その中でトピック分けするのは工夫が必要だ。また、目視でみても特定のトピックにブクマが集中している、ということはなさそうなので、あまりよいやり方ではなさそうだ。
library(topicmodels)
lda <- topicmodels::LDA(dtm, 6)
terms(lda, 10)## Topic 1 Topic 2 Topic 3 Topic 4 Topic 5
## [1,] "電王戦" "Stan" "確率" "Stan" "BUGS"
## [2,] "コンピュータ将棋" "モデリング" "ベイズ" "CAR" "R"
## [3,] "GPS将棋" "統計" "WAIC" "R" "ベイズ"
## [4,] "将棋" "微分方程式" "正規分布" "ガウス過程" "WinBUGS"
## [5,] "勝負" "BUGS" "事後" "JAGS" "階層"
## [6,] "Ponanza" "R" "WBIC" "ベイズ" "時系列"
## [7,] "GPSfish" "パラメータ" "確率変数" "回帰" "AR"
## [8,] "GPS" "回帰" "事前" "パラメータ" "確率"
## [9,] "対局" "RStan" "近似" "状態空間" "勝敗"
## [10,] "名人" "確率過程" "構造" "構造" "統計"
## Topic 6
## [1,] "LDA"
## [2,] "Stan"
## [3,] "パラメータ"
## [4,] "確率"
## [5,] "次元"
## [6,] "MCMC"
## [7,] "グラフィカルモデル"
## [8,] "ハイパーパラメータ"
## [9,] "認識"
## [10,] "R"