読者です 読者をやめる 読者になる 読者になる

廿TT

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

打者の調子の波のモデル化 4

データ

MLB - スポーツナビプロ野球選手ボガーツの打席結果が掲載されています。

例として下表に一部を抜き出します。

日付 打数 安打
6/1 5 1
6/2 5 1
6/3 5 2
6/4 3 0
6/5 4 3
6/6 4 0

このデータをとってくる R のコードは以下のようになります。

# 後で使うパッケージをまとめてロード
library(rvest)
library(dplyr)
library(pipeR)
library(grid)
library(rstan)
library(cowplot)
#Bogaerts
path <- "http://baseball.yahoo.co.jp/mlb/teams/player/fielder/stats/3050735"
tab1 <-read_html(path) %>>%
  html_table()

tmp <-lapply(tab1[4:6], function(x){n <-nrow(x);x[-c((n-1):n),]})
dat1 <-do.call("rbind",tmp)
dat1$Date <-as.Date(dat1$"日付",format = "%m/%d")
N <- nrow(dat1)
dat1 <- dat1[N:1,]
Bogaerts <- list(T=N,n=dat1$"打数",y=dat1$"安打")

モデル

このようなデータに当てはまりそうなのは打数(n_i)を試行回数、安打数(y_i)を成功回数とした二項分布です。

ところで、世間ではよく打者には「調子の波」というものが存在すると言われています。

ヒットが多い「好調」時や、アウトが多い「不調」時があるとして、それをモデル化するにはどうしたらよいでしょうか。

一つのアイデアですが、以下のように考えることで、打率の変化をモデルに取り込むことができます。

y_i \sim {\rm binomial}(p_i,n_i)
\mu_i \sim {\rm Normal}( \mu_{i-1},\sigma^2)

ここで p_i = {\rm logit}^{-1}(\mu_i) です。

このモデルでは打率 p_i が徐々に変化していきます。

しかし打者には本当に「調子の波」があるのでしょうか。偶然ヒットが固まったりアウトが固まったりするだけではないか、と考えることもできます。

調子の波が存在しないという仮説に基づくと以下のようなモデルが考えられます。

y_i \sim {\rm binomial}(p,n_i)

ここでは p は常に一定としています。

調子の波が存在しないモデルと調子の波が存在するモデル、どちらのほうが良いモデルかを WAIC で比べてみよう、というのが本エントリの趣旨です。

コーディング

調子の波が存在しないモデルの Stan コードは以下のようになりました。

data { 
  int<lower=0> T;
  int<lower=0> n[T];
  int<lower=0> y[T];
} 
parameters {
  real<lower=0, upper=1> p;
}
model {
    y ~ binomial(n,p);
}
generated quantities{
  real log_lik[T];
  for(t in 1:T){
   log_lik[t] =binomial_lpmf(y[t]|n[t],p); 
  }
}

一方、調子の波が存在するモデルは以下のようになりました。

data { 
  int<lower=0> T;
  int<lower=0> n[T];
  int<lower=0> y[T];
} 
parameters {
  real mu[T];
#  real beta;
  real<lower=0> sigma;
}
model {
#  sigma ~ exponential(1);
  for(t in 2:T){
    mu[t] ~ normal(mu[t-1],sigma);
  }
  for(t in 1:T){
    y[t] ~ binomial_logit(n[t],mu[t]); 
  }
}
generated quantities{
  real theta[T];
  real log_lik[T];
  for(t in 1:T)
    theta[t] =inv_logit(mu[t]);
  for(t in 1:T)
    log_lik[t] =binomial_logit_lpmf(y[t]|n[t],mu[t]);
}

generated quantities ブロックでは WAIC を算出するための対数尤度を計算しています。

Stan コードを実行して、WAIC を計算する R コードは以下のようになりました。

WAIC のコードは、WAICとWBICを事後分布から計算する - StatModeling Memorandum をまんま使っています。

rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
binomial_logit <-stan_model("binomial_logit.stan")
binomial0 <-stan_model("binomial0.stan")
fit_Bogaerts <-sampling(binomial_logit,Bogaerts,iter=10000,thin=5)
fit_Bogaerts0 <-sampling(binomial0,Bogaerts,iter=10000,thin=5)
#traceplot(fit_Bogaerts,pars=c("sigma","mu[1]","mu[2]",
#                             "mu[3]","mu[4]","mu[5]"))
#traceplot(fit_Bogaerts0,pars="p")

waic <- function(la){
  log_likelihood <- la$log_lik
  training_error <- - mean(log(colMeans(exp(log_likelihood))))
  functional_variance_div_N <- mean(colMeans(log_likelihood^2) - colMeans(log_likelihood)^2)
  waic <- training_error + functional_variance_div_N
  return(waic)
}

round(waic(extract(fit_Bogaerts)),2)
#[1] 1.31
round(waic(extract(fit_Bogaerts0)),2)
#[1] 1.32

結果、調子の波が存在するモデルの WAIC は 1.31、調子の波が存在しないモデルの WAIC は 1.32 でした。

ボガーツには調子の波が存在すると考えたほうがよさそうです。

図示

打率の変化を以下のように図示してみました。

f:id:abrahamcow:20160812222557p:plain

青い折れ線グラフが打率の変化、薄い青の帯が打率の95%ベイズ信頼区間です。

下の積み上げ棒グラフが日ごとのヒットとアウトの回数を示しています。

このグラフを書くためのコードは以下の通りです。

phat1 <-get_posterior_mean(fit_Bogaerts,"theta")[,"mean-all chains"]
lower1 <- apply(extract(fit_Bogaerts)$theta, 2, quantile,0.025)
upper1 <- apply(extract(fit_Bogaerts)$theta, 2, quantile,0.975)

Bogaerts_df  <- data.frame(Date=dat1$Date,out=dat1$"打数"-dat1$"安打",hit=dat1$"安打")
write.csv(Bogaerts_df,file="Bogaerts2016.csv",row.names = FALSE)
dat_g <- tidyr::gather(Bogaerts_df,result,value,-Date)

p1 <- ggplot(dat_g,aes(x=Date,y=value,fill=result))+
  geom_bar(stat = "identity")+
  theme(legend.position="bottom")+
  xlab("")+ylab("")+
  labs(fill="")+
  scale_x_date(date_labels = "%m/%d")+
  scale_fill_grey()
p2 <- ggplot(dat1,aes(x=Date,y=phat1))+
  geom_line(colour="royalblue",size=1)+
  geom_ribbon(aes(ymin=lower1,ymax=upper1),fill="royalblue",alpha=0.2)+
  xlab("")+ylab("")+
  ylim(c(0,1))+
  theme(axis.title.x = element_blank(), axis.text.x = element_blank(),
        axis.ticks = element_blank(),axis.line.x = element_blank())

theme_set(theme_cowplot(font_family = "HiraMaruProN-W4"))
grid.newpage()
grid.draw(rbind(ggplotGrob(p2), ggplotGrob(p1), size = "first"))

付録

成形したデータはこっそり以下において置きます。

Bogaerts2016.csv · GitHub