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

廿TT

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

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

モチベーション

野球にはまったく興味ないんだけど、Albert (2008) "Streaky Hitting in Baseball" を読んでた。

http://citeseerx.ist.psu.edu/viewdoc/download?rep=rep1&type=pdf&doi=10.1.1.150.5808

スポーツデータ解析に関する論文を探すには Journal of Quantitative Analysis in Sports という雑誌が良さそう。

下記はカルロス・ギーエンという選手の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)

年間通しての打率は約 3 割 2 分。

(ybar <- mean(GuillenC)) #batting average
#0.3203593

上記のデータの30打席ごとの移動平均をとってみる。

ma <- function(x, n){
  m <- stats::filter(x, rep(1,n)) / n
  m <- as.numeric(m)
  m[!is.na(m)]
}
m <-ma(GuillenC,30)
####
#moving average graph for Carlos Guillen.
m2 <-c(ybar,m,ybar)
#plot.ts(m2,ylab="moving average", main="2005 Carlos Guillen") 
#abline(h=ybar,lty=2)
library(ggplot2)
df4plot <- data.frame(time=1:length(m2), MA=m2)
df4plot$lower <- ifelse(m2<=ybar,m2,ybar)
ggplot(df4plot,aes(x=time,y=MA)) +
  geom_line() +
  geom_polygon(aes(x=time, ymin=lower,ymax=MA),alpha=0.3) +
  geom_hline(yintercept=ybar) +theme_bw()

f:id:abrahamcow:20160406012937p:plain

こうして見てみるとある時期には 5 割を超えていたり、またある時期には 1 割を下回っていたりする。

これはギーエンという打者に調子の波が存在する証拠だ、と言っていいだろうか。

パラメトリックブートストラップによる検定

調子の波が存在しない選手は、常にコンスタントな打率でヒットを出すから、その打撃成績を上記のように 0 か 1 かに符号化すると、それはベルヌーイ過程になる。

このコンスタントバッターを乱数でシミュレーションしてみる。

n <- length(GuillenC)
set.seed(1)
Chitter <- rbinom(n,size=1,prob=ybar)
m_C <-ma(Chitter,30)
m2_C <-c(ybar,m_C,ybar)
df4plot_C <- data.frame(time=1:length(m2_C), MA=m2_C)
df4plot_C$lower <- ifelse(m2<=ybar,m2,ybar)
ggplot(df4plot_C,aes(x=time,y=MA)) +
  geom_line() +
  geom_polygon(aes(x=time, ymin=lower,ymax=MA),alpha=0.3) +
  geom_hline(yintercept=ybar) +theme_bw()

f:id:abrahamcow:20160406225927p:plain

これはこれで意味ありげな波が出来てしまった。

こうなってくるとちゃんと検定したくなるのが人情だろう。

帰無仮説:
ギーエンの0-1のプロセスがベルヌーイ過程である。”

対立仮説:
ギーエンの0-1のプロセスがベルヌーイ過程でない。”

帰無分布は乱数で生成するので思いつく限りの統計量で仮説検定を試すことができる。

  • アウトの連長の最大値(the length of the longest run of failures)
  • ヒットの連長の最大値(the length of the longest run of successes)
  • アウトの平均連長(the mean length of the lengths of runs of failures)
  • ヒットの平均連長(the mean length of the lengths of runs of success)

ここで 連(run)とは 0 や 1 が連続するひとかたまりのことである。

モチベーションとなった例で出した移動平均からも統計量を作ってみる。

  • 移動平均のレンジ(the range of the moving averages)  R = \max_j m_j − \min_j m_j
  • 移動平均のシーズン平均からの平均変動(the mean variation of the moving averages about the season average)

\displaystyle B =\frac{1}{n-w-1}\sum^{n-w+1}_{j=1}|m_j-\bar y|

この B は冒頭のグラフの黒く塗りつぶした部分の面積に比例する量である。

(R <- diff(range(m)))
#0.5
(B <-sum(abs(ma(GuillenC,30)-mean(GuillenC)))/length(m)) #"black" statistic
#####
runs <-rle(GuillenC)
(lengthRunOf0 <-  max(runs$lengths[runs$value==0]))
#19
(lengthRunOf1 <-  max(runs$lengths[runs$value==1]))
#4
(averageRunOf0 <- mean(runs$lengths[runs$value==0]))
#3.242857
(averageRunOf1 <- mean(runs$lengths[runs$value==1]))
#1.528571
########
#simulation
boot_lengthRunOf0 <- boot_lengthRunOf1 <-
  boot_averageRunOf0 <- boot_averageRunOf1 <-
  boot_B <- boot_R<- numeric(10000)
for(i in 1:10000){
  simv <-rbinom(n,1,ybar)
  runs_sim <-rle(simv)
  boot_lengthRunOf0[i] <-  max(runs_sim$lengths[runs_sim$value==0])
  boot_lengthRunOf1[i] <-  max(runs_sim$lengths[runs_sim$value==1])
  boot_averageRunOf0[i] <- mean(runs_sim$lengths[runs_sim$value==0])
  boot_averageRunOf1[i] <- mean(runs_sim$lengths[runs_sim$value==1])
  boot_m <-ma(simv,30)
  boot_ybar <-mean(simv)
  boot_B[i] <-sum(abs(boot_m-boot_ybar))/length(boot_m)
  boot_R[i] <- diff(range(boot_m))
}
hist(boot_lengthRunOf0)
abline(v=lengthRunOf0,col="red2",lwd=2)
hist(boot_lengthRunOf1)
abline(v=lengthRunOf1,col="red2",lwd=2)
hist(boot_averageRunOf0)
abline(v=averageRunOf0,col="red2",lwd=2)
hist(boot_averageRunOf1)
abline(v=averageRunOf1,col="red2",lwd=2)
hist(boot_B)
abline(v=B,col="red2",lwd=2)
hist(boot_R)
abline(v=R,col="red2",lwd=2)

検定統計量の分布をヒストグラムで示した。

赤い線はギーエンの打席結果から求めた統計量の実現値である。

f:id:abrahamcow:20160409014902p:plain
f:id:abrahamcow:20160409014900p:plain

f:id:abrahamcow:20160409014908p:plain
f:id:abrahamcow:20160409014859p:plain

f:id:abrahamcow:20160409014911p:plain
f:id:abrahamcow:20160409014904p:plain

それぞれの経験 p-値は以下のようになった。

> sum(boot_lengthRunOf0>=lengthRunOf0)/10000
[1] 0.063
> sum(boot_lengthRunOf1>lengthRunOf1)/10000
[1] 0.5318
> sum(boot_averageRunOf0>averageRunOf0)/10000
[1] 0.3266
> sum(boot_averageRunOf1>averageRunOf1)/10000
[1] 0.2612
> sum(boot_R>R)/10000
[1] 0.0133
> sum(boot_B>B)/10000
[1] 0.0069
統計量 p-値
アウトの連長の最大値 0.063
ヒットの連長の最大値 0.5318
アウトの平均連長 0.3266
ヒットの平均連長 0.2612
R 0.0133
B 0.0069

B だけが 5%水準で有意となった。

この結果からギーエンには調子の波があると言えそうである。

分割表の独立性の検定

ギーエンが直前の打席の結果をひきずっているとしたら、アウトの直後にヒットを打つ割合と、ヒットの直後にヒットを打つ割合は変化することになる。

直前\直後 アウト ヒット
アウト 157 70
ヒット 69 37

この分割表に対して独立性の検定を行う

帰無仮説:
”直前の結果と直後の結果は独立である。”

対立仮説:
”直前の結果と直後の結果は独立でない。”

検定してみると p-値は十分に大きく、ギーエンは直前の結果をひきずっているとは考えにくい。

> chisq.test(tab1)

	Pearson's Chi-squared test with Yates' continuity correction

data:  tab1
X-squared = 0.3778, df = 1, p-value = 0.5388

つづく

統計量 B からはギーエンには調子の波があるといえそうなことはわかった。

B は打率に依存するため、B を使って選手間で調子の波の強さを比較したりはできない。

そこでベータ二項モデルを導入してギーエンのヒットの数をモデル化することを考える。

abrahamcow.hatenablog.com

関連エントリ

abrahamcow.hatenablog.com

アマゾンアフィリエイトのコーナー

アルバート先生のご本は邦訳も出ています。

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

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

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

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

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

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

ベルヌーイ過程の連長の分布はフェラーによってイグザクトに求められていた気がする。

確率論とその応用 1 上

確率論とその応用 1 上