## ３章のポアソン回帰モデルA（個人差を表す)
library(rstudioapi)  ##Packages に☑

d <- read.csv('Table3_1.csv')
n <- nrow(d)
Y <- d$Y
data <- list(n=n, Y=Y)

## stanのモデルコード
PoissonRegA <-'
data {
  int<lower=0> n;       // number of customer
  int<lower=0> Y[n];    // number of ryokou
}

parameters {
  real beta0;
  real r[n];
  real<lower=0> s;
}

transformed parameters {
  real<lower=0> theta[n];
  for (i in 1:n) {
    theta[i] = exp(beta0 + r[i]);
  }
}

model {
  for (i in 1:n) {
    Y[i] ~ poisson(theta[i]);
  }
  beta0 ~ normal(0, 100);
  r ~ normal(0, s);
}
'
# Stanの実行
library(rstan)
fita <- stan(model_code=PoissonRegA,data=data,iter=2000,chains=4,
            pars=c('beta0', 'r', 's'), seed=12345)
fita

#　図3.3のトレース・プロット
traceplot(fita, pars=c('beta0','r[1]','s'))

