RstanでMCMCの練習

R MCMC

$ cat poi.stan

data{
int J; // number of data
int y[J];
}
parameters{
real mu;
}

model{
y[J] ~ poisson(mu);
}

平均３０のポアソン分布

###Poison
library("rstan")
y<-rpois(100,lambda=30)
data <- list(J = 100, y =y)
fit <- stan(file = 'poi.stan', data = data,iter = 1000,chains = 4)
print(fit)

出力結果：

Inference for Stan model: poi.
4 chains, each with iter=1000; warmup=500; thin=1;
post-warmup draws per chain=500, total post-warmup draws=2000.

mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
mu 33.61 0.21 5.57 23.95 29.55 33.31 37.09 45.33 683 1
lp__ 81.92 0.02 0.65 80.17 81.77 82.17 82.34 82.38 948 1

Samples were drawn using NUTS(diag_e) at Sat Mar 12 10:44:44 2016.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).

  • 次に平均がランダムに変化するポアソン過程の平均と分散の推定

$ cat poi4.stan

data{
  int<lower=0> J; // number of schools
  real y[J]; // estimated treatment effects
}

parameters {
  real mu1;
  real<lower=0.000> sd1;
}

transformed parameters{
}

model{
  #lambda ~ poisson(mu1);
      for(j in 1:J){
        increment_log_prob(normal_log(lambda1[j],mu1,sd1));
      }
      for(j in 1:J){
         increment_log_prob(poisson_log(y[j],lambda1[j]));

        # increment_log_prob(poisson_log(y[j],mu1));
      }
  #lambda ~ normal(mu1,mu1*k1);
 # y ~ poisson(mu1);
}

＊尤度表示をつかっている．
increment_log_prob は尤度を足し込んでいく関数．

Rコード

###RD process
y<-rpois(100,rnorm(100,mean=100,sd=20))
data <- list(J = 100, y =y)
fit <- stan(file = 'poi_4.stan', data = data,iter = 1000,chains = 4)
traceplot(fit)
plot(fit)

出力

> print(fit)
Inference for Stan model: poi_3.
4 chains, each with iter=1000; warmup=500; thin=1;
post-warmup draws per chain=500, total post-warmup draws=2000.

mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
mu1 103.19 0.05 1.70 99.92 102.04 103.22 104.27 106.49 1049 1.00
sd1 17.74 0.04 1.34 15.42 16.79 17.65 18.52 20.72 1044 1.01
lp__ -425.94 0.04 1.05 -428.84 -426.31 -425.62 -425.23 -424.97 574 1.01

Samples were drawn using NUTS(diag_e) at Sun Mar 13 09:40:51 2016.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).

##RD model2
N<-30
x<-1:30
inits<-list(list(alpha=5,beta=50),list(alpha=5,beta=50),list(alpha=5,beta=50),list(alpha=5,beta=50))
para <- c("alpha", "beta")
a<-rnorm(30,mean=100,sd=100*0.01)
y<-rpois(length(a),lambda=a)
data<-list(y=y,N=N)

ll<-stan(file="test_6.stan",data=data,iter=1000,chain=4,init=inits,par=para)

plot(ll)
traceplot(ll)
ll

data {

int N;
#real x[N];
int y[N];

}


parameters {
real alpha;
real beta;


}

model{
#real mu[N];

alpha ~ normal(1,0.01);
#beta ~ normal(5,0.1);
for(n in 1:N){
#mu[n] ~ normal(alpha+beta*x[n],0.1);
y[n] ~ poisson(alpha*beta);
}


}

出力：

mean se_mean sd 2.5% 25% 50% 75% 97.5%
alpha 1.00 0.00 0.01 0.98 0.99 1.00 1.01 1.02
beta 97.77 0.08 2.05 93.90 96.33 97.75 99.11 101.96
lp__ 10506.72 0.04 1.00 10504.18 10506.34 10507.03 10507.41 10507.69
n_eff Rhat
alpha 835 1.00
beta 672 1.00
lp__ 608 1.01

##RD model6
N<-300
x<-1:30
#inits<-list(list(alpha=5,beta=50),list(alpha=5,beta=50),list(alpha=5,beta=50),list(alpha=5,beta=50))
para <- c("alpha", "beta","sigma","delta")
a<-rnorm(N,mean=1,sd=0.01)
c[1]<-30
for(i in 2:N){
c[i]<-0.8*c[i-1]+rnorm(N,mean=0,sd=30*0.1)+30
}

y<-rpois(length(a),lambda=c*a)
data<-list(y=y,N=N)

#ll<-stan(file="test_7.stan",data=data,iter=1000,chain=4,init=inits,par=para)
ll<-stan(file="test_8.stan",data=data,iter=1000,chain=4,par=para)


plot(ll)
traceplot(ll)

data {

int N;
#real x[N];
int y[N];

}


parameters {
real delta;
real alpha;
real beta;
real sigma;
real c[N];

}

model{
#real mu[N];

delta ~ normal(1,0.01);
for(n in 2:N){
#mu[n] ~ normal(alpha+beta*x[n],0.1);
c[n] ~ normal(alpha+beta*c[n-1],sigma);
y[n] ~ poisson(c[n]*delta);
}


}

出力

Inference for Stan model: test_8.
4 chains, each with iter=1000; warmup=500; thin=1;
post-warmup draws per chain=500, total post-warmup draws=2000.

mean se_mean sd 2.5% 25% 50% 75%
alpha 37.15 0.29 4.81 27.86 33.92 37.11 40.28
beta 0.75 0.00 0.03 0.69 0.73 0.75 0.77
sigma 3.69 0.11 0.83 2.13 3.13 3.61 4.15
delta 1.00 0.00 0.01 0.98 0.99 1.00 1.01
lp__ 179097.66 11.35 62.01 178981.88 179056.96 179096.99 179134.35
97.5% n_eff Rhat
alpha 46.88 274 1.01
beta 0.81 274 1.01
sigma 5.53 57 1.14
delta 1.02 233 1.00
lp__ 179245.10 30 1.17

Samples were drawn using NUTS(diag_e) at Thu Mar 17 22:45:25 2016.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).

data {

int<lower=0> N;

#real x[N];
real c;
int<lower=0> y[N];

}


parameters {

real<lower=0> delta[N];

real alpha;

real beta;

real<lower=0> sigma;
real<lower=0> gamma;
real<lower=0> r[N];
//real r[N];


}

model{

#real mu[N];



for(n in 2:N){
delta[n] ~ normal(1,gamma);
#mu[n] ~ normal(alpha+beta*x[n],0.1);

r[n] ~ normal(alpha+beta*r[n-1],sigma);

y[n] ~ poisson(c*r[n]*delta[n]);

}

###RD model####
N<-300

x<-1:30

#inits<-list(list(alpha=5,beta=50),list(alpha=5,beta=50),list(alpha=5,beta=50),list(alpha=5,beta=50))

para <- c("alpha", "beta","sigma","delta")

a<-rnorm(N,mean=1,sd=0.01)
c1<-NULL
r1<-NULL
r1[1]<-1

for(i in 2:N){

	r1[i]<-0.8*r1[i-1]+rnorm(N,mean=0,sd=0.05)+1

}

c1<-30
y<-rpois(length(a),lambda=c1*r1*a)

data<-list(y=y,N=N,c=30)

#ll<-stan(file="test_7.stan",data=data,iter=1000,chain=4,init=inits,par=para)

#ll<-stan(file="test_8.stan",data=data,iter=1000,chain=4,par=para)
ll<-stan(file="test_8b.stan",data=data,iter=1000,chain=4,par=para)
  • 参考

超要約 Stan Reference
https://rstudio-pubs-static.s3.amazonaws.com/31748_bff68e5646e74bbeb267c61beba9f812.html　

デバッグモードなど設定の切り替え
http://blog.kz-md.net/?p=935

RStanのおさらいをしながら読む 岩波DS 1 http://rpubs.com/uri-sy/iwanami_ds1