ニコニコ動画への投稿数を時系列解析する。
どんなときに投稿数が増減するか、またその予測をする。
事後(というか事前にもわかっていたけど)ボカロの誕生日に投稿数が跳ね上がっていた。
解析結果を見た先輩が
「この2年ほどで初音ミクの誕生日(8月31日)に投稿数も変化点スコアも低下傾向にあるから、オワコン化してきてるんじゃないの?」
と指摘され激怒した。ただ、誕生日での観測に限っていうと、投稿数はなんとなくシグモイドカーブに乗っていて、が2年くらいな気がする。
変化点検出には昔動かせなくて悔しかった、ChangeAnomalyDetectionを使う。
ベイズ変化点検出bcpパッケージも使う。
時系列データでは、データがランダムウォークしているか否か、が重要で、Phillips–Perron testで検定できる。
今回のデータは p=0.01 で、帰無仮説：ランダムウォークである は棄却された。
# 投稿日時の処理 hmiku <- read.delim("hmiku.txt", header=TRUE) submit_date <- as.Date(hmiku$post) birth <- as.Date("07-08-31") # 初音ミクの誕生日 cut_month <- seq(birth, max(submit_date), by="month") # 毎月 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) # plot births <- vector("list", 3) # 初音ミク・鏡音リンレン・巡音ルカの誕生日に相当する日 births[[1]] <- mapply(function(x) grep(paste(x, "-08-31", sep=""), names(music_count)), 8:12) births[[2]] <- mapply(function(x) grep(paste(x, "-12-27", sep=""), names(music_count)), 8:12) births[[3]] <- mapply(function(x) grep(paste(x, "-01-30", sep=""), names(music_count)), 9:12) names(births) <- c("Hatsune", "Kagamine", "Megurine") par(mfrow=c(2, 1), mar=c(5, 4, 4, 2)+0.1) plot(music_count, xlab="", ylab="Posted musics", type="l", xaxt="n") axis(1, at=unlist(mapply(function(x) grep(x, names(music_count)), cut_month)), label=tail(cut_month, -1), las=2) title("Posted musics on Nico Nico video every day") 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("topright", legend=c("time series", "Anomaly detection score", "Bayesian posterior mean"), lty=1:3, col=1:3, bty="n", lwd=3, cex=1.2) points(bcp0$posterior.mean, type="l", lty=2, col=3, lwd=2) par(mar=c(5.1, 4.1, 0, 2.1), bty="l") plot(bcp0$posterior.prob, type="l", xlab="", ylab="Bayesian posterior probability", xaxt="n") axis(1, at=unlist(mapply(function(x) grep(x, names(music_count)), cut_month)), label=tail(cut_month, -1), las=2) mapply(function(x) points(births[[x]], bcp0$posterior.prob[births[[x]]], pch=16, col=c(3, 7, 6)[x], cex=1.5), seq(births)) legend("topright", legend=names(births), pch=16, col=c(3, 7, 6), bty="n", cex=1.5)
PP.test(as.ts(music_count)) Phillips-Perron Unit Root Test data: as.ts(music_count) Dickey-Fuller = -33.9685, Truncation lag parameter = 8, p-value = 0.01
