library(rstan)
library(MASS)
set.seed(1)
y <- rnorm(100, 4, 2)
truehist(y, col="#B2001D")
lines(density(y), col="skyblue", lwd=2)
summary(y)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## -0.4294  3.0120  4.2280  4.2180  5.3830  8.8030

ret <- stanc(file="NormaLDistribution.stan") # Check Stan file
ret_sm <- stan_model(stanc_ret = ret) # Compile Stan code
fit <- sampling(ret_sm, warmup=100, iter=600, seed=1,
                data=list(y, N=length(y)))
stan_trace(fit, inc_warmup = TRUE)
stan_hist(fit)
print(fit, probs=c(0.025, 0.5, 0.975))
## Inference for Stan model: NormaLDistribution.
## 4 chains, each with iter=600; warmup=100; thin=1; 
## post-warmup draws per chain=500, total post-warmup draws=2000.
## 
##          mean se_mean   sd    2.5%     50%   97.5% n_eff Rhat
## mu       4.22    0.01 0.17    3.88    4.22    4.53   771    1
## sigma    1.82    0.00 0.13    1.57    1.81    2.10  2000    1
## y_ppc    4.22    0.04 1.86    0.70    4.21    7.88  2000    1
## lp__  -200.32    0.03 0.97 -203.02 -200.02 -199.41  1016    1
## 
## Samples were drawn using NUTS(diag_e) at Mon Sep 26 07:25:09 2016.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at 
## convergence, Rhat=1).
summary(extract(fit, "y_ppc")[["y_ppc"]])
##   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## -1.483   2.959   4.213   4.222   5.451  10.540 
plot(ecdf(y), main="Posterior predictive check")
lines(ecdf(extract(fit, "y_ppc")[["y_ppc"]]), col="#B2001D")