医薬データ解析のためのベイズ統計学

読んだ。

医薬データ解析のためのベイズ統計学

医薬データ解析のためのベイズ統計学

めっちゃ時間がかかってしかも記事にしておくのも時間かかってた。
 
古典的なp値の頻度論的な話も含んでいて、基礎の確認…にはならない。けっこう大変だった。
しかし、収束の判定法や、変数選択としての使い方とか、他のベイズ本ではあまりしっかり書いていないようなことも書かれているような気がした(読了したのが昔過ぎて記憶があやふやである)。

レプリカ交換法

読んだ。
Bayesian estimation of phase response curves. Neural Netw. 2010 Aug;23(6):752-63.
 
Phase response curve (PRC) という神経細胞の発火の記録を推定したいが、周期のズレや発火タイミングの変化などで普通にやったら推定が収束しないらしい。
レプリカ交換法は現時点でRstan では実装されていないので、コードをコンパイルしたものを使いまわして各iteration の最終サンプリング結果から取り出してきてやる方法がある。
StanとRでレプリカ交換MCMC(parallel tempering) を実行する - StatModeling Memorandum
 
曲線の推定は2階差分としてなめらかになるようにしている。prior に z_{j-1}-2z_j + z_{j+1}=(z_{j-1}-z_j)-(z_j-z_{j+1}) とする。

ギブスサンプリング

MikuHatsune2018-05-29

Metropolis-Hastings サンプリングをやったので、ギブスサンプリングをやってみる。
ある変数X=\{x_1,x_2,\dots,x_{i-1},x_i,x_{i+1},\dots,x_m\} について、i 番目を取り除いたX_{-i}
X_i|X_{-i} \sim x_1,x_2,\dots,x_{i-1},~~~x_{i+1},\dots,x_mi 番目が抜けている)
で順次サンプリングして、その値を入れなおしてまたサンプリングする、を1\dots m について行う。
 
結局、あるひとつi 以外のすべてを固定して、ひとつサンプリングすることになる。多次元正規分布の場合は、こちらを参考に
\mu_{i|-i}\leftarrow \mu_{i}+\Sigma_{-i-i}^{-1}(X_{-i}-\mu_{-i})
\Sigma_{i|-i}\leftarrow \Sigma_{ii}-\Sigma_{-i-i}^{-1}\Sigma_{-ii}
とする。ただし\underbrace{\Sigma_{-i-i}^{-1}}_{1\times (m-1)}=\underbrace{\Sigma-ii}_{1\times(m-1)}\underbrace{\Sigma_{-i-i}}_{(m-1)\times(m-1)} である。なので\mu_{i|-i}\Sigma_{i|-i} はそれぞれ1\times 1 の大きさというかスカラーで、各X_i はひとつずつサンプリングされる。
スクリプトここからパクる。
 
初期値が真の分布からは遠いところにしておいて、収束するまでにちょっと時間がかかるようにする。
Metropolis-Hastings(MH)では、真の分布に近づくまでに少し時間がかかる。サンプリングされた分布は、xy 同時にサンプリングしているので、移動方向は斜めである。
Gibbs (GB)は、各X_i でサンプリングしているので、動き方は各x, y 方向にジグザグである。詳細釣り合い条件や移動確率が1であることはここでは取り上げないが、動きやすいのでMH より収束に向かいやすい。
この2次元の例では、一発目から真の分布に取り込まれている様子がわかる。


 
多次元に拡張するが、図示することを考えて3次元にしてみる。MH もGB もひだり下のマゼンタのところからサンプリングを開始したが、HM は最初の50点の時点ではまだ収束していないが、GB は7点目くらいからもう収束している。

MH の移動方向は斜めだが、GB では各点がサンプリングされるまでに、xyz でそれぞれX_i をサンプリングしているので、カクカクになっている。

# 2D
# 初期値
x0 <- c(-5, 2)
iter <- 1500 # 繰り返し回数
x <- matrix(0, iter, 2) # ランダムウォークの座標
dratio <- rep(0, iter) # 確率密度比
out <- rep(0, iter) # 確率密度比にしたがって確率的に動いたか、動かなかったかの記録
x[1,] <- x0

# 2次元正規分布の相関
rho <- 0.7
sig <- matrix(c(1, rho, rho, 1), 2)
mu <- c(4, 5) # 2次元正規分布の各パラメータの平均

# Metrololis-Hastings
for(i in 2:iter){
  walk <- runif(2, -0.5, 0.5) # 歩く量
  v <- x[i-1,] + walk # 次の候補点
  dratio[i-1] <- dmvnorm(v, mean=mu, sigma=sig)/dmvnorm(x[i-1,], mean=mu, sigma=sig) # 確率密度比
  out[i-1] <- rbinom(1, 1, min(1, dratio[i-1])) # 動いたか、動かなかったか
  x[i,] <- x[i-1,] + out[i-1]*walk
}

kd1 <- kde2d(x[,1], x[,2], c(bandwidth.nrd(x[,1]), bandwidth.nrd(x[,2])), n=1000)
cols <- jet.colors(iter)
plot(x, type="p", pch=16, cex=0.6, col=2, xlab="x1", ylab="x2", main="Metropolis-Hastings sampling", xlim=xl, ylim=yl)
for(i in 2:iter){
  segments(x[i-1,1], x[i-1,2], x[i,1], x[i,2], col=cols[i])
}
points(x[1,1], x[1,2], pch="★", col=6)
points(mu[1], mu[2], pch="★", col=5)
contour(kd1, add=TRUE, col=1)

# Gibbs
X <- matrix(x0, nr=1)
for(j in 2:iter){
  x <- X[j-1,]
  for(i in seq(mu)){
    s <- sig[-i, i] %*% solve(sig[-i, -i])   # Σ_ab Σ_bb ^ -1
    # (PRML 2.81) μ_a|b = μ_a + Σ_ab Σ_bb ^ -1 (x_b - μ_b)
    mu_a_b <- mu[i] + s %*% (x[-i] - mu[-i])
    # (PRML 2.82) Σ_a|b = Σ_aa - Σ_ab Σ_bb ^ -1 Σ_ba
    sig_a_b <- sig[i, i] - s %*% sig[i, -i]
   # [Gibbs] x_a 〜 p(x_a|x_{-a}) = N(μ_a|b, Σ_a|b)
    x[i] <- rnorm(1, mu_a_b, sqrt(sig_a_b))
  }
  X <- rbind(X, x)
}
colnames(X) <- sprintf("V%d", seq(ncol(X)))

kd2 <- kde2d(X[,1], X[,2], c(bandwidth.nrd(X[,1]), bandwidth.nrd(X[,2])), n=1000)
cols <- jet.colors(iter)
plot(X, type="p", pch=16, cex=0.6, col=2, xlab="x1", ylab="x2", main="Gibbs sampling", xlim=xl, ylim=yl)
#lines(x, col=cols)
for(i in 2:iter){
  for(j in seq(ncol(X))){
    y <- rbind(replace(X[i-1,], 0:(j-1), X[i, 0:(j-1)]), replace(X[i-1,], 1:j, X[i, 1:j]))
    segments(y[1,1], y[1,2], y[2,1], y[2,2], col=cols[i])
  }
}
points(X[1,1], X[1,2], pch="★", col=6)
points(mu[1], mu[2], pch="★", col=5)
contour(kd2, add=TRUE, col=1)
# 3D
# 初期値
x0 <- c(-5, 2, -4)

iter <- 1500 # 繰り返し回数
x <- matrix(0, iter, length(x0)) # ランダムウォークの座標
dratio <- rep(0, iter) # 確率密度比
out <- rep(0, iter) # 確率密度比にしたがって確率的に動いたか、動かなかったかの記録
x[1,] <- x0
# 3次元の多次元正規分布を適当に作る
rho <- c(0.7, 0.6, 0.9)
sig <- diag(0, length(x0))
sig[lower.tri(sig)] <- rho
sig <- sig + t(sig)
diag(sig) <- rep(1, length(x0))
mu <- c(4, 5, 3)

# Metropolis-Hastings
for(i in 2:iter){
  walk <- runif(length(x0), -0.5, 0.5) # 歩く量
  v <- x[i-1,] + walk # 次の候補点
  dratio[i-1] <- dmvnorm(v, mean=mu, sigma=sig)/dmvnorm(x[i-1,], mean=mu, sigma=sig) # 確率密度比
  out[i-1] <- rbinom(1, 1, min(1, dratio[i-1])) # 動いたか、動かなかったか
  x[i,] <- x[i-1,] + out[i-1]*walk
}
colnames(x) <- sprintf("V%d", seq(ncol(x)))

cols <- jet.colors(iter)
plot3d(x, type="p", pch=16, size=0.01)
points3d(head(x, 50), size=5)
title3d("Metropolis-Hastings sampling", line=3)
for(i in 2:iter){
  segments3d(x[(i-1):i,], col=cols[i])
}
spheres3d(x[1,], col=6, radius=0.5)

# Gibbs
X <- matrix(x0, nr=1)
for(j in 2:iter){
  x <- X[j-1,]
  for(i in seq(mu)){
    s <- sig[-i, i] %*% solve(sig[-i, -i])   # Σ_ab Σ_bb ^ -1
    # (PRML 2.81) μ_a|b = μ_a + Σ_ab Σ_bb ^ -1 (x_b - μ_b)
    mu_a_b <- mu[i] + s %*% (x[-i] - mu[-i])
    # (PRML 2.82) Σ_a|b = Σ_aa - Σ_ab Σ_bb ^ -1 Σ_ba
    sig_a_b <- sig[i, i] - s %*% sig[i, -i]
   # [Gibbs] x_a 〜 p(x_a|x_{-a}) = N(μ_a|b, Σ_a|b)
    x[i] <- rnorm(1, mu_a_b, sqrt(sig_a_b))
  }
  X <- rbind(X, x)
}
colnames(X) <- sprintf("V%d", seq(ncol(X)))
# ギザギザの描画
Y <- X[1, , drop=FALSE]
for(i in 2:nrow(X)){
  for(j in 1:ncol(X)){
    y <- rbind(replace(X[i-1,], 0:(j-1), X[i, 0:(j-1)]), replace(X[i-1,], 1:j, X[i, 1:j]))
    Y <- rbind(Y, y)
  }
}

cols <- jet.colors(iter)
plot3d(X, size=0.01)
points3d(head(X, 50), size=5)
title3d("Gibbs sampling", line=3)
lines3d(Y, col=rep(cols, each=3))
spheres3d(x[1,], col=6, radius=0.5)

時系列分析と状態空間モデルの基礎: RとStanで学ぶ理論と実装

読んだ。

時系列分析と状態空間モデルの基礎: RとStanで学ぶ理論と実装

時系列分析と状態空間モデルの基礎: RとStanで学ぶ理論と実装

COI:自費で買った。
 
時系列分析とRstan を使った状態空間モデルの勉強がこれ一冊でできる、というのは誇大広告で、いわゆる古典的な時系列解析が前半2/3、ベイズ的な状態空間モデルが後半1/3に書いてある。
最近のはやりで、ベイズ的な状態空間モデルの勉強をしよう、と思ってこの本を手に取ると、少しがっかりする。というのも、状態空間モデルについての話はけっこうスカスカだし、読むなら断然こちらである。 
しかし、この本の素晴らしいところは、前半2/3 で古典的な時系列解析についてわかりやすく説明されており、その流れを踏まえて後半1/3 の状態空間モデルを読むと、状態空間モデルそのものについてはもうほとんどわかった状態で読める、ということである。
前半の時系列解析は、自分でも以前勉強はしたがちょっとあやふやで記憶も薄れていたこともあったが、非常にわかりやすかった、と思う。
 
単位根過程は、原系列が非定常過程であり、差分系列が定常過程であるときにいう。
この検定にはKPSS検定とADF検定があるが、両者の帰無仮説の設定が逆で興味深かった。
KPSS検定は、y_t=\alpha+\beta t + \sum_{i=1}^t u_t + \epsilon_t, u_i\sim iid(0, \sigma_u^2) を仮定して
H_0\sigma_u^2=0帰無仮説は単位根なし)
H_1\sigma_u^2\not{=}0ランダムウォークがあればトレンドを除去しても単位根が残るので、対立仮説は単位根あり)
 
ADF検定は、AR(1) モデルであるy_t=\phi_1y_{t-1}+\epsilon_t, \epsilon_t\sim N(0, \sigma^2) について、
H_0\phi_1=1(ホワイトノイズの累積和になり、ランダムウォークとなりつまり帰無仮説として単位根あり)
H_1|\phi_1|<1(対立仮説は単位根なし)
となる。
 
ふたつ検定があれば単位根のあるなしの判定はほぼ確実にできるのだろうか?

FKが一番上手いのは誰? 「過去5年間の成功率」で見る“ベスト・フリーキッカー”TOP20

MikuHatsune2018-02-13

という記事を見かけた。
FKが一番上手いのは誰? 「過去5年間の成功率」で見る“ベスト・フリーキッカー”TOP20

現在、最もフリーキックが上手い選手は誰なのか……。イギリス誌『FourFourTwo』が、「2013−14シーズンからの5年間で最もフリーキックの成功率が高かった選手TOP20」という記事を掲載した。
FKでのゴール数はよく伝えられるが、その成功率は影に隠れがちだ。
弱小チームの場合、ビッグクラブと比べて好位置でのFK自体が少なく、本当にすごいキッカーが埋もれてしまうこともある。また、優れたキッカーの多いチームに在籍している場合も、蹴るチャンスは少なくなってしまう。
一方、チームで常にキッカーを任されている選手は、蹴る機会が多い分、成功率を維持するのが難しくなってくる。
そういった条件も影響し、「過去5年間のFK成功率TOP20」には意外な選手の名前も登場した。
なお、データは1月23日時点のもので、5年間で最低20本以上のFK蹴った選手だけが対象となっている。また、『FourFourTwo』は意外にもランク外となった“名キッカー”も紹介している。

ということなので、28人分のフリーキック回数とゴール数が載っていたので、これをベイズ的に解析してフリーキック得点率を推定する。
単純に二項分布に当てはめるだけ。

特に何も考えない場合の成功率はこんな感じになる。色のついた部分は80% 信用区間であり、ベイズなので「このフリーキックで得点になる確率」と思ってもらってよい。
クリスティアーノロナウド(p[2]) やメッシ(p[3]) は蹴った回数が多いため、推定の精度があがり、幅は狭くなる。
ミラレム・ピアニッチ(p[17]) は76回蹴って11ゴールなので、他の選手と比べて少し幅が狭い。

 
さて、何も考えずにフリーキックを蹴って、入る確率が0% から100% まで一様分布、と期待するのは、サッカーをしていて非常に不自然である。
というわけで、得点確率p を適当に0% (ゴールが決まりにくい)に偏ったベータ分布からサンプリングすることにしよう。
いま、beta(1, 1) ならば一様分布のピンクだが、beta(1, 10) とすると得点確率0% から低いところに非常にかたよる。

 
この条件で同じようにすると、成功確率は過小評価気味になる。

二項分布はベータ関数x^{\alpha-1}(1-x)^{\beta-1} で表現することができて、これは

qbeta(alpha/2, dat$g+1, dat$n-dat$g+1)

で計算できる。ベータ関数のshape パラメータは -1 なので計算するときには +1 を忘れないこと。
このとき、事前分布が一様ならば、何も乗算しなくてよいが、事前分布をベータ分布で偏らせるならば、事前分布はB(1, 10)=x^{1-1}(1-x)^{10-1}=x^0(1-x)^9 であり、これを尤度関数(上のR コードのところ、つまりx^{\alpha-1}(1-x)^{\beta-1})にかけるだけで事後分布になる。つまり、\alpha\beta のところに各々足す。
 
これ読んでおけばいいと思う(唐突な宣伝、COI なし)。

 
期待としては、クリスティアーノロナウドやメッシといった5-10% くらいがせいぜい達成できるFK成功率だと思っていたが、頑張れば20% もいけるのだな、という印象。

library(rstan)
txt <- "name n g
ズラタン・イブラヒモヴィッチ 91 4
クリスティアーノ・ロナウド 112 7
リオネル・メッシ 165 11
ハカン・チャルハノール 99 10
ケヴィン・デ・ブライネ 35 4
アレクシス・サンチェス 42 5
ネイマール 41 5
ウィリアン 31 4
スソ 22 3
ガエル・カクタ 22 3
ウィサム・ベン・イェデル 22 3
ウェイン・ルーニー 29 4
ロメル・ルカク 21 3
ヤヤ・トゥーレ 28 4
ズラトコ・ユヌゾヴィッチ 35 5
ダヴィド・アラバ 42 6
ミラレム・ピアニッチ 76 11
マルクス・ズットナー 27 4
マウリシオ・レモス 20 3
アンリ・セヴェ 20 3
セヤド・サリホヴィッチ 20 3
マルヴィン・プラッテンハルト 37 6
フェデリコ・ヴィヴィアーニ 30 5
ハメス・ロドリゲス 36 6
マルコス・アロンソ 22 4
フィリペ・コウチーニョ 26 5
パウロ・ディバラ 38 8
フアン・マタ 22 5
"


dat <- read.table(text=txt, header=TRUE)
code <- "
data{
  int N_kicker;
  int N[N_kicker];
  int G[N_kicker];
  vector<lower=0>[2] s;
}
parameters{
  vector<lower=0, upper=1>[N_kicker] p;
}
model{
  p ~ beta(s[1], s[2]);
  G ~ binomial(N, p);
}
"

m <- stan_model(model_code=code)

cols <- c("hotpink", "skyblue")
input <- list(N=dat$n, G=dat$g, N_kicker=nrow(dat), s=c(1, 1))
fit <- sampling(m, data=input, seed=1234)

g <- stan_plot(fit, show_density=T, pars="p", outer_level=0.95, ci_level=0.8, fill_color=cols[1])
g + xlim(0, 0.4)

x <- seq(0, 1, length=1000)
d <- cbind(dbeta(x, 1, 1), dbeta(x, 1, 10))
matplot(x, d, type="l", lwd=4, lty=1, col=cols)

# サンプリングからの分布を考えるならば
ex <- extract(fit)
alpha <- 0.05
ci <- apply(ex$p, 2, quantile, c(alpha/2, 1-alpha/2))

# 一様な事前分布
qbeta(alpha/2, dat$g+1, dat$n-dat$g+1)
# ベータ分布
qbeta(alpha/2, dat$g+1+0, dat$n-dat$g+1+9)

(サッカー解説)2点差は危険なスコアですね ← ???

MikuHatsune2018-01-17

高校サッカーを見ていた。2017年度は前橋育英が初優勝で幕を閉じた。
どの試合だったか忘れてしまったが、2点差がついたときに解説が「2点差は危険」ということを言っていた。
 
調べてみると、やはりよく言われていることのようだが、実際にデータをとってみると、プレミアリーグでは2点差をひっくり返して勝つ確率は1.71%、Jリーグでは2点差からドローが5%、2点差から敗北5%だったらしい。
 
自分はユース年代のファンなので、せっかくJFA が公式に試合記録を出してくれるということもあって、冬の高校選手権の得点時間を抽出して、2点差が危険なのかどうかを解析したい。
 
JFA から公式記録PDF を取得するが、2009年(88回大会)から2017年(96回大会)まで存在していて、各大会47試合ある。ただし、2009年はPDF の都合でデータをパースできなかったので全部で379試合が対象である。
試合記録から、得点時間と得点チームを抽出する。というのも、得点時間の分布と、ある時刻でのリードが最終的に勝率にどのくらい影響するかをみたいからである。
ただし、冬の選手権は決勝までは80分+即PK戦、決勝は90分+15分ハーフの延長戦なので、いい感じに分類したかったが、90分のときのロスタイムでうまく処理できてない気もするが無視。
また、2点差が危険の定義だが、N点差がついていたにもかかわらず、追いつかれて点差が0 になったときとする。というのも、追いつかれたとしてもトーナメント方式なのでPK戦までして勝敗を決めるが、得点時間を取得した都合で、80分内での点差の推移を考えることにする。なので、80分時点での点差が勝ち、引き分け、負けになる。
90分の試合を補正して80分にしてもよかったが、していない。ロスタイムの得点で最大80+6分というのがある。
 
結論から言うと、2点差がつくとたいてい90% 以上の確率でリードを保ったまま試合を終えることができる。
ある試合時刻t のとき、2点差がついている試合数N_{t,2} があり、残りの80-t 時間、リードを保ったまま(点差が0 より大)試合を終えた試合数N_{t,2}^{'} を単純に数えて割ったものが各点であり、指数関数モデルで曲線フィッティングしたのが黒線である。
試合経過時間が短いときに得点しても、追いつくチャンスがまだまだあるので、リードを保ち続けられる確率は低めだが、時間が経てばそのチャンスが減っていくので、リードを保てる確率は増える。
1点差のリードのとき、80-85分のあたりで確率が低くなっており、このせいで推定曲線が全体的に下がったものが推定されているが、ここを除外して推定曲線を引くと、80分時点で1点リードしていればほぼ勝率が100% になる。

 
前半、後半について得点時間の分布を見てみると、特にどの時間帯でピークがある、という感じではなく、一様分布のようである。なのでどの時間帯にトイレにいっても得点シーンを見逃す確率は同じ程度と思われる。

 
一方で、得点が入ってから次の得点がはいるまでの時間の分布は、一様分布ではない。横軸が経過時間、縦軸が確率密度とすれば、ガンマ分布のようである。また、1試合あたりの得点数も、ポアソン分布っぽい。
というわけで、次の得点までの時間分布をガンマ分布で、1試合あたりの両チーム合計の得点数をポアソン分布でモデル化してみると、こんな感じになる。
ガンマ分布のパラメータは(1.33, 0.06) で、最頻5分、中央16分程度、平均20分で点がはいる。
ポアソン分布のパラメータは2.9であり、1試合あたり3点くらいの点がはいる。

 
というわけで、「2点差は危険」というのは、1回のゴールで1点しか点が動かないサッカーのルールの特性上、2点差から追いつくもしくは逆転されるという試合が伝説のように語られる思い出バイアスなので、解説の人は勉強してください()
 
解析

### analysis
library(stringr)
dat <- read.table("soccer.txt", header=TRUE, stringsAsFactors=FALSE)

tmp <- split(dat[,-1], dat$year)
tmp <- mapply(function(z) split(z[,-1], z$no), tmp)

res <- vector("list", length(tmp))
k <- c(1, -1)
for(i in seq(tmp)){
  out <- vector("list", length(tmp[[i]]))
  for(j in seq(tmp[[i]])){
    hoge <- tmp[[i]][[j]]
    foo <- lapply(str_extract_all(hoge$time, "\\d+"), as.numeric)
    time <- mapply(function(z) append(z, rep(0, 2-length(z))), foo)
    team <- k[c(factor(hoge$team, levels=unique(hoge$team)))]
    bar <- rbind(time, team)
    rownames(bar) <- c("time", "extra", "team")
    out[[j]] <- bar
  }
  res[[i]] <- out
}


y0 <- lapply(res, do.call, what=cbind)
y1 <- do.call(cbind, y0)
half <- c("前半", "後半")
idx <- list(y1[1,] <= 40, y1[1,] >  40)
goaltime <- mapply(function(z) colSums(y1[1:2, z]), idx)
yl <- c(0, 4)
par(mfcol=c(2, 2), las=1)
for(i in seq(goaltime)){
  pl1 <- table(goaltime[[i]])/sum(unlist(idx), na.rm=TRUE)
  pl2 <- table(goaltime[[i]])/length(goaltime[[i]])
  barplot(pl1*100, xlab="経過時間", ylab="percentage", ylim=yl, col=i+1, main=sprintf("%sの得点時間分布 (%s時間)", half[i], "全試合"))
  barplot(pl2*100, xlab="経過時間", ylab="percentage", ylim=yl, col=i+1, main=sprintf("%sの得点時間分布 (%s時間)", half[i], half[i]))
}


score <- vector("list", sum(sapply(tmp, length)))
k <- 1
for(i in seq(tmp)){
  for(j in seq(tmp[[i]])){
    hoge <- res[[i]][[j]]
    score[[k]] <- cbind(c(0, colSums(hoge[1:2,,drop=FALSE])), cumsum(c(0, hoge[3,,drop=FALSE])))
    k <- k + 1
  }
}

mat <- NULL
for(i in seq(res)){
  for(j in seq(res[[i]])){
    k <- res[[i]][[j]]
    if(!is.na(k[1])){
      l <- append(rep(NA, max(k[1,])), rep(NA, tail(k[2,], 1)+1))
      elp <- c(1, colSums(k[1:2,,drop=FALSE])+1) # 得点時間のindex
      elp <- append(elp, ifelse(max(k[1,])<=80, max(80, elp), ifelse(max(k[1,])<=90, max(90, elp), 120)))
      s <- c(0, cumsum(k[3,]))
      for(m in 1:(length(elp)-1)){
        l[ elp[m]:elp[m+1] ] <- s[m]
      }
    } else {
      l[1:80] <- 0
      l <- rep(0, 80)
    }
    l <- append(l, rep(NA, 120-length(l)))
    mat <- rbind(mat, l)
  }
}
library(rstan)
# 同時に
code <- "
data{
  int Ngame;
  int M;
  int Ngoal[Ngame];
  vector<lower=-1>[M] G[Ngame];
}
parameters{
  real<lower=0, upper=10> lambda;
  real<lower=0, upper=5> p[2];
}
model{
  Ngoal ~ poisson(lambda);
  for(i in 1:Ngame){
    G[i, 1:Ngoal[i] ] ~ gamma(p[1], p[2]);
  }
}
"

m <- stan_model(model_code=code)
Ng <- sapply(goaltime, length)
game <- t(mapply(function(z) append(z, rep(-1, max(Ng)-length(z))), goaltime))
dat <- list(Ngame=length(Ng), Ngoal=Ng, M=max(Ng), G=game)
fit <- sampling(m, data=dat, iter=3000, warmup=1500, seed=1234)
ex <- extract(fit)


par(mfrow=c(1, 2), las=1, mar=c(5, 4.5, 2, 2), cex.lab=1.5)
goaltime <- mapply(function(z) head(rle(as.numeric(na.omit(mat[z,])))$length, -1), seq(nrow(mat)))
g <- unlist(goaltime)
x <- 1:100
tab <- table(factor(g, x))

dg <- dgamma(x, median(ex$p[,1]), median(ex$p[,2]))
plot(tab/length(g)*100, xlab="直前のゴールからの経過時間(分)", ylab="頻度 [%]")
lines(x, dg*100, col=2, lwd=3)

x <- 0:max(Ng)
Ngd <- table(factor(Ng, x))
dp <- dpois(x, median(ex$lambda))
plot(Ngd/length(Ng)*100, xlab="1試合の総得点数", ylab="頻度 [%]")
lines(x, dp*100, col=2, lwd=3)


# ある時刻での勝利確率を求める
mat80 <- mat[is.na(mat[,90]),]
s <- replicate(2, matrix(0, 2, 90), simplify=FALSE)
for(k in 1:2){
  for(i in 1:nrow(mat80)){
    a <- abs(mat80[i,])
    M <- sum(!is.na(a))
    for(j in 1:M){
      if(a[j] >= k){
        s[[k]][2, j] <- s[[k]][2, j] + 1
        if(all(a[j:M] > 0)){
          s[[k]][1, j] <- s[[k]][1, j] + 1
        }
      }
    }
  }
}

xl <- c(0, 90)
yl <- c(0.5, 1)
cols <- c("blue", "red")
plot(0, type="n", xlim=xl, ylim=yl, xlab="試合経過時間(分)", ylab="試合終了時までリードし続ける確率", cex.lab=1.2, cex.axis=1.2, las=1)
for(i in seq(s)){
  v <- s[[i]][1,]/s[[i]][2,]
  points(v, col=cols[i], pch=16)
  d <- drm(v ~ seq(v), fct=AR.3())
  lines(d$data[,1], d$predres[,1], lwd=3)
}
legend("bottomright", legend=sprintf("%d 点差", 1:2), col=cols, pch=16, cex=2)

 
データ取得と前処理

# 2017-2015 年までの96-93 大会
wd <- "/soccer/"
dir.create(wd)
y <- 2017:2014
m <- 1:47

for(i in y){
  urls <- sprintf("http://www.jfa.jp/match/alljapan_highschool_%d/match_report/m%d.pdf", i, m)
  dir.create(sprintf("%s%d", wd, i))
  for(j in urls){
    command <- sprintf("wget %s -P %s%d", j, wd, i)
    system(command)
  }
}

for(i in 2013){
  urls <- sprintf("https://www.jfa.or.jp/match/matches/2014/0113koukou_soccer/match_report/m%d.pdf", m)
  dir.create(sprintf("%s%d", wd, i))
  for(j in urls){
    command <- sprintf("wget %s -P %s%d", j, wd, i)
    system(command)
  }
}

y <- 2012:2008
d1 <- c("0114", "0110", "0110", "0110", "")
for(i in seq(y)){
  urls <- sprintf("http://www.jfa.or.jp/match/matches/%d/%skoukou_soccer/schedule_result/pdf/m%02d.pdf", y[i]+1, d1[i], m)
  dir.create(sprintf("%s%d", wd, y[i]))
  for(j in urls){
    command <- sprintf("wget %s -P %s%d", j, wd, y[i])
    system(command)
  }
}
system(sprintf("mv %s%d/* %s%d/; rm -r %s%d", wd, 2008, wd, 2009, wd, 2008))

###
library(pdftools)
library(stringr)
year <- list.files(wd)
year <- as.character(2017:2013)
# 2013-2017
output <- NULL
for(j in seq(year)){
  print(sprintf("%s", year[j]))
  files <- list.files(sprintf("%s%s", wd, year[j]), pattern="pdf", recursive=TRUE, full.names=TRUE)
  for(f in seq(files)){
  text <- pdf_text(files[f])
  text <- strsplit(text, "\n")[[1]]
  i <- 1
  while(i < length(text)){
    if(any(strsplit(text[i], " ")[[1]] == "得点時間")){
      i <- i + 1
      res <- NULL
      while(length(grep("PK戦の経過", text[i])) < 1){
        tmp <- strsplit(gsub(" +", " ", text[i]), " ")[[1]][2:3] if(length(grep("分", tmp)) > 0){
          res <- rbind(res, tmp)
        }
        i <- i + 1
      }
      tmp <- table(res[,2])
      if(length(tmp) > 1){
        if(tmp[1] == tmp[2]){
          i <- i + 1
          a1 <- strsplit(gsub(" +", " ", text[i]), " ")[[1]]
          a1 <- a1[a1 != " "]
          i <- i + 1
          a2 <- strsplit(gsub(" +", " ", text[i]), " ")[[1]]
          a2 <- a2[a2 != " "]
          a3 <- c(a1[1], a2[1])[which.max(c(length(grep("○", a1)), length(grep("○", a2))))]
          res <- cbind(res, (res[,2] == a3)+0)
        } else {
          i <- i + 100
          a3 <- names(tmp)[which.max(tmp)]
          res <- cbind(res, (res[,2] == a3)+0)
        }
      } else {
        res <- cbind(res, 1)
      }
    }
    i <- i + 1
  }
  res <- switch(1 + (length(res)==1), res, matrix(c(NA, NA, res), 1))
  res <- cbind(year[j], f, res)
  output <- rbind(output, res)
  }
}

year <- as.character(2012:2009)
for(j in seq(year)){
  files <- list.files(sprintf("%s%s", wd, year[j]), pattern="pdf", recursive=TRUE, full.names=TRUE)
  for(f in seq(files)){
  text <- pdf_text(files[f])
  text <- strsplit(text, "\n")[[1]]
  if(length(text) > 0){
    i <- 1
    while(i < length(text)){
      pktext <- str_extract(text[i], "先.*先")
      if(!is.na(pktext)){
        pk <- as.numeric(str_extract_all(pktext, "\\d+")[[1]])
      }
      if(any(strsplit(text[i], " ")[[1]] == "得点チーム")){
        i <- i + 1
        res <- NULL
        while(length(grep("備考欄|戦評者氏名", text[i])) < 1){
          tmp <- strsplit(gsub(" +", " ", text[i]), " ")[[1]]
          tmp <- tmp[grep("分", tmp)]
          if(length(grep("分", tmp)) > 0){
            tmp <- strsplit(gsub("分", "分 ", tmp), " ")[[1]]
            res <- rbind(res, tmp)
          }
          i <- i + 1
        }
        tmp <- table(res[,2])
        if(length(tmp) > 1){
          if(tmp[1] == tmp[2]){
            i <- i + 1
            a1 <- strsplit(gsub(" +", " ", text[i]), " ")[[1]]
            a1 <- a1[a1 != " "]
            i <- i + 1
            a2 <- strsplit(gsub(" +", " ", text[i]), " ")[[1]]
            a2 <- a2[a2 != " "]
            a3 <- c(a1[1], a2[1])[which.max(c(length(grep("○", a1)), length(grep("○", a2))))]
            res <- cbind(res, (res[,2] == a3)+0)
          } else {
            i <- i + 100
            a3 <- names(tmp)[which.max(tmp)]
            res <- cbind(res, (res[,2] == a3)+0)
          }
        } else {
          res <- cbind(res, 1)
        }
      }
      i <- i + 1
    }
    res <- switch(1 + (length(res)==1), res, matrix(c(NA, NA, res), 1))
    res <- cbind(year[j], f, res)
    output <- rbind(output, res)
    } else {
      print(files[f])
    }
  }
}

きららフェスタ2017 に出演していた声優たちの集客力をRstan で推定する

MikuHatsune2017-12-03

この記事は
RStudio Advent Calendar 2017 - Qiita
まんがタイムきらら Advent Calendar 2017
ごちうさ Advent Calendar 2017
Stan Advent Calendar 2017 - Qiita
R Advent Calendar 2017 - Qiita
の3日目の配当記事です。
 
声優統計第9号で、きららフェスタ2016に出演した声優の集客力について検討したが、2017年度もキララフェスタに参加したのでまったく同じことをやってみる。
 
2017年では4作品のアニメが出て、東山奈央大久保瑠美が総合司会ということで出演していた。この二人は4作品には出演していなかった。
各声優のアニメ出演パターンは以下の通りである。茅野愛衣が3作品に出ている。水瀬いのりも意外にサクラクエストに出演していた。

 
各声優の集客力の推定結果は以下の通りである。サクラクエスト声優陣の推定がグダグダである。
2016年度ほどに、いろいろなアニメにほどよく散らばって出演している声優が少なくて推定結果はグダグダになったようである。

 
さで、Rstudio アドベントカレンダーでもあるので、Rstan の推定結果をRstudio からRpubs へ出すときの注意点を述べる。Rstudio でknitr すると、library の読み込み時に関連メッセージがでてきて、これがhtml に紛れ込んでしまう。
library のオプションでquietly=FALSE にしても、チャンクオプションでinclude=FALSE にしても出てきてしまうのでイラッ☆とする。
ここで、library で読み込まずに、namespace::function として呼び出すとこの関連メッセージを呼び起こすことなく関数を実行できるっぽいので、こうする。

さて、rstan の実行結果を使ってhtml に計算結果や図を出したいが、rstan のsampling はいちいち実行メッセージを出してしまうので、これもすべてhtml に紛れ込んでしまう。なので、sampling した結果をsave 関数によって保存して、knitr するときにはload することによってknitr でコンパイルするときにはrstan による演算はしないことで回避できるっぽい。

ということで、Rstudio でknitr して、Rpubs にpublish ボタンでアップロードすれば、html やrevealjs がいろいろ作成できる。
JapanR#5 お疲れ様でした。こんな感じでrevealjs のスライドができます。
RPubs - 20171202japanr#5

anime,bd,安済知佳,上田麗奈,大久保瑠美,茅野愛衣,久保ユリカ,小松未可子,高田憂希,竹尾歩美,田中ちえ美,東山奈央,徳井青空,戸田めぐみ,七瀬彩夏,原田彩楓,日笠陽子,本渡楓,水瀬いのり,村川梨衣,山口愛,佳村はるか
ご注文はうさぎですか,11038,0,0,0,1,0,0,0,0,0,0,1,0,0,0,0,0,1,1,0,0
NEWGAME,6452,0,0,0,1,0,0,1,1,0,0,0,1,0,0,1,0,0,0,1,0
サクラクエスト,1165,1,1,0,0,0,1,0,0,1,0,0,0,1,0,1,0,1,0,0,0
うらら迷路貼,694,0,0,0,1,1,0,0,0,0,0,0,0,0,1,0,1,0,0,0,1