########## R script: mitsubBayesStan ########## # For Bayesian simple linear regression analysis # of price of Mitsubishi cars versus their age. # This version is done using Stan. # Last changed: 02 SEP 2022 # Load required package and functions: library(rstan) source("summMCMC.r") # Set flag for code re-compilation: reCompileCode <- TRUE # Set MCMC parameters: nWarm <- 1000 # Length of warm-up. nIter <- 1000 # Size of the kept sample. nThin <- 1 # Thinning factor. # Put data into arrays: age <- c(13,14,14,12,9,15,10,14,9,14,13,12,9,10,15, 11,15,11,7,13,13,10,9,6,11,15,13,10,9,9,15, 14,14,10,14,11,13,14,10) price <- c(2950,2300,3900,2800,5000,2999,3950,2995,4500, 2800,1990,3500,5100,3900,2900,4950,2000,3400, 8999,4000,2950,3250,3950,4600,4500,1600,3900, 4200,6500,3500,2999,2600,3250,2500,2400,3990, 4600,450,4700) # Plot the data: plot(age,price,bty="l",col="dodgerblue") readline("Hit Enter to continue.\n") # Standardise the continuous predictor birweight: meanAge <- mean(age) sdAge <- sd(age) x <- (age-meanAge)/sdAge meanPrice <- mean(price) sdPrice <- sd(price) y <- (price-meanPrice)/sdPrice # Specify model in Stan: slrModel <- 'data { int n; vector[n] y; matrix[n,2] X; real sigmabeta; real A; real B; } parameters { vector[2] beta; real sigsq; } model { y ~ normal(X*beta,sqrt(sigsq)); beta ~ normal(0,sigmabeta); sigsq ~ inv_gamma(A,B); }' # Store data in a list in format required by Stan: X <- cbind(1,x) allData <- list(n=39,y=y,X=X,sigmabeta=1e5,A=0.01,B=0.01) # Compile code for model if required: if (reCompileCode) { stanCompilObj <- stan(model_code=slrModel,data=allData, iter=1,chains=1) save(stanCompilObj,file="slrModel.rda") } if (!reCompileCode) load("slrModel.rda") # Obtain MCMC samples for each parameter using Stan: stanObj <- stan(model_code=slrModel,data=allData, warmup=nWarm,iter=(nWarm+nIter), chains=1,thin=nThin,refresh=100,fit=stanCompilObj) betaMCMC <- NULL for (j in 1:2) { charVar <- paste("beta[",as.character(j),"]",sep="") betaMCMC <- rbind(betaMCMC,extract(stanObj,charVar,permuted=FALSE)) } beta0MCMC <- betaMCMC[1,] beta1MCMC <- betaMCMC[2,] sigsqMCMC <- as.vector(extract(stanObj,"sigsq",permuted=FALSE)) # Transform to original units: beta1OrigMCMC <- sdPrice*beta1MCMC/sdAge beta0OrigMCMC <- meanPrice + sdPrice*(beta0MCMC - beta1MCMC*meanAge/sdAge) sigsqOrigMCMC <- (sdPrice^2)*sigsqMCMC # Summarise the results: summMCMC(list(cbind(beta0OrigMCMC,beta1OrigMCMC,sigsqOrigMCMC)), parNames=list(c(expression(beta[0])),c(expression(beta[1])), c(expression(sigma^2)))) readline("Hit Enter to continue.\n") # Plot fit to the data: plot(age,price,bty="l",type="n") ng <- 1001 xg <- seq(min(age),max(age),length=ng) Xg <- cbind(rep(1,ng),xg) mugMCMC <- Xg%*%rbind(beta0OrigMCMC,beta1OrigMCMC) mug <- apply(mugMCMC,1,mean) credLowerg <- apply(mugMCMC,1,quantile,0.025) credUpperg <- apply(mugMCMC,1,quantile,0.975) polygon(c(xg,rev(xg)),c(credLowerg,rev(credUpperg)), col="palegreen",border=FALSE) lines(xg,mug,lty=1,lwd=2,col="darkgreen") points(age,price,col="dodgerblue") ########## End of mitsubBayesStan ##########