Gamma Exponential with Stan

Gamma_Exponential_Stan.R

library(rstan)
stanmodelcode <- "
data {
  int<lower=0> N;
  int<lower=0> y[N];
}
parameters {
  real<lower=0.00001> Theta;
}
model {
  Theta ~ gamma(4, 1000);
  for (n in 1:N)
    y[n] ~ exponential(Theta);
}
generated quantities{
  real y_pred;
  y_pred = exponential_rng(Theta);
}
"
stanDso <- stan_model(model_code = stanmodelcode) 
claims.obs <- c(100, 950, 450)
N <- length(claims.obs)
dat <- list(N = N, y = claims.obs); 

fit <- sampling(stanDso, data = dat, iter = 10000, warmup=200) 
fit
## Inference for Stan model: stanmodelcode.
## 4 chains, each with iter=10000; warmup=200; thin=1; 
## post-warmup draws per chain=9800, total post-warmup draws=39200.
## 
##          mean se_mean     sd   2.5%    25%    50%    75%   97.5% n_eff Rhat
## Theta    0.00    0.00   0.00   0.00   0.00   0.00   0.00    0.01 13756    1
## y_pred 416.86    2.85 492.28   8.83 106.34 262.70 542.47 1730.55 29773    1
## lp__   -48.65    0.01   0.70 -50.68 -48.83 -48.38 -48.20  -48.15 14203    1
## 
## Samples were drawn using NUTS(diag_e) at Tue May 19 06:06:08 2015.
## 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).

Theta <- extract(fit, 'Theta')
Theta <- unlist(Theta, use.names=FALSE)
y_pred <- extract(fit, 'y_pred')
y_pred <- unlist(y_pred, use.names=FALSE)

## Analytical hyper-parameters
prior.alpha <- 4
prior.beta <- 1000
posterior.alpha <- prior.alpha + N
posterior.beta <- prior.beta + sum(claims.obs)

op <- par(mfrow=c(1,2))
# Simulated posterior parameter
plot(density(Theta),
     xlab=expression(Theta), col=grey(0, 0.8),
     main="Parameter distribution")
# Analytical posterior parameter
curve(dgamma(x, posterior.alpha, posterior.beta), 
      add=TRUE, col=4, lty=2, lwd=1.5)
# Analytical prior parameter
curve(dgamma(x, prior.alpha, prior.beta), 
      add=TRUE, col=2)
legend(x="topright", col=c(2, 4, grey(0, 0.8)), lty=c(1,2,1), bty="n",
       legend=c("Prior", "Analytical posterior", "Sampling posterior"))

# Simulated posterior predictive
plot(density(y_pred), xlim=c(1,2000),
     xlab="Loss", col=grey(0, 0.8),
     main="Predicitive distribution")
# Analytical posterior predictive
library(actuar) # Required for pareto distribution
curve(dpareto(x, posterior.alpha, posterior.beta), 
      add=TRUE, col=4, lwd=1.5, lty=2)
# Analytical prior predictive
curve(dpareto(x, prior.alpha, prior.beta), 
      add=TRUE, col=2)
legend(x="topright", col=c(2, 4, grey(0, 0.8)), lty=c(1,2,1), bty="n",
       legend=c("Prior","Analytical posterior", "Sampling posterior"))
par(op)

# Review quantiles and probabilities
qpareto(0.75, prior.alpha, prior.beta)
qpareto(0.75, posterior.alpha, posterior.beta)
quantile(y_pred, 0.75)
ppareto(950, posterior.alpha, posterior.beta)
ecdf(y_pred)(950)