japanr.r

setwd("E:/TokyoR/JapanR2015")

library(rstan)

#stanファイル名
scr <- "mitaronstandraft.stan"

d <- read.csv("sample.csv")

#ユニークな個人の情報
d1 <- d[!duplicated( data.frame(d$id) ), ]
#ファイル内の変数をリストにまとめる
data = list(N=nrow(d1),
            NT=nrow(d),
            y = d$離反,
            x = cbind(rep(1,length=nrow(d)),matrix(scale(d$価格),nrow=nrow(d), ncol=1)),
            w = cbind(
                      matrix(scale(d$ucm),nrow=nrow(d), ncol=1),
                      matrix(scale(d$scm),nrow=nrow(d), ncol=1),
                      matrix(scale(d$qcm),nrow=nrow(d), ncol=1),
                      matrix(scale(d$kcm),nrow=nrow(d), ncol=1),
                      matrix(scale(d$pcm),nrow=nrow(d), ncol=1),
                      matrix(scale(d$usns),nrow=nrow(d), ncol=1),
                      matrix(scale(d$ssns),nrow=nrow(d), ncol=1),
                      matrix(scale(d$qsns),nrow=nrow(d), ncol=1),
                      matrix(scale(d$ksns),nrow=nrow(d), ncol=1),
                      matrix(scale(d$psns),nrow=nrow(d), ncol=1)
                    ),
            id = d$id,
            z = cbind(
　　　　　　　　　　　matrix(scale(d1$年齢),nrow=nrow(d1), ncol=1),
　　　　　　　　　　　matrix(scale(d1$性別),nrow=nrow(d1), ncol=1),
                      matrix(scale(d1$未既婚),nrow=nrow(d1), ncol=1),
                      matrix(scale(d1$世帯年収),nrow=nrow(d1), ncol=1),
                     )
            )


#サンプリング結果を出力する変数名
par <- c("alpha","beta")

war <- 10 #バーンイン期間
ite <- 100 #繰り返し回数
see <- 1234 #乱数の種類
dig <- 3 #有効数字
cha <- 1 #連鎖構成数

#並列計算
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())

#モデル
fit <- stan(file="mitaronstandraft.stan",
            data=data,
            pars=c("alpha","beta"),
            verbose=F,
            seed=71,
            chains=1,
            warmup=1000,
            iter=10000)

#shinystan
library(shinystan)
sso_normal <- launch_shinystan(fit)
  ##DIAGNOSE ：収束診断
  ##ESTIMATE ：推定結果
  ##EXPLORE ：素敵なグラフ
  ##MORE：Stanコード・単語の説明など?

#結果の出力

fit.summary <- data.frame(summary(fit)$summary)
write.csv(fit.summary, "fit_summary1.csv", quote=F)

#summary <- summary(fit)
#write.csv(summary, "summary.csv",quote=FALSE)

pdf("fit_plot1.pdf", width=600/72, height=600/72)
plot(fit)
dev.off()

pdf("fit_traceplot1.pdf", width=600/72, height=600/72)
traceplot(fit)
dev.off()

library("coda")
fit.coda<-mcmc.list(lapply(1:ncol(fit),function(x) mcmc(as.array(fit)[,x,])))
pdf("fit_density1.pdf", width=600/72, height=600/72)
plot(fit.coda)
dev.off()

#library (reshape)
#library (ggplot2)
#post   <- extract (fit, permuted = F)
#m.post <- melt (post)
#graph  <- ggplot (m.post, aes(x = value))
#graph  <- graph + geom_density () + facet_grid(. ~ parameters, scales = "free") + theme_bw() #facet_wrapのほうがいいかも
#pdf("fit_density2.pdf", width=600/72, height=600/72)
#plot (graph)