library(rstan)
library(shinystan)

model_file = 'C:/User/balajir/Documents/Research/Bayes-Modeling-R/example_models/examples/4-multiple-linear-regression/mlr.stan'
    # sampler iterations
iterations = 5000
warmup = 2000
rng_seed = 1111

df = read.csv("C:/User/balajir/Documents/Research/Bayes-Modeling-R/example_models/examples/4-multiple-linear regression/colorado_tmin_ave.csv")
predictors = as.matrix(df[,c('lon','lat','elev')])
response = df$tmin
N = nrow(df)
K = ncol(predictors)

fit = lm(tmin~lon+lat+elev,data=df)

initf = function()list(
    alpha = coefficients(fit)[1],
    betas = coefficients(fit)[-1],
    sigma = sd(residuals(fit)))

stan_data = list(y=response, x=predictors, N=N, K=K)

# setup the model
model = stan(model_file, data = stan_data, chains = 0)
stanfit = stan(fit = model, data = stan_data, 
                chains = 1, iter=iterations, 
                seed=rng_seed, init=initf, 
                warmup = warmup)

### get the output stats - mean, median standard deviation from the posterior
### of the paramters
print(stanfit,digits=2)

# look at the results - in terms of convergence
traceplot(stanfit)

### Plot the posterior
dev.new()
par(mfcol=c(1,2))
hist(extract(stanfit)$mu)
hist(extract(stanfit)$sigma)


### use shinystan
launch_shinystan(stanfit)