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

廿TT

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

[KFAS]0-1データの状態空間モデル(打者の調子の波のモデル化 3)

とりあえずプロット

  1. Albert (2008) 打者の調子の波のモデル化 - 廿TT
  2. 打者の調子の波のモデル化(幾何分布編) - 廿TT

と同じく、カルロス・ギーエンの打撃成績のデータを使います。

y <- 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)

ヒットを 1、アウトを 0 としています。

以下のように loess でスムージングすると、

library(cowplot)
n <-length(y)
ggplot()+
  geom_point(shape=3,aes(x=1:n,y=y))+
  geom_smooth(aes(x=1:n,y=y))

こうなりました。

だいたい 30 打席が一週間くらいだそうです。

f:id:abrahamcow:20160612040948p:plain

ヒットが固まっている時期やアウトが固まっている時期があるように見えます。

このような調子の波のモデルを考えます。

KFAS パッケージで状態空間モデル

R の KFAS パッケージでは、ポアソン分布、二項分布など正規分布以外の分布を扱えます。

冒頭の 0-1 データ y に対して、以下のようなモデルを仮定します。

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

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

打率がロジットスケールでランダムウォークすることにしています。

以下のように推定すると、

library(KFAS)
model1 <- SSModel(y~SSMtrend(degree = 1, Q = list(NA)), 
                  u=1,distribution="binomial")

fit1 <- fitSSM(model1,inits = 1, method="Brent", lower=-100,upper=100)
out1 <- KFS(fit1$model, smoothing = c("state", "mean"))

ggplot()+
  geom_point(aes(x=1:n,y=y),shape=3)+
  geom_line(aes(x=1:n,y=fitted(out1)),size=1)+
  geom_hline(aes(yintercept=mean(y)),linetype=2)

こうなりました。

f:id:abrahamcow:20160612042319p:plain

点線は年間を通じた平均の打率です

95%信頼区間も合わせてプロットするにはこうします。

out1_2 <-predict(fit1$model,interval ="confidence")

library(cowplot)
ggplot()+
  geom_point(aes(x=1:n,y=y),shape=3)+
  geom_line(aes(x=1:n,y=out1_2[,"fit"]),size=1)+
  geom_ribbon(aes(x=1:n,ymin=out1_2[,"lwr"],ymax=out1_2[,"upr"]),alpha=0.3)+
  geom_hline(aes(yintercept=mean(y)),linetype=2)

f:id:abrahamcow:20160612115447p:plain

議論

教えてください。

予測

このモデルは打率の予測に使えるかもしれません。

一期先の打率は次のように求めます。

> predict(fit1$model,n.ahead = 1)
Time Series:
Start = 335 
End = 335 
Frequency = 1 
           fit
[1,] 0.3219524

ランダムウォークなので点推定値は直前の打率と一緒ですが、これが年間を通じた打率を用いるより、良い予測値になっていたらいいなと思いました。

調子の波の定量化

 \sigma の推定値は、

> exp(fit1$optim.out$par)
[1] 0.004501138

です。

この値の大小によって、選手ごとに調子の波の激しさを比較できないでしょうか?

モデル比較

そもそも調子の波というようなものは本当に存在するのでしょうか。

打率が一定と仮定したモデル(単純なベルヌーイ試行)と、AIC かなにかで較べてこのモデルが勝ったら、「調子の波は存在する」と言えそうです。

この場合 AIC はどうやって求めるのでしょうか。

自由パラメータの数は分散パラメータの数だけでいいのかな?

model0 <- SSModel(y~SSMtrend(degree = 1, Q = list(0)), 
                  u=1,distribution="binomial")

-2*logLik(fit1$model)+2
#[1] 420.7668

-2*logLik(model0)
#[1] 421.3783

追記

2階のトレンドモデルもやってみた。対数尤度はローカルレベルモデルのほうが大きいようです。

model2 <- SSModel(y~SSMcustom(Z=matrix(c(1,0),1,2),T=matrix(c(2,-1,1,0),2,2,byrow = TRUE),Q=NA),
                  u=1,distribution ="binomial")
fit2 <- fitSSM(model2,inits = 1, method="Brent", lower=-100,upper=100)
out2 <- KFS(fit2$model, smoothing = c("state", "mean"))
out2_2 <-predict(fit2$model,interval ="confidence")
ggplot()+
  geom_point(aes(x=1:n,y=y),shape=3)+
  geom_line(aes(x=1:n,y=out2_2[,"fit"]),size=1)+
  geom_ribbon(aes(x=1:n,ymin=out2_2[,"lwr"],ymax=out2_2[,"upr"]),alpha=0.3)+
  geom_hline(aes(yintercept=mean(y)),linetype=2)

f:id:abrahamcow:20160622030544p:plain