setwd("C:/xxx")		#作業ディレクトリの指定

mtopdata <- read.csv("mtopdata.csv")
mtopdata_lst <- list(T=nrow(mtopdata),Y=mtopdata$kouryo,sqGRP=mtopdata$sqrt_GRP,c=61)

stancode <- '
data {
	int T;				//データの時点数
	real Y[T];			//目的変数（ブランド考慮）
	real sqGRP[T];			//説明変数（GRPのルート）
	int c;				//商品リニューアル後の最初の時点
}

parameters {
	real m[T];			//1時点前に依存して変化するm
	real b[T];			//GRPのルートの係数
	real<lower=0> s_e;		//yの観測方程式の標準偏差
	real<lower=0> s_rm;		//mの状態方程式の標準偏差
	real<lower=0> s_rb;		//bの状態方程式の標準偏差
}

model {
	for (t in 1:T) {
		Y[t] ~ normal(m[t]+b[t]*sqGRP[t], s_e);	//尤度関数
		}
	for (t in 2:T) {
		m[t] ~ normal(m[t-1], s_rm);		//mの事前分布
		}
	for (t in 2:c-1) {
		b[t] ~ normal(b[t-1], s_rb);		//bの事前分布（商品リニューアル前)
		}
	for (t in c+1:T) {
		b[t] ~ normal(b[t-1], s_rb);		//bの事前分布（商品リニューアル後)
		}
	m[1] ~ normal(0,100);				//mの1時点目の事前分布
	b[1] ~ normal(0,10);				//bの1時点目の事前分布
	b[c] ~ normal(0,10);				//bの商品リニューアル後の最初の時点の事前分布
	s_e ~ uniform(0,100);				//s_eの事前分布
	s_rm ~ uniform(0,100);				//s_rmの事前分布
	s_rb ~ uniform(0,10);				//s_rbの事前分布
}
'

library(rstan)
rstan_options(auto_write=TRUE)
options(mc.cores=parallel::detectCores())

stanmodel <- stan_model(model_code=stancode)

fit <- sampling(stanmodel, data=mtopdata_lst,
		init=function(){
			list(m=runif(nrow(mtopdata),0,100), b=runif(nrow(mtopdata),0,1),
			s_e=runif(1,0,10),s_rm=runif(1,0,10),s_rb=runif(1,0,0.1))
		},
		iter=20500, warmup=500, chains=4, seed=12345)

write.table(summary(fit)$summary,file="outputsummary.txt",sep="\t",quote=FALSE,col.names=NA)
traceplot(fit,pars=c("b[1]","s_rb"),inc_warmup=TRUE,window=c(1,500),nrow=2)

save(fit, file="mtop_fit.dat")
mcsample <- extract(fit,permuted=TRUE)
write.table(mcsample,file="mcsample.txt",sep="\t",quote=FALSE,col.names=NA)
