## ３章のポアソン回帰モデルB（ゆとり度Xを加える)
library(rstudioapi)  ##Packages に☑

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

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

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

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

model {
  for (i in 1:n) {
    Y[i] ~ poisson(theta[i]);
  }
  beta0 ~ normal(0, 100);
  beta1 ~ normal(0, 100);
  r ~ normal(0, s);
}
'
# Stanの実行
library(rstan)
fitb <- stan(model_code=PoissonRegB,data=data,iter=2000,chains=4,
            pars=c('beta0','beta1', 'r', 's'), seed=23456,control=list(adapt_delta=0.9))
fitb
traceplot(fitb, pars=c('beta0','beta1'))