#matplot(x,ZB_new,type="l",lwd=2,ylim=c(-2,2),xlab="x",ylab="")
matplot(x,ZB_new-5,type="l",lwd=2,add=TRUE)
abline(h=-4)
dev.off()
res_mu
res
ls(res)
?extract
x
mu
setwd("/Users/zhenkewu/Dropbox/ZW/professional/teaching/graphical\ models/notes/lecture14/Stan")
fit <- stan('model.stan',
data=c("M","C","ZB","N"),
iter=10000,warmup=5000,thin=5, chains=1)
print(fit)
res <- extract(fit,permuted=TRUE)
res_mu <- mu
png(file.path(result_folder,"estimated_mean_curve.png"),
width=10,height=8,units = "in",res=72)
par(mar=c(5,5,3,3))
matplot(x,t(res_mu),type="l",col=2,main=paste0("Mean curve"),
xlab="x",ylab="y",
ylim=c(-6,4))
points(x,y)
points(x,mu,col="black",lty=1,lwd=4,type="l")
#par(mar=c(4,5,0,3))
beta_postmean <- colMeans(res$beta)
ZB_new <- sweep(ZB,MARGIN = 2,STATS = beta_postmean,"/")
#matplot(x,ZB_new,type="l",lwd=2,ylim=c(-2,2),xlab="x",ylab="")
matplot(x,ZB_new-5,type="l",lwd=2,add=TRUE)
abline(h=-4)
dev.off()
res
matplot(x,res_mu,type="l",col=2,main=paste0("Mean curve"),
xlab="x",ylab="y",
ylim=c(-6,4))
res <- extract(fit,permuted=TRUE)
res_mu <- mu
png(file.path(result_folder,"estimated_mean_curve.png"),
width=10,height=8,units = "in",res=72)
par(mar=c(5,5,3,3))
matplot(x,res_mu,type="l",col=2,main=paste0("Mean curve"),
xlab="x",ylab="y",
ylim=c(-6,4))
points(x,y)
points(x,mu,col="black",lty=1,lwd=4,type="l")
#par(mar=c(4,5,0,3))
beta_postmean <- colMeans(res$beta)
ZB_new <- sweep(ZB,MARGIN = 2,STATS = beta_postmean,"/")
#matplot(x,ZB_new,type="l",lwd=2,ylim=c(-2,2),xlab="x",ylab="")
matplot(x,ZB_new-5,type="l",lwd=2,add=TRUE)
abline(h=-4)
dev.off()
res <- extract(fit,permuted=TRUE)
res_mu <- mu
png(file.path(result_folder,"estimated_mean_curve.png"),
width=10,height=8,units = "in",res=72)
par(mar=c(5,5,3,3))
matplot(x,res_mu,type="l",col=2,main=paste0("Mean curve"),
xlab="x",ylab="y",
ylim=c(-6,4))
points(x,y)
points(x,mu,col="black",lty=1,lwd=4,type="l")
#par(mar=c(4,5,0,3))
beta_postmean <- colMeans(res$beta)
ZB_new <- sweep(ZB,MARGIN = 2,STATS = beta_postmean,"/")
#matplot(x,ZB_new,type="l",lwd=2,ylim=c(-2,2),xlab="x",ylab="")
matplot(x,ZB_new-5,type="l",lwd=2,add=TRUE)
abline(h=-4)
dev.off()
setwd("/Users/zhenkewu/Dropbox/ZW/professional/teaching/graphical\ models/notes/lecture14/Stan")
fit <- stan('model.stan',
data=c("M","C","ZB","N"),
iter=10000,warmup=5000,thin=5, chains=1)
print(fit)
res <- extract(fit,permuted=TRUE)
res_mu <- res$mu
png(file.path(result_folder,"estimated_mean_curve.png"),
width=10,height=8,units = "in",res=72)
par(mar=c(5,5,3,3))
matplot(x,t(res_mu),type="l",col=2,main=paste0("Mean curve"),
xlab="x",ylab="y",
ylim=c(-6,4))
points(x,y)
points(x,mu,col="black",lty=1,lwd=4,type="l")
#par(mar=c(4,5,0,3))
beta_postmean <- colMeans(res$beta)
ZB_new <- sweep(ZB,MARGIN = 2,STATS = beta_postmean,"/")
#matplot(x,ZB_new,type="l",lwd=2,ylim=c(-2,2),xlab="x",ylab="")
matplot(x,ZB_new-5,type="l",lwd=2,add=TRUE)
abline(h=-4)
dev.off()
setwd("/Users/zhenkewu/Dropbox/ZW/professional/teaching/graphical\ models/notes/lecture14/Stan")
fit <- stan('model.stan',
data=c("M","C","ZB","N"),
iter=10000,warmup=5000,thin=5, chains=1)
print(fit)
res <- extract(fit,permuted=TRUE)
res_mu <- res$mu
png(file.path(result_folder,"estimated_mean_curve.png"),
width=10,height=8,units = "in",res=72)
par(mar=c(5,5,3,3))
matplot(x,t(res_mu),type="l",col=2,main=paste0("Mean curve"),
xlab="x",ylab="y",
ylim=c(-6,4))
points(x,y)
points(x,mu,col="black",lty=1,lwd=4,type="l")
#par(mar=c(4,5,0,3))
beta_postmean <- colMeans(res$beta)
ZB_new <- sweep(ZB,MARGIN = 2,STATS = beta_postmean,"/")
#matplot(x,ZB_new,type="l",lwd=2,ylim=c(-2,2),xlab="x",ylab="")
matplot(x,ZB_new-5,type="l",lwd=2,add=TRUE)
abline(h=-4)
dev.off()
setwd("/Users/zhenkewu/Dropbox/ZW/professional/teaching/graphical\ models/notes/lecture14/Stan")
fit <- stan('model.stan',
data=c("M","C","ZB","N"),
iter=10000,warmup=5000,thin=5, chains=1)
print(fit)
res <- extract(fit,permuted=TRUE)
res_mu <- res$mu
png(file.path(result_folder,"estimated_mean_curve.png"),
width=10,height=8,units = "in",res=72)
par(mar=c(5,5,3,3))
matplot(x,t(res_mu),type="l",col=2,main=paste0("Mean curve"),
xlab="x",ylab="y",
ylim=c(-6,4))
points(x,y)
points(x,mu,col="black",lty=1,lwd=4,type="l")
#par(mar=c(4,5,0,3))
beta_postmean <- colMeans(res$beta)
ZB_new <- sweep(ZB,MARGIN = 2,STATS = beta_postmean,"/")
#matplot(x,ZB_new,type="l",lwd=2,ylim=c(-2,2),xlab="x",ylab="")
matplot(x,ZB_new-5,type="l",lwd=2,add=TRUE)
abline(h=-4)
dev.off()
setwd("/Users/zhenkewu/Dropbox/ZW/professional/teaching/graphical\ models/notes/lecture14/Stan")
fit <- stan('model.stan',
data=c("M","C","ZB","N"),
iter=10000,warmup=5000,thin=5, chains=1)
print(fit)
res <- extract(fit,permuted=TRUE)
res_mu <- res$mu
png(file.path(result_folder,"estimated_mean_curve.png"),
width=10,height=8,units = "in",res=72)
par(mar=c(5,5,3,3))
matplot(x,t(res_mu),type="l",col=2,main=paste0("Mean curve"),
xlab="x",ylab="y",
ylim=c(-6,4))
points(x,y)
points(x,mu,col="black",lty=1,lwd=4,type="l")
#par(mar=c(4,5,0,3))
beta_postmean <- colMeans(res$beta)
ZB_new <- sweep(ZB,MARGIN = 2,STATS = beta_postmean,"/")
#matplot(x,ZB_new,type="l",lwd=2,ylim=c(-2,2),xlab="x",ylab="")
matplot(x,ZB_new-5,type="l",lwd=2,add=TRUE)
abline(h=-4)
dev.off()
setwd("/Users/zhenkewu/Dropbox/ZW/professional/teaching/graphical\ models/notes/lecture14/Stan")
fit <- stan('model.stan',
data=c("M","C","ZB","N"),
iter=10000,warmup=5000,thin=5, chains=1)
print(fit)
res <- extract(fit,permuted=TRUE)
res_mu <- res$mu
png(file.path(result_folder,"estimated_mean_curve.png"),
width=10,height=8,units = "in",res=72)
par(mar=c(5,5,3,3))
matplot(x,t(res_mu),type="l",col=2,main=paste0("Mean curve"),
xlab="x",ylab="y",
ylim=c(-6,4))
points(x,y)
points(x,mu,col="black",lty=1,lwd=4,type="l")
#par(mar=c(4,5,0,3))
beta_postmean <- colMeans(res$beta)
ZB_new <- sweep(ZB,MARGIN = 2,STATS = beta_postmean,"/")
#matplot(x,ZB_new,type="l",lwd=2,ylim=c(-2,2),xlab="x",ylab="")
matplot(x,ZB_new-5,type="l",lwd=2,add=TRUE)
abline(h=-4)
dev.off()
setwd("/Users/zhenkewu/Dropbox/ZW/professional/teaching/graphical\ models/notes/lecture14/Stan")
fit <- stan('model.stan',
data=c("M","C","ZB","N"),
iter=10000,warmup=5000,thin=5, chains=1)
print(fit)
res <- extract(fit,permuted=TRUE)
res_mu <- res$mu
png(file.path(result_folder,"estimated_mean_curve.png"),
width=10,height=8,units = "in",res=72)
par(mar=c(5,5,3,3))
matplot(x,t(res_mu),type="l",col=2,main=paste0("Mean curve"),
xlab="x",ylab="y",
ylim=c(-6,4))
points(x,y)
points(x,mu,col="black",lty=1,lwd=4,type="l")
#par(mar=c(4,5,0,3))
beta_postmean <- colMeans(res$beta)
ZB_new <- sweep(ZB,MARGIN = 2,STATS = beta_postmean,"/")
#matplot(x,ZB_new,type="l",lwd=2,ylim=c(-2,2),xlab="x",ylab="")
matplot(x,ZB_new-5,type="l",lwd=2,add=TRUE)
abline(h=-4)
dev.off()
setwd("/Users/zhenkewu/Dropbox/ZW/professional/teaching/graphical\ models/notes/lecture14/Stan")
fit <- stan('model.stan',
data=c("M","C","ZB","N"),
iter=10000,warmup=5000,thin=5, chains=1)
print(fit)
res <- extract(fit,permuted=TRUE)
res_mu <- res$mu
png(file.path(result_folder,"estimated_mean_curve.png"),
width=10,height=8,units = "in",res=72)
par(mar=c(5,5,3,3))
matplot(x,t(res_mu),type="l",col=2,main=paste0("Mean curve"),
xlab="x",ylab="y",
ylim=c(-6,4))
points(x,y)
points(x,mu,col="black",lty=1,lwd=4,type="l")
#par(mar=c(4,5,0,3))
beta_postmean <- colMeans(res$beta)
ZB_new <- sweep(ZB,MARGIN = 2,STATS = beta_postmean,"/")
#matplot(x,ZB_new,type="l",lwd=2,ylim=c(-2,2),xlab="x",ylab="")
matplot(x,ZB_new-5,type="l",lwd=2,add=TRUE)
abline(h=-4)
dev.off()
setwd("/Users/zhenkewu/Dropbox/ZW/professional/teaching/graphical\ models/notes/lecture14/Stan")
fit <- stan('model.stan',
data=c("M","C","ZB","N"),
iter=10000,warmup=5000,thin=5, chains=1)
print(fit)
res <- extract(fit,permuted=TRUE)
res_mu <- res$mu
png(file.path(result_folder,"estimated_mean_curve.png"),
width=10,height=8,units = "in",res=72)
par(mar=c(5,5,3,3))
matplot(x,t(res_mu),type="l",col=2,main=paste0("Mean curve"),
xlab="x",ylab="y",
ylim=c(-6,4))
points(x,y)
points(x,mu,col="black",lty=1,lwd=4,type="l")
#par(mar=c(4,5,0,3))
beta_postmean <- colMeans(res$beta)
ZB_new <- sweep(ZB,MARGIN = 2,STATS = beta_postmean,"/")
#matplot(x,ZB_new,type="l",lwd=2,ylim=c(-2,2),xlab="x",ylab="")
matplot(x,ZB_new-5,type="l",lwd=2,add=TRUE)
abline(h=-4)
dev.off()
setwd("/Users/zhenkewu/Dropbox/ZW/professional/teaching/graphical\ models/notes/lecture14/Stan")
fit <- stan('model.stan',
data=c("M","C","ZB","N"),
iter=10000,warmup=5000,thin=5, chains=1)
print(fit)
res <- extract(fit,permuted=TRUE)
res_mu <- res$mu
png(file.path(result_folder,"estimated_mean_curve.png"),
width=10,height=8,units = "in",res=72)
par(mar=c(5,5,3,3))
matplot(x,t(res_mu),type="l",col=2,main=paste0("Mean curve"),
xlab="x",ylab="y",
ylim=c(-6,4))
points(x,y)
points(x,mu,col="black",lty=1,lwd=4,type="l")
#par(mar=c(4,5,0,3))
beta_postmean <- colMeans(res$beta)
ZB_new <- sweep(ZB,MARGIN = 2,STATS = beta_postmean,"/")
#matplot(x,ZB_new,type="l",lwd=2,ylim=c(-2,2),xlab="x",ylab="")
matplot(x,ZB_new-5,type="l",lwd=2,add=TRUE)
abline(h=-4)
dev.off()
setwd("/Users/zhenkewu/Dropbox/ZW/professional/teaching/graphical\ models/notes/lecture14/Stan")
fit <- stan('model.stan',
data=c("M","C","ZB","N"),
iter=10000,warmup=5000,thin=5, chains=1)
print(fit)
res <- extract(fit,permuted=TRUE)
res_mu <- res$mu
png(file.path(result_folder,"estimated_mean_curve.png"),
width=10,height=8,units = "in",res=72)
par(mar=c(5,5,3,3))
matplot(x,t(res_mu),type="l",col=2,main=paste0("Mean curve"),
xlab="x",ylab="y",
ylim=c(-6,4))
points(x,y)
points(x,mu,col="black",lty=1,lwd=4,type="l")
#par(mar=c(4,5,0,3))
beta_postmean <- colMeans(res$beta)
ZB_new <- sweep(ZB,MARGIN = 2,STATS = beta_postmean,"/")
#matplot(x,ZB_new,type="l",lwd=2,ylim=c(-2,2),xlab="x",ylab="")
matplot(x,ZB_new-5,type="l",lwd=2,add=TRUE)
abline(h=-4)
dev.off()
library(shinystan)
my_sso <- launch_shinystan(fit)
rm(list=ls())
library(baker)
# set up parameter values:
N       <- 500 # no. of cases/controls.
seed_start <- 20150923
set.seed(seed_start)
p <- 1
X0 <- matrix(runif(N*p,0,1),nrow=N,ncol=p)
X  <- apply(X0,2,scale)
f1 <- function(x) {sin(8*pi*(x-0.5)/7)}
f2 <- function(x) {-0.5*sin(2*pi*x)/(pi*x)}
f3 <- function(x) {-0.5*sin(8*pi*(x-0.5)/7)}
f4 <- function(x) {0}
f5 <- function(x) {4*(exp(3*x)/(1+exp(3*x))-0.5)}
f6 <- function(x) {0}
f  <- f1
x  <- sort(X[,1])
mu <- f(x)
sigma <- 1
y  <- mu + rnorm(N,0,sigma)
plot(x,y)
points(x,mu,col="blue",lwd=2)
# include stratification information in file name:
dated_strat_name    <- file.path("/Users/zhenkewu/Dropbox/ZW/professional/teaching/graphical\ models/notes/lecture14/code",paste0("JAGS"))
# create folder
dir.create(dated_strat_name)
fullname <- dated_strat_name
# for finer scenarios, e.g., different types of analysis applicable to the
# same data set. Here we just perform one analysis:
result_folder <- fullname
dir.create(result_folder)
# prepare data:
M <- y
C0 <- 10
C  <- C0
myknots2 <- quantile(x,seq(0,1,length = (C0 - 2))[-c(1,(C0 - 2))])
ZB       <- matrix(splines::bs(x,knots= myknots2,intercept=TRUE),nrow=N)
# #par(mfrow=c(2,1))
# matplot(x,ZB)
# abline(v=myknots)
# matplot(x,ZB2,add=TRUE)
# abline(v=myknots2,col="red")
matplot(x,ZB,main="B-spline basis used here.",type="l",lwd=2)
diff_ord <- 1
D  <- diff(diag(C),differences = diff_ord)
PM <- t(D)%*%D
# in_data:
in_data <- c("M","C","ZB","N")
# out_parameter:
out_parameter <- c("beta","mu","taubeta","sigma")
# in_init:
in_init <- function(){
list(beta  = rep(0,C))
}
# fit Bayesian model:
mcmc_options <- list(debugstatus     = TRUE,
n.chains   = 1,
n.itermcmc = 10000,
n.burnin   = 5000,
n.thin     = 5,
result.folder = result_folder,
bugsmodel.dir = result_folder,
winbugs.dir   = "C:\\Program Files\\WinBUGS14\\",
jags.dir      = ""#"C:\\Program Files\\JAGS\\JAGS-3.4.0\\x64\\bin\\jags"
)
filename <- file.path(result_folder,"model.bug")
# gs <- R2WinBUGS::bugs(data     = in_data,
#                       inits    = in_init,
#                       parameters.to.save = out_parameter,
#                       model.file = filename,
#                       working.directory=mcmc_options$result.folder,
#                       bugs.directory  = mcmc_options$winbugs.dir,  #<- special to WinBUGS.
#                       n.iter         = mcmc_options$n.itermcmc,
#                       n.burnin       = mcmc_options$n.burnin,
#                       n.thin         = mcmc_options$n.thin,
#                       n.chains       = mcmc_options$n.chains,
#                       DIC      = FALSE,
#                       debug    = mcmc_options$debugstatus)
lapply(in_data, dump, append = TRUE, file = file.path(mcmc_options$result.folder,"jagsdata.txt"))
library(rjags)
load.module("glm")
time1 <- Sys.time()
# jags require more meticulous initial values specifications:
gs <- R2jags::jags2(data   = in_data,
inits  = in_init,
parameters.to.save = out_parameter,
model.file = filename,
working.directory = mcmc_options$result.folder,
n.iter         = as.integer(mcmc_options$n.itermcmc),
n.burnin       = as.integer(mcmc_options$n.burnin),
n.thin         = as.integer(mcmc_options$n.thin),
n.chains       = as.integer(mcmc_options$n.chains),
DIC            = FALSE,
clearWD        = FALSE#,              #<--- special to JAGS.
);
time2 <- Sys.time()
time2-time1
#
#
#
#
#
#
#
#
#
# Visualize Posterior Results:
#
#
#
#
#
#
#
#
#
#
#
# # WinBUGS:
# DIR_NPLCM <- result_folder
# bugs.dat  <- dget(file.path(DIR_NPLCM,"data.txt"))
# res_nplcm <- coda::read.coda(file.path(DIR_NPLCM,"coda1.txt"),
#                              file.path(DIR_NPLCM,"codaIndex.txt"),
#                              quiet=TRUE)
# JAGS:
DIR_NPLCM <- result_folder
new_env <- new.env()
source(file.path(DIR_NPLCM,"jagsdata.txt"),local=new_env)
bugs.dat <- as.list(new_env)
rm(new_env)
res <- coda::read.coda(file.path(DIR_NPLCM,"CODAchain1.txt"),
file.path(DIR_NPLCM,"CODAindex.txt"),
quiet=TRUE)
print_res <- function(x) plot(res[,grep(x,colnames(res))])
get_res   <- function(x) res[,grep(x,colnames(res))]
# get specific results:
ind_mu <- grep("^mu",colnames(res))
res_mu <- res[,ind_mu]
#
# plot linear predictors:
#
png(file.path(result_folder,"estimated_mean_curve.png"),
width=10,height=8,units = "in",res=72)
#layout(matrix(1:2,nrow=2,ncol=1),width=10,height=c(6,2))
par(mar=c(5,5,3,3))
matplot(x,t(res_mu),type="l",col=2,main=paste0("Mean curve"),
xlab="x",ylab="y",
ylim=c(-6,4))
points(x,y)
points(x,mu,col="black",lty=1,lwd=4,type="l")
#par(mar=c(4,5,0,3))
beta_postmean <- colMeans(get_res("beta\\["))
ZB_new <- sweep(ZB,MARGIN = 2,STATS = beta_postmean,"/")
#matplot(x,ZB_new,type="l",lwd=2,ylim=c(-2,2),xlab="x",ylab="")
matplot(x,ZB_new-5,type="l",lwd=2,add=TRUE)
abline(h=-4)
dev.off()
rm(list=ls())
library('rstan')
library('Rcpp')
# set up parameter values:
N       <- 500 # no. of cases/controls.
seed_start <- 20150923
set.seed(seed_start)
p <- 1
X0 <- matrix(runif(N*p,0,1),nrow=N,ncol=p)
X  <- apply(X0,2,scale)
f1 <- function(x) {sin(8*pi*(x-0.5)/7)}
f2 <- function(x) {-0.5*sin(2*pi*x)/(pi*x)}
f3 <- function(x) {-0.5*sin(8*pi*(x-0.5)/7)}
f4 <- function(x) {0}
f5 <- function(x) {4*(exp(3*x)/(1+exp(3*x))-0.5)}
f6 <- function(x) {0}
f  <- f1
x  <- sort(X[,1])
mu <- f(x)
sigma <- 1
y  <- mu + rnorm(N,0,sigma)
plot(x,y)
points(x,mu,col="blue",lwd=2)
# include stratification information in file name:
dated_strat_name    <- file.path(paste0("Stan"))
# create folder
dir.create(dated_strat_name)
fullname <- dated_strat_name
# for finer scenarios, e.g., different types of analysis applicable to the
# same data set. Here we just perform one analysis:
result_folder <- fullname
dir.create(result_folder)
# prepare data:
M <- y
C0 <- 10
C  <- C0
myknots2 <- quantile(x,seq(0,1,length = (C0 - 2))[-c(1,(C0 - 2))])
ZB       <- matrix(splines::bs(x,knots= myknots2,intercept=TRUE),nrow=N)
setwd("/Users/zhenkewu/Dropbox/ZW/professional/teaching/graphical\ models/notes/lecture14/code/Stan")
fit <- stan('model.stan',
data=c("M","C","ZB","N"),
iter=10000,warmup=5000,thin=5, chains=1)
print(fit)
res <- extract(fit,permuted=TRUE)
res_mu <- res$mu
png(file.path(result_folder,"estimated_mean_curve.png"),
width=10,height=8,units = "in",res=72)
par(mar=c(5,5,3,3))
matplot(x,t(res_mu),type="l",col=2,main=paste0("Mean curve"),
xlab="x",ylab="y",
ylim=c(-6,4))
points(x,y)
points(x,mu,col="black",lty=1,lwd=4,type="l")
#par(mar=c(4,5,0,3))
beta_postmean <- colMeans(res$beta)
ZB_new <- sweep(ZB,MARGIN = 2,STATS = beta_postmean,"/")
#matplot(x,ZB_new,type="l",lwd=2,ylim=c(-2,2),xlab="x",ylab="")
matplot(x,ZB_new-5,type="l",lwd=2,add=TRUE)
abline(h=-4)
dev.off()
# to explore the posterior:
library(shinystan)
my_sso <- launch_shinystan(fit)
