伊庭先生から教わったことを超速で復習するよ@Stan

超速復習CAR＠酔っ払い

空間構造とは系列に並んだデータが前後（時間的，地理的）のデータに影響を受けるモデルだよ

require( dplyr)

## Loading required package: dplyr
## 
## Attaching package: 'dplyr'
## 
## The following object is masked from 'package:stats':
## 
##     filter
## 
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

require( rstan)

## Loading required package: rstan
## Loading required package: Rcpp
## Loading required package: inline
## 
## Attaching package: 'inline'
## 
## The following object is masked from 'package:Rcpp':
## 
##     registerPlugin
## 
## rstan (Version 2.4.0, packaged: 2014-07-21 17:16:17 UTC, GitRev: c995bc8b9d67)

まずは正解の生成

N <- 100
r <- rnorm(N,0,0.2) %>%
  cumsum()
r <- r - min(r)
lambda <- exp(0.1 + r)

y <- rpois(N,lambda)
dat.list <- list(N=N, y=y)

plot(lambda,type='l')

シンプルなアイディアでは前後の影響を受ける場合 1回階差のシステムノイズを考えるよ．

stancode <- '
data{
  int N;
  int y[N];
}
parameters{
  real<lower=0> beta;
  real<lower=0> sigma;
  vector<lower=0>[N] r;
}
transformed parameters{
  vector[N] lambda;
  lambda <- exp(beta + r);
}
model{
  beta ~ normal(0,10);
  r[1] ~ normal(0,10);
  for(n in 2:N){
    r[n] ~ normal(r[n-1],sigma);
  }

  y ~ poisson(lambda);
}
'

fit <- stan(model_code = stancode, data = dat.list, chain = 1,iter=1000)

## 
## TRANSLATING MODEL 'stancode' FROM Stan CODE TO C++ CODE NOW.
## COMPILING THE C++ CODE FOR MODEL 'stancode' NOW.
## 
## SAMPLING FOR MODEL 'stancode' NOW (CHAIN 1).
## 
## Iteration:   1 / 1000 [  0%]  (Warmup)
## Iteration: 100 / 1000 [ 10%]  (Warmup)
## Iteration: 200 / 1000 [ 20%]  (Warmup)
## Iteration: 300 / 1000 [ 30%]  (Warmup)
## Iteration: 400 / 1000 [ 40%]  (Warmup)
## Iteration: 500 / 1000 [ 50%]  (Warmup)
## Iteration: 501 / 1000 [ 50%]  (Sampling)
## Iteration: 600 / 1000 [ 60%]  (Sampling)
## Iteration: 700 / 1000 [ 70%]  (Sampling)
## Iteration: 800 / 1000 [ 80%]  (Sampling)
## Iteration: 900 / 1000 [ 90%]  (Sampling)
## Iteration: 1000 / 1000 [100%]  (Sampling)
## #  Elapsed Time: 10.2326 seconds (Warm-up)
## #                5.70777 seconds (Sampling)
## #                15.9404 seconds (Total)

lambda_p <- extract(fit,pars='lambda')$lambda %>%
  apply(2,median)
plot(lambda_p,type='l')

でもCARでは2回の階差，つまり変化率が推移するモデルを 考えていているよ

stancode2 <- '
data{
  int N;
  int y[N];
}
parameters{
  real<lower=0> beta;
  real<lower=0> sigma;
  vector<lower=0>[N] r;
}
transformed parameters{
  vector[N] lambda;
  lambda <- exp(beta + r);  
}
model{

  r[1] ~ normal(0,10);
  r[2] - r[1] ~ normal(0,10);
  for(n in 3:N){
    r[n] - r[n-1] ~ normal(r[n-1]-r[n-2],sigma);
  }
  y ~ poisson(lambda);

  beta ~ normal(0,10);
}
'

fit2 <- stan(model_code = stancode2, data = dat.list, chain = 1,iter=1000)

## 
## TRANSLATING MODEL 'stancode2' FROM Stan CODE TO C++ CODE NOW.
## COMPILING THE C++ CODE FOR MODEL 'stancode2' NOW.
## 
## SAMPLING FOR MODEL 'stancode2' NOW (CHAIN 1).
## 
## Iteration:   1 / 1000 [  0%]  (Warmup)
## Iteration: 100 / 1000 [ 10%]  (Warmup)
## Iteration: 200 / 1000 [ 20%]  (Warmup)
## Iteration: 300 / 1000 [ 30%]  (Warmup)
## Iteration: 400 / 1000 [ 40%]  (Warmup)
## Iteration: 500 / 1000 [ 50%]  (Warmup)
## Iteration: 501 / 1000 [ 50%]  (Sampling)
## Iteration: 600 / 1000 [ 60%]  (Sampling)
## Iteration: 700 / 1000 [ 70%]  (Sampling)
## Iteration: 800 / 1000 [ 80%]  (Sampling)
## Iteration: 900 / 1000 [ 90%]  (Sampling)
## Iteration: 1000 / 1000 [100%]  (Sampling)
## #  Elapsed Time: 182.874 seconds (Warm-up)
## #                150.867 seconds (Sampling)
## #                333.74 seconds (Total)

lambda_p2 <- extract(fit2,pars='lambda')$lambda %>%
  apply(2,median)

plot(lambda_p2,type='l')

この場合推移は1回階差に比べて滑らかになるよ．

どちらがいいかは，モデリングする対象次第だから適宜えらぼうね．

個人的に受け取ったメッセージはこんな感じ

だれかタイムシフト予約した方，見直したいいいいい． @piroyoungまで！！！！