Skip to content

Instantly share code, notes, and snippets.

@mages mages/NormalDistribution.R
Last active Sep 26, 2016

Embed
What would you like to do?
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")
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.