読者です 読者をやめる 読者になる 読者になる

廿TT

譬如水怙牛過窓櫺 頭角四蹄都過了 因甚麼尾巴過不得

[RStan]変曲点のある単回帰分析

豊田秀樹『マルコフ連鎖モンテカルロ法』5章8節の例題。

マルコフ連鎖モンテカルロ法 (統計ライブラリー)

マルコフ連鎖モンテカルロ法 (統計ライブラリー)

変曲点のある回帰直線、

 \displaystyle y_i =   \begin{cases} \alpha + \beta_1(x_i - x_{\rm change}) & x_i \le x_{\rm change}\\ \alpha + \beta_2(x_i - x_{\rm change}) & x_i > x_{\rm change} \end{cases}

を attenu データに当てはめる。

Stan のコードはこんな感じ。

data {
	int<lower=0>N;
	real<lower=0>x[N];
	real<lower=0>y[N];
}
parameters{
	real<lower=0> alpha;
	real beta1;
	real beta2;
	real<lower=x[1], upper=x[N]> changepoint;
	 real<lower=0> sigma;
}
transformed parameters  {
   real mu[N];
   for (i in 1:N) {
   	if(x[i] <= changepoint){
   	mu[i] <- alpha + beta1*(x[i] - changepoint); 
   	}
   	else{
   	mu[i] <- alpha + beta2*(x[i] - changepoint);	
   	}
   }
}
model{
	for (i in 1:N) {
	y[i] ~ normal(mu[i], sigma); 
	}
}

R のコードはこう。

library(rstan)
data(attenu)
model1 <- stan_model(file="changepoint.stan")
dat4stan <- list(N=nrow(attenu),x=attenu$dist,y=attenu$accel)
fit <- sampling(model1, data=dat4stan,iter=10000)

ex <-extract(fit)
changepoint <- mean(ex$changepoint)
beta1 <- mean(ex$beta1)
beta2 <- mean(ex$beta2)
alpha <- mean(ex$alpha)
f1 <-function(x){
  alpha + ifelse(x<=changepoint,beta1,beta2)*(x-changepoint)
}
plot(accel~dist,data=attenu,col="grey20")
curve(f1(x),add=TRUE,lwd=3,col="red",n=1000)

traceplot(fit,pars=c("alpha","beta1","beta2","changepoint","sigma"))

f:id:abrahamcow:20150908172213p:plain

f:id:abrahamcow:20150908172202p:plain

各パラメータの95%信用区間。

> quantile(ex$alpha,c(0.025, 0.975))
      2.5%      97.5% 
0.05914026 0.13664998 
> quantile(ex$beta1,c(0.025, 0.975))
       2.5%       97.5% 
-0.02881408 -0.01019488 
> quantile(ex$beta2,c(0.025, 0.975))
         2.5%         97.5% 
-0.0007918284 -0.0001460663 
> quantile(ex$sigma,c(0.025, 0.975))
      2.5%      97.5% 
0.09199037 0.11367734 
> quantile(ex$changepoint,c(0.025, 0.975))
    2.5%    97.5% 
13.21038 29.26548 

abrahamcow.hatenablog.com