Help us understand the problem. What is going on with this article?

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

More than 3 years have passed since last update.

サンプルデータ

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

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

モデル

神様の方法を真似して、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         
Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
Comments
No comments
Sign up for free and join this conversation.
If you already have a Qiita account
Why do not you register as a user and use Qiita more conveniently?
You need to log in to use this function. Qiita can be used more conveniently after logging in.
You seem to be reading articles frequently this month. Qiita can be used more conveniently after logging in.
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
ユーザーは見つかりませんでした