DTM解析

DTMを行うためにコメントを文章化して lexicalize し、データ加工する。
Rに移る。w コメントは異常に多いため、w が2つ以上並ぶものは w 1つに処理する。
1つの動画に付いているコメントをつなぎあわせて文章化する。コメント 1000 以上ある動画を対象にした。というのも、10万も動画があったら解析が終わらないので、適当にコメント数を設定したら 1000 以上で 10560動画(約10%)になったからこれでいっかって。

wd <- "/niconico/output/"

### あるコメント数を超えている動画だけ対象とする
metafile <- read.delim(paste(wd, "meta.txt", sep=""))
idx <- which(metafile$comment_counter > 1000)
commentfile <- paste(metafile$video_id[idx], ".txt", sep="")

library(foreach)
library(doSNOW)
cl <- makeCluster(12, "SOCK")
registerDoSNOW(cl)
doc1 <- foreach(i = seq(commentfile), .combine=c) %dopar% {
	tmp <- read.delim(paste(wd, "mecab/", commentfile[i], sep=""))
	paste(gsub("w{2,}", "w", tmp$comment), collapse=" ") # w が2個以上連続するものは1個にする
}

lexicalize するが、

library(lda)
lex0 <- lexicalize(doc1) # 5h

だと6時間くらいかかるので並列化すると

v0 <- unique(unlist(strsplit(doc1, " ")))
lexvec <- foreach(i = seq(doc1), .combine=c) %dopar% {
	t0 <- table(strsplit(doc1[i], " "))
	t1 <- rbind(match(names(t0), v0) - 1, t0)
	c(t1)
}
uw <- foreach(i = seq(doc1), .combine=c) %dopar% {
	t0 <- length(unique(strsplit(doc1[i], " ")[[1]]))
}
s0 <- unlist(mapply(rep, seq(uw), uw*2))
lex0$documents <- list(documents=lapply(split(lexvec, s0), matrix, nr=2), vocab=v0)

このまま解析すると、あまりにも登場頻度の高すぎる単語や、たった一度しか使われない単語などのせいでゴミみたいなのが多くなってしまうので、TFIDFっぽいスコアを計算しようと思ったのだが、これでいくと「かわいい」が消えてしまうので、登場頻度が30以上の単語を使用することにした。

library(foreach)
library(doSNOW)
library(utils)
library(iterators)
library(doParallel)
cl <- makeCluster(10, "SOCK")
registerDoSNOW(cl)
# なんかよくわからないが 294696 番目の "ぁぁぁぁって" が消える

wordset <- mapply(function(x) x[1,], lex0$documents) # documents中の単語リスト
allfreq <- matrix(unlist(lex0$documents), nr=2)
wordfreq1 <- tapply(allfreq[2,], allfreq[1,], sum) # すべての単語の、全documents中の出現頻度
wordfreq <- replace(rep(1, length(v0)), as.numeric(names(wordfreq1)) + 1, wordfreq1)
res <- matrix(0, nr=length(wordfreq), nc=4)
dimnames(res) <- list(v0, c("documents", "count", "freq", "score"))
res[, "documents"] <- length(lex0$documents)

mn <- 1000 # 一度に並列計算するとCPUがなぜか止まるので分割
mn0 <- ceiling(length(lex0$vocab)/mn)
pb <- txtProgressBar(min=1, max=mn0, style=3)
for(i in 1:(mn0-1)){
	setTxtProgressBar(pb, i)
	res[mn*(i-1)+seq(mn), "freq"] <- foreach(v = mn*(i-1)+seq(mn), .combine=c) %dopar% { # vocab と i は 1 ずれているので注意
		count_docs <- sum(sapply(lapply(wordset, "==", v-1), any)) # その単語が出現する文章の数
	}
}
rest <- (mn*(mn0-1)):length(lex0$vocab)
res[rest, "freq"] <- foreach(v = rest, .combine=c) %dopar% { # vocab と i は 1 ずれているので注意
	count_docs <- sum(sapply(lapply(wordset, "==", v-1), any)) # その単語が出現する文章の数
}

res[, "count"] <- wordfreq
res[, "score"] <- log(res[, "count"]) * log(res[, "documents"]/res[, "freq"])
res <- as.data.frame(res)

DTMで解析できるようにファイルを作る。
ある期間内の動画数が適当な数はあるように時系列を設定する。
今回は14区間、約4~5ヶ月分くらいに設定すると、1区間あたり600〜1000動画になった。

# upload time
up <-  as.Date(substr(as.character(metafile$upload_time[idx]), 1, 10))
time5 <- seq(min(up), max(up), length=15) # 4~5 ヶ月くらい
ordidx <- order(up)
timeidx <- cut(up[ordidx], time5, right=FALSE, include.lowest=TRUE)

これを時系列データ seq.dat として保存する。

write.table(c(length(table(timeidx)), table(timeidx)), "niconico-seq.dat", row.names=FALSE, col.names=FALSE, quote=FALSE)

文章の品詞分解データをDTM用に mult.dat とする。

dat <- mapply(function(x) paste(uw1[x], paste(lex1$documents[[x]][1,], lex1$documents[[x]][2,], sep=":", collapse=" ")), ordidx)
# write.table(dat, "niconico-mult.dat", row.names=FALSE, col.names=FALSE, quote=FALSE)

DTM用のスクリプト ./main があるディレクトリで以下の様に実行した。
入力ディレクトリは corpus_prefix, 出力ディレクトリは outname を適宜変更する。
トピック数 ntopics は20、サンプリングのパラメータ\alphaや収束回数などはデフォルトのままとした。

# Run the dynamic topic model.
./main \
  --ntopics=20 \
  --mode=fit \
  --rng_seed=0 \
  --initialize_lda=true \
  --corpus_prefix=/niconico/output/dtm/niconico \
  --outname=/niconico/output/dtm/dtm_output/ \
  --top_chain_var=0.005 \
  --alpha=0.01 \
  --lda_sequence_min_iter=6 \
  --lda_sequence_max_iter=20 \
  --lda_max_em_iter=10

各時系列でのトピックの存在割合を示す。

# topic distribution
theta <- matrix(scan(paste(dtmwd, "gam.dat", sep="")), nc=20, byrow=TRUE)
theta <- theta / rowSums(theta)
ntime <-  table(timeidx)
theta0 <- apply(theta, 2, tapply, tidx, sum)/ntime
pidx <- c(2, 5:6, 10:12, 7) # わかりやすいやつ

# 原稿
par(mar=c(5, 4, 4, 4))
matplot(theta0[, pidx], type="l", lty=seq(pidx), col=1, xaxt="n", lwd=3, ylab="probability")
#axis(1, seq(nrow(theta0)), head(time5, -1), las=2, srt=45)
axis(1, seq(nrow(theta0)), label=NA)
text(seq(nrow(theta0)), par()$usr[3], head(time5, -1), xpd=TRUE, adj=c(1.1, 1.6), srt=45)
text(par()$usr[2], theta0[14, pidx], paste("Topic", pidx), xpd=TRUE, pos=4)
title("Topic distribution probability")

# 他のトピック
pidx <- c(1,3:4,8,16:17,20)
par(mar=c(5, 4, 4, 4))
matplot(theta0[, pidx], type="l", lty=seq(pidx), col=seq(pidx), xaxt="n", lwd=4, ylab="probability")
#axis(1, seq(nrow(theta0)), head(time5, -1), las=2, srt=45)
axis(1, seq(nrow(theta0)), label=NA)
text(seq(nrow(theta0)), par()$usr[3], head(time5, -1), xpd=TRUE, adj=c(1.1, 1.6), srt=45)
ty <- theta0[14, pidx] +c(-0.0025,0,0,0,0,0,0)
text(par()$usr[2], ty, paste("Topic", pidx), xpd=TRUE, pos=4)
title("Topic distribution probability")

タグ数解析を行う。

# tag 数
cv <- as.character(unlist(read.csv("cvlist.txt", header=FALSE)))
tagtable <- table(factor(unlist(strsplit(as.character(metafile$tag), ",")), cv))
# 期間ごとの tag 数
up1 <-  as.Date(substr(as.character(metafile$upload_time), 1, 10))
time5 <- seq(min(up), max(up), length=15) # 4~5 ヶ月くらい
ordidx1 <- order(up1)
cut0 <- cut(up1, time5)

tmp <- mapply(function(x) table(factor(x, cv)), lapply(split(strsplit(as.character(metafile$tag), ","), cut0), unlist))
tmp[order(rowSums(tmp), decreasing=TRUE),]
tmp1 <- sweep(tmp, 2, colSums(tmp), "/")

idx2 <- head(order(apply(apply(tmp1, 1, range), 2, diff), decreasing=TRUE), 5)
par(mar=c(5, 4, 4, 5))
matplot(t(tmp1[idx2,]), type="l", ylab="proportion in tags", col=1, lwd=3, xaxt="n")
axis(1, seq(nrow(theta0)), label=NA)
text(seq(nrow(theta0)), par()$usr[3], head(time5, -1), xpd=TRUE, adj=c(1.1, 1.6), srt=45)
texty <- tmp1[idx2,14] + c(0, 0, 0, -0.002, 0)
text(par()$usr[2], texty, rownames(tmp1)[idx2], xpd=TRUE, pos=4)
title("Dynamic change of the proportion of CVs in tags")

# 伸びてきた声優
idx3 <- apply(tmp1[,1:5]<0.01, 1, all) & apply(tmp1[,10:14]>0.01, 1, all)
idx3[c(13,59,152)] <- FALSE # キタエリはプロットしない
par(mar=c(5, 4, 4, 5))
lty0 <- c(2, 1, 5, 4, 3)
tmp2 <- t(tmp1[idx3,])
matplot(tmp2, type="l", ylab="proportion in tags", col=1, lwd=3, lty=lty0, xaxt="n")
axis(1, seq(nrow(theta0)), label=NA)
text(seq(nrow(theta0)), par()$usr[3], head(time5, -1), xpd=TRUE, adj=c(1.1, 1.6), srt=45)
texty <- tmp1[idx3,14] + c(0, 0.0005, 0, 0, -0.0005)
text(par()$usr[2], texty, rownames(tmp1)[idx3], xpd=TRUE, pos=4)
title("Dynamic change of the proportion of CVs in middle period")