初音ミクの投稿数をバースト解析

初音ミクの投稿数の時系列解析とバースト解析を組み合わせただけ。
前回は19943曲だったが、今回はVOCALONOBISというサイトから2013年5月22付で94334曲のIDを頂いたので、これでやってみる。
wgetの段階で94322曲になった。
　
直近1年間くらいでオワコン化が進んできているのではないかという指摘があったが、このデータで見ると投稿数はじわじわと増加傾向で、しかもこの直近半年ほどでちょっとだけバーストが起きているようだ。
デビュー直後はみんながこぞって購入、投稿し始めたのは想像に難くないが、この直近半年の動きはなんなんだろう…

data1 <- read.delim("20130522count_info.txt")
submit_date <- as.Date(data1[,5])
music_count <- table(submit_date)

library(ChangeAnomalyDetection)
cad <- changeAnomalyDetection(music_count)
cad0 <- cad / max(cad) * max(music_count) # plot 用に値を補正

library(bcp)
bcp0 <- bcp(c(music_count), w0 = 0.2, p0 = 0.2)
plot(bcp0)

library(bursts)
kl0 <- kleinberg(cumsum(table(music_count)))

# プロット用のいろいろ
births <- vector("list", 3) # 初音ミク・鏡音リンレン・巡音ルカの誕生日に相当する日
births[[1]] <- mapply(function(x) grep(paste(x, "-08-31", sep=""), names(music_count)), 7:12)
births[[2]] <- mapply(function(x) grep(paste(x, "-12-27", sep=""), names(music_count)), 7:12)
births[[3]] <- mapply(function(x) grep(paste(x, "-01-30", sep=""), names(music_count)), 9:12)
names(births) <- c("Hatsune", "Kagamine", "Megurine")

# プロットを始める
par(mar=c(6.5, 4, 4, 2))
plot(music_count, xlab="", ylab="Posted musics", type="l", xaxt="n", ylim=c(0, 200))
axis(1, at=unlist(mapply(function(x) grep(x, names(music_count)), cut_month)), label=cut_month, las=2)
title(paste("Posted musics on Nico Nico video every day", paste(range(submit_date), collapse=" ~ ")), cex.main=2)
mapply(function(x) points(births[[x]], music_count[births[[x]]], pch=16, col=c(3, 7, 6)[x], cex=1.5), seq(births))
points(cad0, lty=2, col=2, type="l", lwd=2)
legend("topleft", legend=c("time series", "Anomaly detection score", "Bayesian posterior mean"), lty=1:3, col=1:3, bty="n", lwd=3, cex=1.4)
legend(par()$usr[1], 160, legend=paste(names(births), "birthday"), pch=16, col=c(3, 7, 6), bty="n", cex=1.5)
legend("topright", "Bayesian posterior probability > 0.9", col=4, pch="↓", pt.cex=2, bty="n")
points(bcp0$posterior.mean, type="l", lty=2, col=3, lwd=3)

postprob <- bcp0$posterior.prob # 事後確率が 0.9 以上のものだけプロットする
#segments(which(postprob > 0.9), 0, y1=-10, lwd=2, col=3)
arrows(which(postprob > 0.9), music_count[postprob > 0.9]+25, y1=music_count[postprob > 0.9]+10, lwd=3, col=4, length=0.1)
par(xpd=TRUE)
segments(subset(kl0, level  == 2)$start, -50, subset(kl0, level  == 2)$end, lwd=3, col=2)
text(rowMeans(subset(kl0, level==2)[, c("start", "end")]), -50, "Burst", cex=1.5, adj=c(NA, 1.5))

ファイル処理

wd = "/20130522/"
import codecs
import unicodedata
import re
from progressbar import *
import time
info = [u"再生:<strong>", u"コメント:<strong>", u"マイリスト:<strong>", u"</strong> 投稿", u"<title>"]
res = [""] * len(info)
date = re.compile("\d{2}/\d{2}/\d{2}") # 正規表現
digit = re.compile("[,\d]+") # 正規表現
title = re.compile("<title>.*</title>")

IDs = [f0.rstrip().split("/")[-1] for f0 in open("20130522list.txt", "rU")]
w0 = codecs.open("20130522count_info1.txt", "w", "utf-8")

widgets = ["progress:", Percentage(), Bar()]
maxval = len(IDs)
pbar = ProgressBar(maxval=maxval, widgets=widgets).start()
for f1 in IDs:
	pbar.update(pbar.currval + 1)
	if len(f1) > 0:
		g0 = codecs.open(wd + f1, "rU", "utf-8")
		for tmp in g0:
			tmp = unicodedata.normalize('NFKC', tmp)
			infoTF = map(lambda x: x in tmp, info)
			if any(infoTF):
				if infoTF.index(True) == 3:
					m0 = date.search(tmp)
					res[infoTF.index(True)] = m0.group(0)
				elif infoTF.index(True) == 4:
					m0 = title.search(tmp)
					res[infoTF.index(True)] = m0.group(0)[7:-8]
				else:
					m0 = digit.search(tmp)
					res[infoTF.index(True)] = "".join(m0.group(0).split(","))
		
		w0.write("\t".join([f1] + res) + "\n")

w0.close()