この記事はごちうさ住民 Advent Calendar 2015 の第26日目の記事です(2年連続の勝手に参加。
ネタと解析がなやましい これが私の最大限です
11月の始めころになんとなしにググっていたところ、上記のアドベントカレンダーが今年も行われることに気づいた。
ごちうさに関連してそうなことって何にしようかとネタを考えてると、ふとツイッターの解析でもしてみようと思った。
青い鳥と黄色い蛇
要はツイッターからPython でデータ収集をする話です。
残念なことに現在のツイッターAPI では1週間しかさかのぼってツイート収集ができません。この解析を始めたのが11月半ばだったので、11月9日までしかさかのぼってデータが取れませんでした。
1週間経たないうちにこまめに収集スクリプトを回します。しかしここでも残念ながら、2回目からいきなり忘れていて11月16日あたりのツイートが収集できませんでした。
スクリプトはここからパクりました。1回のAPI で100しか取れず、15分制限があるので #gochiusa もしくは #ごちうさ のハッシュタグを含むツイートを1週間分取得するには1回で4時間くらいかかりました。
技術的には、ツイッターデベロッパーでAPI を使えるようにすることと、日本語のツイート内容も出力していたら改行とか絵文字とかでファイル構造がグダグダになってしまっていたので、ツイート内容は解析しないならば出力に含まないほうがよさそう。一般公開しているユーザー名やGIS データも今後解析に使えそうなので取ります。倫理? なにそれ美味しいの?
#!/usr/bin/python # -*- coding: utf-8 -*- #↑これがないと日本語でのコメントは出来ない import twitter import time import csv CONSUMER_KEY = "あなたの CONSUMER_KEY" CONSUMER_SECRET = "あなたの CONSUMER_SECRET" ACCESS_TOKEN_KEY = "あなたの ACCESS_TOKEN_KEY" ACCESS_TOKEN_SECRET = "あなたの ACCESS_TOKEN_SECRET" api = twitter.Api(consumer_key=CONSUMER_KEY, consumer_secret=CONSUMER_SECRET, access_token_key=ACCESS_TOKEN_KEY, access_token_secret=ACCESS_TOKEN_SECRET) #指定ハッシュタグツイートを検索 maxcount=1000 maxid =0 terms=["#gochiusa" ,u"#ごちうさ"] # ここを好きに変えよう search_str=" OR ".join(terms) tweets = api.GetSearch(term=search_str, count=100, result_type='recent') #今回は ['id', 'date', 'follow', 'follower', 'tweet']の順に格納 wait_time = 15*60/180 + 0.01 rate = api.GetRateLimitStatus() print "Limit %d / %d" % (rate['resources']['search']['/search/tweets']['remaining'],rate['resources']['search']['/search/tweets']['limit']) tm = time.localtime(rate['resources']['search']['/search/tweets']['reset']) print "Reset Time %d:%d" % (tm.tm_hour , tm.tm_min) print "-----------------------------------------\n" fout = open('tweet 内容も含む.txt', 'w') fout0 = open('tweet 内容は含まない.txt', 'w') found = api.GetSearch(term=search_str, count=100, result_type='recent') time.sleep(wait_time) # 引っかからないように待ちます i = 0 maxid = 0 maxcount = 100000000 while True: for f in found: if maxid > f.id or maxid == 0: maxid = f.id i = i + 1 if len(found) == 0: break if maxcount <= i: break print maxid,i found = api.GetSearch(term=search_str, count=100, result_type='recent', max_id=maxid-1) time.sleep(wait_time) for tweet in found: #呟き集を一つずつfor文で回す a0 = [ tweet.user.id, # ユーザーid tweet.user.created_at, # ユーザー登録時刻 tweet.created_at, # 呟き時刻 tweet.user.friends_count, # フォロイー tweet.user.followers_count, # フォロワー tweet.id, # id tweet.user.screen_name, # screen name tweet.user.time_zone, # time zone tweet.text # 呟き内容 ] list_mining = map(lambda x: unicode(x).encode("utf-8"), a0) # csv出力 fout.write("\t".join(list_mining)+"\n") fout0.write("\t".join(list_mining[:-1])+"\n") fout.close() fout0.close()
ビッグデータ処理伝説 data.table
がんばってツイート収集しまくると、1週間あたり20MB くらいでした。これくらいならば普通にR で読み込んでも問題ないですが、容量のでかいデータを読み込むのに最適なdata.table パッケージのfread を使うとかなりストレスが減ります。
RjpWiki 先輩の優雅なPOSIXlt チュートリアル
時間のオブジェクトにはDate とPOSTXlt があります。Date は
"2015-12-26"
のように日付だけですが、POSTXlt は
"2015-12-26 22:30:00 JST"
のように時間とタイムゾーンがつきます。演算は秒単位でできるが、プロットしたときにはたぶんPC 上のカウント時間? になるのでPOSTXlt オブジェクトでx軸をプロットするときには注意したほうがいい。
Date とPOSTXlt の変換は
strptime("2015/12/26 22:30:00", "%Y/%m/%d %H:%M:%S", tz="JST")
format(x, "%Y/%m/%d") format(x, "%H:%M:%S")
でゴリ押しできる。
ひと解析で普通のエラーだと見抜いたよ
なお、POSTXlt への変換のときに、下のスクリプトを一回行わないと変なことになる。
lct <- Sys.getlocale("LC_TIME"); Sys.setlocale("LC_TIME", "C") x <- c("1jan1960", "2jan1960", "31mar1960", "30jul1960") z <- as.Date(x, "%d%b%Y")
時系列データ攻略完了(みっしょんあんこんぷりーてぃっど)
収集開始し始めた日+tweet 内容は含まない.txt みたいなファイル名が複数できるので、これらを読み込みます。
重複した日の重複したツイートがあるっぽいけど、ID がなんかうまくいっていないことと、消されたかなにかで数が変動しているので、各日で最大数収集しているデータを採用することにします。
library(data.table) wd <- "/gochiusaadv/" setwd(wd) txt <- list.files(wd, pattern="tweet 内容は含まない.txt") dat <- mapply(fread, txt, SIMPLIFY=FALSE) # データの一括読み込み dat2 <- mapply(function(x) strsplit(as.matrix(x$V3), " "), dat) # 投稿時間を取得 z <- mapply(function(y) mapply(function(x) strptime(paste(x[c(3, 2, 6, 4)], collapse=""), "%d%b%Y%H:%M:%S", "JST")+ 9*60*60, y, SIMPLIFY=FALSE), dat2) # R のPOSIXt 形式にする s <- mapply(function(x) sapply(x, format, "%y/%m/%d"), z, SIMPLIFY=FALSE) # Date オブジェクト t0 <- lapply(s, table) # 各日のツイート数 # 各ファイルの各日のツイート数にわける lab <- unique(unlist(lapply(t0, names))) m <- matrix(0, length(txt), length(lab), dimnames=list(txt, lab)) for(i in seq(txt)){ m[i, names(t0[[i]])] <- t0[[i]] } # 最大数の集計 res <- NULL for(i in 2:ncol(m)){ j <- which(s[[which.max(m[,i])]]==lab[i]) res <- c(res, z[[which.max(m[,i])]][j]) }
甘えん坊なあのデータはシャボン玉のように儚く消える
12月24日夕方までデータを収集して2147467 ツイートのうち、重複した日をのぞいて1315764 ツイートが解析対象です。
スニーキングストーキングストーカープロット
プロットして可視化できるようにデータを整えていきます。
titles <- c("笑顔とフラッシュがやかましい これが私の自称姉です", "灰色兎と灰かぶり姫", "回転舞踏伝説アヒル隊", "ココア先輩の優雅なお茶会チュートリアル", "ひと口で普通のもちもちだと見抜いたよ", "木組みの街攻略完了(みっしょんこんぷりーと)", "甘えん坊なあの子はシャボン玉のように儚く消える", "スニーキングストーキングストーカーストーリー", "毛玉は特攻し無慈悲なボタンは放たれる", "Eを探す日常", "スターダスト・マイムマイム", "宝物は君の決定的瞬間") titles <- paste(paste("第", seq(titles), "羽", sep=""), titles, sep="\t") d <- seq(as.Date("2015/11/9"), as.Date("2015/12/31"), "day") # 収集し始めた日から vline <- seq(as.Date("2015/10/10"), as.Date("2015/12/26"), "week") # 放送日 ymd <- lapply(res, as.Date) # yyyy/mm/dd 形式に変換 ymd_diff <- sapply(ymd, "-", d[1]) + 1 # d オブジェクトのindex にしたい x <- table(d[ymd_diff]) # 各日のツイート数
ロリコンは特攻し無慈悲な結果は放たれる
データを収集した期間の1日毎のツイート数をプロットします。
1周間単位では、最速放送日の土曜日にピークがあり、その3日前にちいさなピークがあります。
第9羽の放送日(12月5日)の前日が異常に伸びています。これは何かあったのでしょうか?
そうですね、チノちゃんの誕生日ですね。
# 横にベタッとプロット plot(d[match(names(x), as.character(d))], x, type="o", xlab="Day", ylab="Tweets/day #gochiusa", lwd=3, pch=16) axis(2) abline(v=vline) text(vline, par()$usr[4], paste("#", seq(titles) ,sep=""), xpd=TRUE, pos=3)
P(ピーク)を探す日常
1分ごとのツイートに集計する。
s01 <- strptime(c("2015/11/09 00:00:00", "2015/11/16 00:00:59"), "%Y/%m/%d %H:%M:%S") x0 <- head(seq.POSIXt(s01[1], s01[2], by="min"), -1) m <- array(NA, c(60, 24, length(d))) dimnames(m) <- list(0:59, 0:23, as.character(d)) for(i in seq(max(ymd_diff))){ for(j in which(ymd_diff==i)){ hr <- as.numeric(format(res[[j]], "%H")) + 1 minute <- as.numeric(format(res[[j]], "%M")) + 1 m[minute, hr, i] <- sum(m[minute, hr, i], 1, na.rm=TRUE) } }
各週ごとにツイート数の推移をみてみよう。土曜日の日付が変わる瞬間がもっともツイート数が多い。これはTOKYO MX, サンテレビ、KBS京都の22:30-23:00の放送枠である。最速はAT-X の21:30-22:00 だが、このときのツイート数は前述の1/10 程度である。
他、BS11 の水曜1:00-1:30 もちらほら。
cols <- rainbow(7) yl <- c(0, max(m, na.rm=TRUE)) w <- c("Mon", "Tue", "Wed", "Thu", "Fri", "Sat", "Sun") tmpx <- seq(s01[1], by = "day", length.out = 7) labelx <- c(tmpx, tail(tmpx, 1) + 24*60*60) xl <- range(labelx) y0 <- c(m[,,as.character(tmpx)]) plot(x0, y0, type="l", xaxt="n", xlim=xl, ylim=yl, xlab="Weekdays", ylab="Tweets/minute #gochiusa", col=cols[1], lwd=2) abline(v=labelx, lty=3) axis(1, at=labelx, label=NA) text(head(labelx, -1)+12*60*60, par()$usr[3], w, xpd=TRUE, pos=1) for(i in 2:7){ tmpx <- seq(tail(labelx, 1), by = "day", length.out = 7) labelx <- c(tmpx, tail(tmpx, 1) + 24*60*60) y0 <- c(m[,,as.character(tmpx)]) lines(x0, y0, col=cols[i], lwd=2) } legend("topleft", legend=tail(titles, length(cols)), lty=1, col=cols, lwd=3, bg="white")
個別にはこうなる。
par(mfrow=c(7,1) ,mar=c(2, 4, 2, 2)) tmpx <- seq(s01[1], by = "day", length.out = 7) labelx <- labelx0 <- c(tmpx, tail(tmpx, 1) + 24*60*60) xl <- range(labelx) y0 <- c(m[,,as.character(tmpx)]) plot(x0, y0, type="l", xaxt="n", xlim=xl, ylim=yl, xlab="Weekdays", ylab="Tweets/minute #gochiusa", col=cols[1], lwd=2) abline(v=labelx, lty=3) axis(1, at=labelx, label=NA) text(head(labelx, -1)+12*60*60, par()$usr[3], w, xpd=TRUE, pos=1) title(tail(titles, length(cols))[1]) for(i in 2:7){ tmpx <- seq(tail(labelx, 1), by = "day", length.out = 7) labelx <- c(tmpx, tail(tmpx, 1) + 24*60*60) y0 <- c(m[,,as.character(tmpx)]) plot(x0, y0, type="l", xaxt="n", xlim=xl, ylim=yl, xlab="Weekdays", ylab="Tweets/minute #gochiusa", col=cols[i], lwd=2) abline(v=labelx0, lty=3) axis(1, at=labelx, label=NA) text(head(labelx0, -1)+12*60*60, par()$usr[3], w, xpd=TRUE, pos=1) title(tail(titles, length(cols))[i]) abline(v=labelx0, lty=3) axis(1, at=labelx0, label=NA) text(head(labelx0, -1)+12*60*60, par()$usr[3], w, xpd=TRUE, pos=1) }
データダスト・マイムマイム
22:30-23:00の放送枠を中心にみていこう。
うむ、なるほどわからん。が、やはりAT-X の21:30-22:00 の放送枠はあまりツイートが多くなく、22:30 直前からツイートがすごい速度であがり、放送中はなんやかんやつぶやき、22:27 くらいから放送が終わるので急速に減っている。
tmpx <- seq(s01[1], by = "day", length.out = 7) labelx <- c(tmpx, tail(tmpx, 1) + 24*60*60) xl <- labelx[weekdays(labelx)=="Saturday"]+60*60*22.5 + c(-1, 1)*60*60*1 y0 <- c(m[,,as.character(tmpx)]) plot(x0, y0, type="l", xaxt="n", xlim=xl, ylim=yl, xlab="Weekdays", ylab="Tweets/minute #gochiusa", col=cols[1], lwd=2) axis(1, at=labelx, label=NA) text(head(labelx, -1)+12*60*60, par()$usr[3], w, xpd=TRUE, pos=1) for(i in 2:7){ tmpx <- seq(tail(labelx, 1), by = "day", length.out = 7) labelx <- c(tmpx, tail(tmpx, 1) + 24*60*60) y0 <- c(m[,,as.character(tmpx)]) lines(x0, y0, col=cols[i], lwd=2) } abline(v=seq(xl[1], xl[2], by="5 min"), lty=3) # 5分ごとに破線 t30 <- seq(xl[1], xl[2], by="30 min") text(t30, par()$usr[3], format(t30, "%H:%M"), pos=1, xpd=TRUE) # 30分ごとにテキスト legend("topleft", legend=tail(titles, length(cols)), lty=1, col=cols, lwd=3, bg="white")
時系列解析を行ってツイート数の真の推移を推定して、実際にもっとも盛り上がっていたシーンはなんだったのかを検証しよう。
真面目にカルマンフィルターを用いたらまったくいい感じには推定できなかったので、適当なsplineで補完する。
もっとも盛り上がっていそうなのは、第7羽の22:50 あたりのシーン。
いま見返しているけれどもたぶんここらあたりじゃないのだろうか。
やっぱりチノちゃんだね！(ロリコン乙
# 時系列解析 # 10 秒ごとにする s01 <- strptime(c("2015/11/09 22:20:00", "2015/11/09 23:10:00"), "%Y/%m/%d %H:%M:%S") x0 <- head(seq.POSIXt(s01[1], s01[2], by="10 sec"), -1) x1 <- format(x0, "%H:%M:%S") m <- matrix(NA, length(x1), length(vline)) dimnames(m) <- list(x1, as.character(vline)) for(i in tail(vline, 7)-d[1]+1){ postx <- strptime(paste(d[i], x1), "%Y-%m-%d %H:%M:%S", tz="JST") for(j in which(ymd_diff==i)){ lag <- difftime(res[[j]], postx, units="secs") flag <- format(postx[which(0<=lag & lag<10)], "%H:%M:%S") m[flag, as.character(d[i])] <- sum(m[flag, as.character(d[i])], 1, na.rm=TRUE) print(c(flag, m[flag, as.character(d[i])])) } }
yl <- range(m, na.rm=TRUE) sp <- mapply(function(x) smooth.spline(m[,x], df=50)$y, seq(ncol(m))) plot(x0, sp[,1], type="n", ylim=yl, xlab="22:20 - 23:10", ylab="Tweets/10 seconds #gochiusa") for(i in 6:11){ lines(x0, sp[,i], col=cols[i-5], lwd=2, lty=1) # spline #lines(x0, m[,i], col=cols[i-5], lwd=2, lty=1) # 生データ } abline(v=seq(s01[1], s01[2], by="min"), lty=3) abline(v=seq(s01[1], s01[2], by="5 min"), lty=1, lwd=1.2) legend("topleft", legend=tail(titles, length(cols)), lty=1, col=cols, lwd=3, bg="white")
宝物は君の決定的瞬間
最終話の放送が終わって27日未明に再度データ収集をして、最終版の結果をまた書きます(続くのかよ
