11
9

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 5 years have passed since last update.

{rstan} 欠損値あり時系列データのモデル推定

Last updated at Posted at 2017-03-02

サンプルデータ

適当なサンプルデータを用意します。

ランダムにデータの半分をブッコヌキます。

モデル

神様の方法を真似して、2階差分を仮定します。

\begin{eqnarray}
&\mu_t - \mu_{t-1} &\sim& normal(\mu_{t-1} - \mu_{t-2},\ s_{mu})\\
\leftrightarrow &\mu_t &\sim& normal(2*\mu_{t-1} - \mu_{t-2},\ s_{mu})
\end{eqnarray}

stan code

えーと、方法は2つあります。

・個々のIDごとの補完を行いたい場合
  →  神の書の9.5.2 欠測値を参照。
・全体の傾向を抽出したい場合
  → 大元のマニュアルのChapter 15を参照。

ここでは後者。
まんまなので、コードを載せて終わらせますが、$segment$を使うのがポイントです。

系列idごとのtimeが揃っていないという困難を回避する為の細工です。
これは、いわゆる「縦長データ」(R業界ではtidyなデータと呼びますが)を作っておいて、グループごとのサンプルサイズ(変数Sに格納)と共に処理するものです。

データYはベクトルで与えられ、Sに従って1:N_timeに分解されるという理屈です。

  data{
    int N_sample;
    int N_time;

    vector[N_sample] Y;
    int S[N_time];
  }

  parameters{
    real<lower = 0> s;
    real<lower = 0> s_mu;
    real mu[N_time];
  }

  model{
    int pos;
    pos = 1;
    for(k in 1:N_time){
      segment(Y, pos, S[k]) ~ normal(mu[k], s);
      pos = pos + S[k];
    }

    for(k in 3:N_time){
      mu[k] ~ normal(2 * mu[k-1] - mu[k-2], s_mu);
    }
  }

結果

スクリーンショット 2017-03-02 18.14.53.png

Appendix

サンプルデータ

set.seed(111)

N_time <- 24
N_id <- 10
SD <- 0.3

x <- 1:N_time
y <- sin(x/N_time * 2 * pi)

dat <- data.frame(time = x, id = factor(rep(1:N_id, each = N_time)),
                  value = rnorm(N_time * N_id, y, SD))

# 一様乱数を使ってデータを半分にしています。
dat1 <- dat %>% 
  dplyr::mutate(index = runif(N_time * N_id, 0, 1)) %>% 
  dplyr::filter(index > 0.5)

stan キックコード

datastan <- list(N_sample = nrow(dat1), N_time = N_time,
                 Y = dat1$value[order(dat1$time)], 
                 S = as.numeric(table(dat1$time)))

# stanコードはテキストにしてstanmodelという変数に格納しておきます。
fit <- stan(model_code = stanmodel, data = datastan, seed = 123)

お絵かき

サンプルデータの図

ggplot(dat, aes(time, value))+
  theme_classic()+
  geom_path(aes(color = id), alpha = 0.5)+
  geom_point(aes(color = id))

stan の結果の図

色々な描き方があるので、お好みで。

dat_g <- tidy(summary(fit)$summary) %>% 
  filter(.rownames %in% paste("mu[", 1:N_time, "]", sep = "")) %>% 
  mutate(time = as.numeric(str_sub(.rownames, start = 4, end = -2)))


ggplot(dat_g, aes(x = time, y = mean))+
  theme_classic()+
  geom_ribbon(aes(ymax = X97.5., ymin = X2.5.), fill = "lightgrey", alpha =0.5)+
  geom_ribbon(aes(ymax = X75., ymin = X25.), fill = "darkgrey", alpha =0.5)+
  geom_line()

環境

> sessionInfo()
R version 3.3.1 (2016-06-21)
Platform: x86_64-apple-darwin13.4.0 (64-bit)
Running under: OS X 10.10.5 (Yosemite)

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] stringr_1.2.0        broom_0.4.1          ggmcmc_1.1          
[4] tidyr_0.6.0          rstan_2.14.1         StanHeaders_2.14.0-1
[7] ggplot2_2.2.1        dplyr_0.5.0         
11
9
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
11
9

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?