廿TT

譬如水怙牛過窓櫺 頭角四蹄都過了 因甚麼尾巴過不得

Albert (2008) 打者の調子の波のモデル化(後編)

Albert (2008) 打者の調子の波のモデル化(前編) - 廿TT の続きです。

以降の分析の目的は、よく言われる野球選手の「調子の波」を選手間で比較可能な指標にするにはどうするか、ということです。

分析対象のデータはカルロス・ギーエンという選手の2005年の打撃成績のデータで、ヒットを 1、アウトを 0 とコード化してある。

GuillenC <- c(0, 1, 0, 1, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0,
1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 1, 0, 0, 1, 1, 0,
1, 0, 1, 1, 0, 1, 0, 1, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 1, 0, 1, 0, 0, 1, 1, 0, 0,
0, 1, 0, 1, 0, 0, 0, 1, 1, 1, 0, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 0, 0, 1, 0, 0, 1, 0, 1,
0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 1, 0,
0, 0, 0, 1, 1, 1, 1, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 1, 0,
0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0,
0, 1, 0, 0, 0, 0, 1, 1, 0, 1, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0,
1, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0,
0, 1, 1, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 1, 1, 1, 0, 1, 0, 0,
0, 0, 0, 1, 0, 1, 1, 0, 0, 0, 1, 0, 0, 0, 1)

この 0, 1 を表のように集計する。

20 はおよそ 4 試合ごとの打席数とのこと。

ヒットの数 打席数
5 20
5 20
7 20
10 20
10 20
10 20
6 20
9 20
4 20
4 20
6 20
7 20
4 20
2 20
6 20
12 34

ヒットの数は二項分布に従うと仮定する。

さらに二項分布の成功確率 p はベータ分布に従うとする。

ここでの工夫はベータ分布を以下のようにパラメタライズすること。

\displaystyle \frac{1}{B(K\eta,K(1-\eta))}p^{K\eta-1}(1-p)^{K(1-\eta)-1}

こうすると η は打率に対応するパラメータ、K は打率の精度に関するパラメータと解釈できる。

K が大きいほどばらつきが小さくなる。

Albert (2008) では ηK に対して以下の無情報的事前分布を仮定している。

\displaystyle g(\eta, K) \propto \frac{1}{\eta (1-\eta)}\frac{1}{(1+K)^2}

ちなみにこの事前分布を外して一様分布を仮定すると、MCMC がまったく収束しなかった。

下図のような具合。

f:id:abrahamcow:20170413193142p:plain

事前分布を入れるとこんな感じで収束しました。

f:id:abrahamcow:20170413193211p:plain

推定には Stan を使った。

Stan のコードはこう。

data {
  int<lower=0> m;
  int<lower=0> n[m];
  int<lower=0> x[m];
  real<lower=0,upper=1> avg;
}
parameters {
  real<lower=0, upper=1> p[m];
  real<lower=0> K;
  real<lower=0, upper=1> eta;
}
model {
  x ~ binomial(n, p);
  p ~ beta(K*eta, K*(1-eta));
  target += -(log(eta)+log(1-eta))-2*log(1+K);
}
generated quantities{
  int<lower=0> pred1[m];
  int<lower=0> pred2[m];
  for(i in 1:m){
    pred1[i] = binomial_rng(n[i], p[i]);
    pred2[i] = binomial_rng(n[i], avg);
  }
}

R のコードはこう。

library(rstan)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
len <-length(GuillenC)
m <- len %/% 20
x <- numeric(m)
j<-0
for(i in 1:(m-1)){
  x[i] <- sum(GuillenC[seq(j+1,j+20,by=1)])
  j <- j+20
}
x[m] <- sum(GuillenC[(j+1):len])
n=c(rep(20,m-1),len-j)
avg=mean(GuillenC)
dat4stan <-list(m=m,x=x,n=n,avg=avg)
fitbetabinom <- stan("betabinom.stan",data = dat4stan,iter=10000,thin=4)

事後期待値はそれぞれ \hat K= 52.98, \hat \eta= 0.32 だった.

20 打席ごとに区切った打率と \hat p の推移を見るとこんな感じである。

phat <-get_posterior_mean(fitbetabinom,pars="p")[,"mean-all chains"]
eta <-get_posterior_mean(fitbetabinom,pars="eta")[,"mean-all chains"]
lower <- summary(fitbetabinom,pars="p")$summary[,"2.5%"]
upper <- summary(fitbetabinom,pars="p")$summary[,"97.5%"]
dat <- data.frame(time=1:m,x=x,n=n,phat=phat,lower=lower,upper)
library(cowplot)
ggplot(dat,aes(x=time))+
  geom_pointrange(aes(y=phat,ymin=lower,ymax=upper,colour="phat"))+
  geom_point(aes(y=x/n,colour="avg"),size=2)+
  geom_hline(aes(yintercept=avg,colour="avg"),linetype=2)+
  geom_hline(aes(yintercept=eta,colour="phat"),linetype=2)+
  ylim(c(0,1))+xlab("")+ylab("")+
  theme(legend.title=element_blank())

f:id:abrahamcow:20170413193805p:plain

青い棒は 95% ベイズ信頼区間を表す。

MCMC を回すのと同時に generated quantities ブロックでヒット数の予測分布を生成していた。

ベータ二項モデルの予測分布のヒストグラムに、実際のヒット数からカーネル密度推定した値を重ねたのが次のプロット、

ex_pred1 <-extract(fitbetabinom,"pred1")
ex_pred2 <-extract(fitbetabinom,"pred2")
hist(ex_pred1$pred1,freq=FALSE,main="",xlab="")
lines(density(x),lwd=2,col="blue")

f:id:abrahamcow:20170413194321p:plain

普通の二項分布を使った予測分布のヒストグラムに、実際のヒット数からカーネル密度推定した値を重ねたのが次のプロットである。

hist(ex_pred2$pred2,freq=FALSE,main="",xlab="")
lines(density(x),lwd=2,col="blue")

f:id:abrahamcow:20170413194427p:plain

普通の二項分布を使うよりもベータ二項モデルのほうが、当てはまりがましになっていることが見て取れる。

Rで学ぶベイズ統計学入門

Rで学ぶベイズ統計学入門

メジャーリーグの数理科学〈下〉 (シュプリンガー数学リーディングス)

メジャーリーグの数理科学〈下〉 (シュプリンガー数学リーディングス)