状態空間モデルにはどれほどデータが必要か

この記事はStan Advent Calendar2017のエントリー記事です。

Stan盛り上がってますかー。盛り上がってますね（断定）。
私は本当に，Stanのおかげで生き返ったように楽しくモデリング出来ています。

先日の日本心理学会でもベイジアンモデリングのチュートリアルが大にぎわいでした。
それともう一つ人気だったな，と思うのが「時系列データ分析」とか「経験サンプリング」の話だったような印象です。
最近は加速度計とかモーションキャプチャ，アイトラックレコーダなどで，大量のデータが取れるようになってきているので，心理学者的には「わぁい隅々まで行動を計測できる！」というあたりが嬉しかったりするわけです。

でもまあ，大量に，時系列的に取れたデータはそれにあった分析が必要。特に

• 自己相関が高い；一時点，あるいはそれ以上前の数値に関係のある現時点のデータ
• 誤差まみれのデータ；無意味な情報も含んでいる大量のデータから意味のある部分だけ取り出す必要がある

という二つの問題から，「適切なモデリング」が必要になってきます。
そうすると，データ生成メカニズムを考えているベイジアンモデリングとの相性がいい。特に状態空間モデルは，ベイジアンアプローチをした方が，素直にモデリングできるのでお勧めです。

ということで，「いよいよ心理学にも空間統計とかはいってくるようになるかぁ」なんて思ったりしますけど，ちょっと待って。
心理学業界の一般的なデータは，上記のような器具を使うのでなければ，2，3時点の銃弾的データがほとんどで，1,000や10,000点以上あるでーたなんて（今のところ）あんまり見かけない。
どれぐらいデータがあればちゃんとした推定ができるんですかね？とちょっと気になりました。

と言うことで，一次元の状態空間モデル（時系列的なデータ）の仮想データを作り，そのデータサイズがどれぐらいだったら推定値として使えるのかのチェックをしてみようと思い立ちました。感度分析っていうのかな？こういうのも。

状態空間モデルについては過去記事もご参照いただくとして，味噌は目に見えない「状態」をモデリングし，それに観測誤差がついてデータになる，と考えること。すなわち，状態をS,観測値をXとすると，

$$S^{t+1}=S^t + \sigma$$
$$X^{t} = S^t + \phi$$

とします。あとは適切な事前分布をおいてやればオッケー。Stanのコードも簡単ですよ。

data{
  int<lower=1> L;//data length
  vector[L] X;
}

parameters{
  vector[L] r;
  real<lower=0> sig;
  real<lower=0> ER;
}

model{
  //obj
  X ~ normal(r,ER);
  //state
  r[2:L] ~ normal(r[1:(L-1)],sig);
  //prior
  ER ~ cauchy(0,5);
  sig ~ cauchy(0,5);
}

ほとんどイメージのまま書けますし，なめら化とかもしなくていいので楽です（わざとです。念のため。）。
さてさて，どれほどのデータがいるのかを検証するために，データサイズLをいろいろ変えたサンプルデータを作り，このコードで推定させて見ましょう。シミュレーションなので各試行を100回やってその平均値を考えることにします。

データセットを作り，MED推定値を返してくる関数を作りました。データサイズNと変動sig，測定誤差ERを渡すと

# 一次元状態空間モデルのデータを作る関数
extract.fit <- function(N,sig,ER){
  r <- rep(0,N)
  r[1] <- 10
  for(i in 2:N){
    r[i] = r[i-1] +rnorm(1,0,sig)
  }
  err <- rnorm(N,0,ER)
  datastan <- list(L=N,X=r+err)
  fit <- rstan::sampling(model,datastan)
  fit.sig <- median(rstan::extract(fit,pars="sig")$sig)
  fit.ER  <- median(rstan::extract(fit,pars="ER")$ER)
  return(list(fit.sig=fit.sig,fit.ER=fit.ER))
}

これをコンパイルして走らせます。データサイズは指数的に増えるようにコード化しました。ここでは，Nは3,4,5,7,9,12,16,22,30,40,55,74,99,134,181,245,330,446,602,812,1097と増えて行きます。

# サンプルサイズ
check <- c(round(exp(seq(1,7,0.3))))
# 検証回数
N <- 100

# 検証
result.matrix <- c()
for(P in check){
  results <- list()
  for(n in 1:N){
    pc <- proc.time()[1]
    ret <- extract.fit(P,3,5)
    time <- proc.time()[1]-pc
    result.matrix=rbind(result.matrix,c(P=P,n=n,sig=ret$fit.sig,ER=ret$fit.ER,time=time))
  }
}

これでシャラーンと走らせた結果がこんな感じ。

状態の変動分，sigの真値は3で，観測誤差ERの真値は5です。実際のスコアで見るとこんな感じになります。

P	key	mean	sd
3	sig	4.468877	1.6121974
4	sig	4.052790	1.3807442
5	sig	3.839792	1.2260976
7	sig	3.789735	1.3442228
9	sig	3.543699	1.3567688
12	sig	3.566722	1.3576163
16	sig	3.345746	1.4561492
22	sig	3.323947	1.3837509
30	sig	3.113396	1.1905917
40	sig	3.142804	0.9679240
55	sig	3.101363	0.9104340
74	sig	3.176283	0.7182003
99	sig	3.017480	0.5714553
134	sig	2.997657	0.5463738
181	sig	2.974278	0.5059406
245	sig	2.983619	0.3664567
330	sig	3.006930	0.2936964
446	sig	2.997985	0.2609966
602	sig	3.069371	0.2744286
812	sig	3.001113	0.1897300
1097	sig	3.022984	0.1875880

だいたい3桁（99か134)あたりで真値に辿り着いている感じ。40，55点ぐらいだと真値+0.1ぐらいで，SDも1.0ぐらいあるので当たってるけどちょっと精度悪いと言う感じでしょうか。30以下はちょっと外れすぎかな。もちろんそれぞれの領域においてどの程度の精度が必要か，と言うのは変わってくると思いますけどね。ちなみに測定誤差の方についてもほぼ同様。

P	key	mean	sd
3	ER	4.214914	1.6879550
4	ER	4.045415	1.6517713
5	ER	3.827050	1.3597881
7	ER	4.267874	1.5982695
9	ER	4.256110	1.3975502
12	ER	4.378329	1.2951766
16	ER	4.406097	1.2772903
22	ER	4.751362	1.0473943
30	ER	4.736674	0.9558007
40	ER	4.907256	0.8717817
55	ER	4.811194	0.7367819
74	ER	4.859863	0.5464713
99	ER	4.949720	0.5134745
134	ER	4.942098	0.4262254
181	ER	5.004040	0.4054680
245	ER	4.991647	0.3249265
330	ER	4.955008	0.2834097
446	ER	4.970939	0.2508226
602	ER	4.961162	0.1931686
812	ER	4.991295	0.1552512
1097	ER	4.980203	0.1451442

それぞれの変動分を大きくしたり小さくしたりして，試して見ることができます。

いや，実は自分の持っているデータでなんとかしようと思ったのですが，測定が35点ぐらいしかなくて，うまくいってなかったので「じゃあ何点ならできるんだよ！」と思って作ったコードでした。
100点を超えるデータって，例えばアンケート調査とかではほぼ無理だろうと思います。心理屋さんに状態空間モデルが流行するのがいつになるでしょうか。

最後に手前味噌ながら，宣伝を。こちらの本「Doing Bayesian Data Analysis」は翻訳・監修で関わらせていただきました。StanのメカニズムやMCMCのアルゴリズムについて，おそらく日本で一番簡単な図入りで説明している本です。お子様のクリスマスプレゼントにいかがでしょうか。

犬4匹本サポートサイトでは正誤表や関係イベントなどの情報を発信しております。

もしちょっと専門性の高いませたお子様，お友達がいると言うのであれば，こちらの通称「コワい本」を贈りましょう。これを枕の下に入れて眠ると，赤くて怖い人たちがパラメータを推定せよと迫ってくる夢が見られますよ。素敵！

ツイッタのハッシュタグ，「#犬4匹本」「#コワい本」でも情報発信しておりますのでぜひご活用ください。