Adding quantile noise example.

m5-estimate-hierarchy.stan

@@ -22,7 +22,8 @@ model {
   for (i in 1:(M*N)) { 
     real mu_local;
     mu_local <- mu + team_effects[team[i]];
-    increment_log_prob(normal_log(y[i], mu_local, sigma));
+    increment_log_prob(normal_log(y[i], 
+      mu_local, sigma));
   }
 }
 

m6-estimate-nc-hierarchy.stan

@@ -20,15 +20,17 @@ parameters {
 transformed parameters {
   real batch_mean[M];
   for (m in 1:M) {
-    batch_mean[m] <- mu + team_effects_raw[batch_team_index[m]]*team_sigma;
+    batch_mean[m] <- mu + 
+      team_effects_raw[batch_team_index[m]] * 
+      team_sigma;
   }
 }
 
 model {
-  mu ~ normal(0,5);
+  mu ~ normal(0,2);
   team_effects_raw ~ normal(0,1);
   team_sigma ~ gamma(2,1);
-  sigma ~ gamma(2,.1);
+  sigma ~ gamma(2,1);
   for (i in 1:(M*N)) { 
     real mu_local;
     mu_local <- batch_mean[widget_batch_index[i]];

m7-estimate-nc-t-hierarchy.stan

@@ -13,14 +13,16 @@ data {
 parameters {
   real mu;
   real<lower=0> sigma;
-  real<lower=-pi()/2, upper=pi()/2> team_effects_raw[K];
+  real<lower=-pi()/2, upper=pi()/2> 
+    team_effects_raw[K];
   real<lower=0> team_sigma;
 }
 
 transformed parameters {
   real batch_mean[M];
   for (m in 1:M) {
-    batch_mean[m] <- mu + tan(team_effects_raw[batch_team_index[m]])*team_sigma;
+    batch_mean[m] <- mu + 
+      tan(team_effects_raw[batch_team_index[m]])*team_sigma;
   }
 }
 

odsc-10-quantile-estimation.R

@@ -0,0 +1,35 @@
+## Quantiles are often used to calculate estimates from
+## Bayesian samples and they are often easily applied.
+## Quantiles also have a few properties that make them
+## a little touchy to use.  This script addresses one
+## particular issue: commonly used 2.5%/97.5% quantile
+## estimates are very noise when estimated from samples
+## making it critical that end-users understand the estimate
+## they are producing.  In general the tail weight of the
+## distribution giving rise to the samples is critical 
+## in determining the noise level, but here we do a 
+## simulation study for the normal distribution and 
+## the usual 2.5%/97.5% quantiles.
+
+library(rstan); library(shinystan)  # 
+library(dplyr); library(tidyr)
+
+# We are simulating M samples {m:1<=M}, with N=m^2 draws
+# per sample from a standard normal.
+M <- 10^3
+mu <- 0
+sigma <- 1
+
+s15 <- list()
+for (m in 1:M) {
+  s15[[m]] <- rnorm(m^2,0,1)
+}
+
+e15 <- sapply(s15, quantile, prob=0.975)
+
+plot( x=(1:M)^2, y=e15/qnorm(0.975), ylim=c(0.9,1.1), pch='.');
+abline(h=1); abline(v=10^3, col='red'); abline(v=10^4, col='red')
+
+
+
+