この記事はStan Advent Calendar 2017の12月8日の記事です。
Stanの盛り上がりに一役も二役も買った，「StanとRでベイズ統計モデリング」の，読書会があちこちで行われておりますが，去る11月18日に関西学院大学で行われたOsaka.stanも読書会。この本の11，12章でした。この最後の二つの章は，ちょっと複雑で難しいモデリングでもあって，上手に解説していただいたおかげで自分でもできるかも，という気になりますね。
さて，12章の状態空間モデル，これについては私も去年のアドカレや今年のアドカレでもやってるんですが，まあ手習い程度で複雑なモデリングはしてないのです。で，読書会でスピーカーの @masashikomori 先生が「まあ一階の状態空間による予測って，フエルトも減るとも言わずにわからねーっていう信用区間を示しますから，あんまり意味がないですよね。二階差分ぐらいしないとね」というわけです。ま，たしかにそうだと。式を見てもそんなに難しくないわけです。だからやってみようかな，と思ってこのエントリーもしてみました。
二階差分は「差分の差分」なので，「傾きを滑らかにする」ようにモデリングするから，方向性（上昇，下降）をもった予測ができるよ，と。三階になると加速度になるので現実的な空間モデルではよくわかんないとのこと。たしかにそうです。式としてはややこしく見えるんだけど，
$$ u^t – u^{t-1} = u^{t-1} – u^{t-2} + e $$
から
$$u^t = u^{t-1} + u^{t-1} – u^{t-2} + e$$
$$ u^t = 2u^{t-1} – u^{t-2} + e$$
となるだけで，簡単にできるのです。
さて，私は毎朝体重と体脂肪を測定するという趣味を持っていますから，そのデータを今回も使います。データそのものはかなり長くあるんだけど，これの今年度部分を使いたいと思います。2017年1月1日からの体重変化をモデリング！
(今回使ったデータも含めて公開中です。おヒマな方はどうぞ。)
コードはこちら。
観測が欠損値や予測の箇所は，データに9999というありえない数字を入れて，if文で分岐させるようにしています。
data {
int<lower=0> n; //データ総数
vector[n] W; //体重データ
int<lower=0> Nmiss; //欠損値の数
}
parameters {
real muZero; //左端
vector[n] mu; //確率的レベル
vector<lower=0>[Nmiss] Miss_W; //欠測値はパラメータ
real<lower=0> sdV; //観測誤差の大きさ
real<lower=0> sdW; //過程誤差の大きさ
}
model {
// 状態方程式の部分
// 左端から初年度の状態を推定する
mu[1] ~ normal(muZero, sdW);
mu[2] ~ normal(muZero, sdW);
// 観測方程式の部分
{
int nmiss;
nmiss = 0;
for(i in 1:n) {
if(W[i]!=9999){
W[i] ~ normal(mu[i], sdV);
}else{
nmiss = nmiss+1;
Miss_W[nmiss] ~ normal(mu[i],sdV);
}
}
}
// 状態の遷移
for(i in 3:n){
mu[i] ~ normal(2*mu[i-1]-mu[i-2], sdW);
}
sdV ~ cauchy(0,5);
sdW ~ cauchy(0,5);
}
これで予測したのが次の図。
帯の濃い色は50%信用区間，薄い部分は95%信用区間。
なるほど，予測線が下を向いてます。これでできてるし，体重は下降傾向。よしよし。
・・・ちょっと待てよ。
薄々気づいてはいたんだけど，最近体重がちょっと増え気味じゃないかなぁ。ベスト体重は66kgですよ。そして結婚してからこっち，それを超えて76kgで安定しているのは知っていますよ。それにしても80kgの大台に乗るのはよろしくないなあ。これ，どこかで「デブモード(Dev Mode?)」に入ってしまったのではないか。これを考えるモデルはないだろうか。
ということで，色々考えて見ました。
1. デブモードモデル
読書会11章で読んだ知識を活用します。俺には「普通モード」と「デブモード」の二種類があって，デブモードのときは普通モードよりも体重が高い水準で推移するとします。普通モードになるか，デブモードになるかは確率$$\pi$$のベルヌーイ分布(0/1)になるとします。そんなモデルがこちら。
dataブロックは同じなので省略
parameters {
...中略...
real<lower=0,upper=1> pi; //モード変換率
real<lower=0> beta;
}
model {
...ここまで同じなので省略...
// 状態の遷移
for(i in 3:n){
real tmp;
tmp = (2*mu[i-1]-mu[i-2]);
target += log_sum_exp(
bernoulli_lpmf(0|pi) + normal_lpdf(mu[i]|tmp,sdW),
bernoulli_lpmf(1|pi) + normal_lpdf(mu[i]|tmp+beta,sdW)
);
}
sdV ~ cauchy(0,5);
sdW ~ cauchy(0,5);
beta ~ cauchy(0,5);
}
betaは下限0，つまりデブモードは増えている前提で，どれぐらい増えているかわからないので裾の思いコーシー分布を置いてみました。結果はこちら。
Inference for Stan model: dietDouble_mode.
4 chains, each with iter=10000; warmup=5000; thin=1;
post-warmup draws per chain=5000, total post-warmup draws=20000.
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
sdW 0.02 0.00 0.01 0.01 0.01 0.02 0.02 0.03 179 1.03
sdV 0.41 0.00 0.02 0.37 0.40 0.41 0.43 0.46 991 1.01
beta 25.32 4.15 532.99 0.00 0.40 3.21 9.25 95.70 16512 1.00
pi 0.04 0.00 0.14 0.00 0.00 0.00 0.01 0.54 2403 1.00
Samples were drawn using NUTS(diag_e) at Wed Nov 22 09:48:17 2017.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).
ちょっと収束が悪いところもありますが，betaの中央値が3.21kg,EAPが25kgですな。半コーシーだからだいぶん歪んでます。まあデブモードが平均25キロ増ってのは言い過ぎでしょうよ，このモードに入ったら3キロ太るのか。なるほど。で，このモードに入る確率は4%ぐらい。うーん，滅多にないことなんだけど，入ったら3キロ太るのかー，気をつけないと。
とは思ったんだけど，なんかシックリ来ないんですよねえ。
2. ただ切り替わるんじゃなくて，日付の関数だよね
このモデルだと，ただの確率で二つのモードがコロコロ切り替わるわけです。が，俺が考えたかったのは「あれ？ある日を基準にデブモードに入ってるんじゃない？」ということなのでした。これはどう考えたらいいか。
変化点モデルとか，隠れマルコフモデルとかが該当するんだろうか，とちょっと悩んだんですが，もう少し単純に表現できるのではないかと。
そこで考えたのが，ベルヌーイ分布の確率$$\pi$$に構造を入れること。イメージはIRTの1PLモデルで，要するに逆ロジット（ロジスティック）関数の中に困難度母数bを入れてS字カーブを右に左にずらす，あんな感じです。つまりある日を境に変化率$$\pi$$が１の方に上がっていくイメージ。
例えばXデーが300日目だったらこんなカーブになるわけです。
ということで，そんなモデルを書いてみました。ここで二階差分モデルを捨てたことに注意。
dataブロックは同じなので省略
parameters {
...中略...
real<lower=2,upper=n> tau; //変化日
real<lower=0> beta;
}
model {
...ここまで同じなので省略...
// 状態の遷移
for(i in 3:n){
target += log_sum_exp(
bernoulli_logit_lpmf(0|i-tau) + normal_lpdf(mu[i]|mu[i-1],sdW),
bernoulli_logit_lpmf(1|i-tau) + normal_lpdf(mu[i]|mu[i-1]+beta,sdW)
);
}
sdV ~ cauchy(0,5);
sdW ~ cauchy(0,5);
beta ~ cauchy(0,5);
}
これでやってみた結果がこちら。
Inference for Stan model: dietDouble_mode2.
4 chains, each with iter=10000; warmup=5000; thin=1;
post-warmup draws per chain=5000, total post-warmup draws=20000.
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
sdW 0.17 0.00 0.03 0.12 0.15 0.17 0.19 0.23 740 1.00
sdV 0.35 0.00 0.02 0.30 0.34 0.35 0.37 0.40 1224 1.00
beta 2.51 0.20 7.18 0.01 0.38 1.54 3.33 9.90 1282 1.01
tau 310.75 2.34 60.67 106.95 324.20 331.12 338.25 343.53 674 1.00
Samples were drawn using NUTS(diag_e) at Tue Nov 21 18:54:28 2017.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).
大変結構。tauの50%HDIが324から338というところ。どうやらデブモードに入ったのは10月20日〜11月3日の頃だ！ということがわかったわけです。
3.じゃあもうinv_logitの関数でbetaをかけてやったらいいんじゃね？
で，ここまでやって思いましたけど，これ普通に変化日までのS字カーブに増加率betaをかけるモデルにした方がシンプルなんじゃないかな？と思ってやって見ました。
...ここまで同じなので省略...
model {
...ここまで同じなので省略...
// 状態の遷移
for(i in 2:n){
mu[i] ~ normal(mu[i-1]+beta*inv_logit(i-tau),sdW);
}
beta ~ cauchy(0,5);
}
こちらで分析したのがこちら。
Inference for Stan model: dietDouble_mode3.
4 chains, each with iter=10000; warmup=5000; thin=1;
post-warmup draws per chain=5000, total post-warmup draws=20000.
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
sdW 0.17 0.00 0.03 0.12 0.15 0.17 0.19 0.23 551 1.01
sdV 0.35 0.00 0.03 0.31 0.34 0.35 0.37 0.40 1286 1.00
beta 14.07 1.55 56.84 0.45 2.14 4.72 10.74 91.60 1337 1.00
tau 330.83 0.66 6.66 322.63 324.25 329.81 336.44 343.19 101 1.05
Samples were drawn using NUTS(diag_e) at Wed Nov 22 10:12:18 2017.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).
ちょっとtauのRhatが悪いので，微調整が必要かもですが，やはりXデーは330日目の付近なわけです。うーん，そうかぁ。。。
ということがわかったからといって，なんの解決にもならないことに，数時間自分の体重を見つめ続けたあとで気づいたような次第です。
来年はデブモードを脱出するぞ！なんならヘルスサイドに堕ちるぞ！
みなさんMerry Christmas and Enjoy Stan!