循環する変数の統計モデリング

周期性のある変数・循環する変数を含むモデリングを実践しましたので紹介します。スライドは埋め込んで、ソースコードのコピペ＆解説をメインにします。

とある病んだ院生の体内時計（サーカディアンリズム） from . .

使用したデータは以下。自由に使ってください。

元データ: data.txt
起床時刻だけ抜き出したもの: data_sleepout.txt
起床時刻から起床時刻の時間: data_out2out.txt
起床時刻の累積時刻: data_cum.txt

まずはじめに起床時間の分布のパラメータを求めたStanコードは以下になります（model1.stan）。

data {
   int<lower=0> N;
   real Sleep_out[N];
}
parameters {
   real<lower=0, upper=2*pi()> mu;
   real<lower=0, upper=100> sigma;
}
model {
   for (i in 1:N)
      Sleep_out[i] ~ von_mises(mu, sigma);
}

7行目・11行目、Stanでvon Mises分布を使うときはsigmaをある程度小さくしましょう。

kickするRコードは以下になります。

library(rstan)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())

d <- read.delim("input/data_sleepout.txt", sep="\t")

data <- list(
   N = nrow(d),
   Sleep_out = d$sleep.out/24*2*pi
)

fit <- stan(file='model/model1.stan', data=data, iter=1000)
save.image("output/result1.Rdata")

2～3行目： Stan 2.7.0からこの記述をするだけでstan関数が並列実行されます。
9行目: von Mises分布に渡すために0-24を0-2πに変換しています。

次に起床時刻から起床時刻への時間のモデルです。Stanのhello worldともいうべき簡潔さ。

data {
   int<lower=0> N;
   real<lower=0> Out2out[N];
}
parameters {
   real<lower=0> mu;
   real<lower=0, upper=100> sigma;
}
model {
   for (i in 1:N)
      Out2out[i] ~ normal(mu, sigma);
}

これの変化点検出バージョン。

data {
   int<lower=0> N;
   real<lower=0> Out2out[N];
}
parameters {
   real<lower=0> mu;
   real<lower=0, upper=100> sigma;
   real<lower=-20, upper=20> m_d[2];
   real<lower=-20, upper=20> s_d[2];
   real<lower=30, upper=80> c1;
   real<lower=80, upper=180> c2;
   real<lower=350, upper=450> c3;
}
model {
   for (i in 1:N){
      real m;
      real s;
      m <- mu    + int_step(i - c1) * int_step(c2 - i) * m_d[1]
                 + int_step(i - c3) * m_d[2];
      s <- sigma + int_step(i - c1) * int_step(c2 - i) * s_d[1]
                 + int_step(i - c3) * s_d[2];
      Out2out[i] ~ normal(m, s);
   }
}

8行目: 工事期間やスッキリ期間は周期が変わると思われるのでその差分です。
9行目:　同じく分散が変わると思われるのでその差分です。
10行目: 工事期間の始まり。
11行目:　工事期間の終わり。
12行目: スッキリ期間の始まり。この変化点たちはもう少し範囲を狭めないと推定が厳しいようです。
18～21行目: int_step()はその引数が0より大きい時1、それ以外0をとるStanの関数です。工事が終わると元に戻ると仮定してこんな形になっています。

最後に累積時刻と周期を同時にモデル化したStanコードです（model4.stan）。

data {
   int<lower=0> N;
   real Cum_hour[N];
   real Z_lower[N];
   real Z_upper[N];
}
parameters {
   real<lower=0, upper=10> s_noise;
   real<lower=0, upper=50> s_z;
   real<lower=0, upper=50> mu;
   real<lower=-5, upper=5> z_pre[N];
}
transformed parameters{
   real<lower=0> z[N];
   for (i in 1:N)
      z[i] <- Z_lower[i] + (Z_upper[i] - Z_lower[i]) * inv_logit(z_pre[i]);
}
model {
   for (i in 1:N)
      Cum_hour[i] ~ normal(z[i], s_noise);
   for (i in 2:N)
      z[i] ~ normal(z[i-1] + mu, s_z);
}

4～5行目: 潜在変数zの範囲はデータとして与えます。 11行目: z_preもある程度範囲を狭めておかないと、たまにとんでもない大きな（もしくは小さな）値をとって抜け出せなくなります。BUGSの時も同様の小手先のテクニックがありました。
13～17行目: z_pre[i] を [Z_lower[i], Z_upper[i]] の値に変換しています。インデックスによって範囲が異なる場合はこのように設定するしか方法がないと思います。
21～22行目: 潜在変数zはARモデルで時間変化します。muが周期に相当します。

実行するRコードは以下になります。

library(rstan)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())

d <- read.delim("input/data_cum.txt", sep="\t")
cum_hour <- d$cum.hour

data <- list(
   N = nrow(d),
   Cum_hour = cum_hour,
   Z_lower = cum_hour - 24,
   Z_upper = cum_hour + 24
)

fit <- stan(file='model/model4.stan',data=data, iter=60000, thin=30)
save.image("output/result4.Rdata")

11～12行目: 潜在変数zの範囲をデータとして与えています。さすがに前後24時間取ればそれ以上ずれることないという常識からです。もっと狭めてもよいと思います。

StatModeling Memorandum

StanとRとPythonでベイズ統計モデリングします. たまに書評.

循環する変数の統計モデリング