2010年声優言及数 作業メモ(4)データ集計からグラフ出力まで、のコード

JSONデータの取得まではできたので、あとはRでデータを読み込んで集計してグラフにするよ。

とりあえず先にコードを貼っておく。何だかんだで毎年フルスクラッチしている気がする・・・楽しいからいいけど。R書く機会って少ないし。

library(rjson)
library(ggplot2)
library(cairoDevice)

## functions

# linear regression
sy.stats <- function(df, df.range=11:25, df.x=df.range){
  ddply(df, .(keyword), function(x){
    df.coef <- lm(as.numeric(x[df.range]) ~ df.x)$coefficients
    data.frame(a=df.coef[2], b=df.coef[1])
  })
}

# percentile (cf. http://phi.med.gunma-u.ac.jp/swtips/percentile.R)
percentile.table <- function(dat, p=2) {
  pt1  <- quantile(dat, probs=seq(10^(-p), 1, 10^(-p)), type=3) # SAS仕様でパーセンタイルを求める
  unique(as.data.frame(pt1), fromLast=T)            # convert into data.frame to preserve name
}
percentile.conv  <- function(dat, table, p=2) {
  if(missing(table)){
    table <- percentile.table(dat, p)
  }
  pt4  <- as.numeric(strsplit(rownames(table), "%"))
  pt4[as.integer(cut(dat, c(min(dat, table[,1])-1, table[,1]), labels=seq(pt4)))]
}

# create ggplot2 graph
#
# input data format:
#
#    keyword 2009-01-01 2009-02-01 2009-03-01
# 160 早見沙織      5      3      3
#
# NOTE:
# 1. locale should be set as "C"       (e.g. Sys.setlocale("LC_TIME", "C"))
# 2. font family "MSGothic" must be defined  (e.g. windowsFonts(MSGothic=windowsFont("MS Gothic")))
#
plot.sy <- function(d, sy.range, outdir=".")
{
  keyword <- as.character(d$keyword)

  # moving average
  n  <- 3
  ma <- filter(ts(as.integer(d[-1]), freq=1), rep(1/n, n))

  # plot
  x <- as.Date(colnames(d)[sy.range])
  y <- as.integer(d[sy.range])
  q <- qplot(x, y, geom=c("point", "path"), colour="#1E5692", xlab="", ylab="", main=sprintf("言及数:%s", keyword))

  # plot moving average
  s.x <- x
  s.y <- ma[sy.range-1]  # first element is already deleted
  q <- q + stat_smooth(aes(x=s.x[!is.na(s.y)], y=s.y[!is.na(s.y)]), col="#DA2025", fill="#2578C5", alpha=0.1)

  # set X scale
  q <- q + scale_x_date(major="3 months", minor="months", format = "%b%y")

  # misc settings
  q <- q + theme_bw() + structure(list(
                     plot.title  = theme_text(family = "MSGothic", size=9),
                     axis.text.x = theme_text(family = "MSGothic", size=9),
                     axis.text.y = theme_text(family = "MSGothic", size=9)),
                  class="options")
  q <- q + scale_colour_identity()
  q <- q + opts(legend.position = "none")

  # output to PNG file
  ggsave(q, file=sprintf("%s\\%s_%s.png", outdir, rownames(d), keyword),
       device=Cairo_png, pointsize=8, width=6, height=4)
  return(q)
}

## main

## 1. read keyword data (JSON format)
#
# format of sy.data:
#      keyword     from       to         X1 
# 1   いのくちゆか 2009-01-01 2010-12-31 13 
# 2 おみむらまゆこ 2009-01-01 2010-12-31  0 
#
sy.data <- ldply(list.files("data", "\\.json", full=T), function(f){
  d <- fromJSON(readLines(f))[[1]]
  r <- as.integer(sapply(d$refer_count, function(x) ifelse(is.null(x), 0, x)))
  with(d, data.frame(keyword, from, to, t(r)))
})

## 2. merge into month-wise data
#
# format of sy.data.m:
#      keyword 2009-01-01 2009-02-01 
# 1   いのくちゆか     46     13 
# 2 おみむらまゆこ      9     15 
#
sy.data.m <- ddply(sy.data, .(keyword), function(x){
  label  <- cut(seq(as.Date(x$from), as.Date(x$to), by=1), "month")
  tapply(as.integer(x[4:ncol(x)]), label, sum)
})

## 3. convert sy.data.m into percentile
# sy.data.mp has the same format as sy.data.m
sy.data.mp <- data.frame(sy.data.m$keyword, apply(sy.data.m[,-1], 2, percentile.conv, p=3))
colnames(sy.data.mp) <- colnames(sy.data.m)


## 4. calc statistics
sy.range     <- 11:25  ## 1 quarter of last year + this year
sy.stats.m   <- sy.stats(sy.data.m,  sy.range)
sy.stats.mp  <- sy.stats(sy.data.mp, sy.range)
sy.stats.sum <- apply(sy.data.m[sy.range], 1, sum)

sy.stats.all <- data.frame(
               keyword = sy.stats.m$keyword,
               sum     = sy.stats.sum,
               a       = sy.stats.m$a,
               pa      = sy.stats.mp$a,
               a1pa    = sy.stats.m$a^1 * sy.stats.mp$a,
               a2pa    = sy.stats.m$a^2 * sy.stats.mp$a,
               a3pa    = sy.stats.m$a^3 * sy.stats.mp$a
)

# pick up "increasing" keywords
sy.stats.inc <- sy.stats.all[((sy.stats.all$a > 0) & (sy.stats.all$pa > 0)),]
# head(sy.stats.inc[sort.list(sy.stats.inc$a1pa, dec=T),], 30)

## 5. print graphs

# set locale (for abb names of month)
t <- Sys.getlocale("LC_TIME")
Sys.setlocale("LC_TIME", "C")
windowsFonts(MSGothic=windowsFont("MS Gothic"))

# print graphs (retain plot data in case for further use)
plots <- daply(sy.data.m, .(keyword), plot.sy, sy.range, "plot", .progress="text")

# reset locale (e.g. "Japanese_Japan.932")
Sys.setlocale("LC_TIME", t)

あとは個別にコメントというかメモ。

パーセンタイル変換

今年のコンセプトとしては、言及数の増加だけでなく%の増加も考慮したいなと。去年少し考えた順序ベースのヤツとコンセプトは似たような感じになるのかな?

2009年声優言及数 作業メモ 反省会(1) - XXXannex

で。XX分点のデータを作るのは簡単だけど、じゃあ元のデータがどの部分に入ってくるのか?という処理をどうしよう。まさか配列を全部なめる、まあバイナリーサーチでもいいけど、いちいちそういうことをするのか?って思って若干憂鬱になっていたのですが、こちらにズバリなものがありました。

http://phi.med.gunma-u.ac.jp/swtips/percentile.R

中澤先生マジ神。今までピンとこなかったcutのすごさを思い知りました。

パーセンタイルの表を別なデータにも使いまわしたいので、若干変更しました。あと、些細なバグがあったので直しておきました。

percentile.table <- function(dat, p=2) {
  pt1  <- quantile(dat, probs=seq(10^(-p), 1, 10^(-p)), type=3) # SAS仕様でパーセンタイルを求める
  unique(as.data.frame(pt1), fromLast=T)            # convert into data.frame to preserve name
}
percentile.conv  <- function(dat, table, p=2) {
  if(missing(table)){
    table <- percentile.table(dat, p)
  }
  pt4  <- as.numeric(strsplit(rownames(table), "%"))
  pt4[as.integer(cut(dat, c(min(dat, table[,1])-1, table[,1]), labels=seq(pt4)))]
}

移動平均

毎年回帰直線を引いてたけど、考えてみれば別にいらないか?と思えてきたので今年は移動平均

ma <- filter(ts(as.integer(d[-1]), freq=1), rep(1/n, n))

ここ数日勉強していた時系列分析の成果が!・・・1行で終わってしまった。まあいいけど。

ggplot2で日本語が出力できた!けど・・・

ggplot2でアンチエイリアス&アルファブレンド&日本語出力なグラフを描く方法をずっと考えてたんだけど、それなりに実現できそうです。

まずテーマに日本語フォントを指定します。

  q <- q + theme_bw() + structure(list(
                     plot.title  = theme_text(family = "MSGothic", size=9),
                     axis.text.x = theme_text(family = "MSGothic", size=9),
                     axis.text.y = theme_text(family = "MSGothic", size=9)),
                  class="options")

ggsaveにはdeviceオプションで好きなグラフィックデバイスを渡すことができるようです。今回はCairo_png

ggsave(q, file=sprintf("%s\\%s_%s.png", outdir, rownames(d), keyword), device=Cairo_png)

無事日本語の出力ができた!けど、なぜか明朝体になってしまう。ウインドウに出力させるとフォントの設定が効いてるようだし、日本語フォントを設定しないと文字化けするので、フォントの設定かcairoDeviceの問題だと思う。

追いかけるのも面倒なので、出ただけでもよしとします。