Thank you again Bob. I tried (my (mis)interpretation of) your suggestion to estimate the normal mu and sigma based on sufficient statistics. Bottom line - I am still seeing a difference in the estimate of sigma relative to that obtained based on the raw data. Likely I have incorrectly expressed the corresponding reference prior when using sufficient statistics. My worry is that it may not be possible, when estimating from sufficient statistics, to identify a reference prior for sigma that corresponds to the usual Jeffreys reference prior used when estimating from the raw data.
Below I provide my R and Stan codes (separate models for each method) as well as the estimates (you will see there is a difference in sigma posterior mean). Any thoughts greatly appreciated!!
Here is the Stan code used for estimation based on the full raw data vector:
// estimation of mu and sigma based on Y
data {
int<lower=0> nY;
vector[nY] Y;
}
parameters {
real mu; // implied reference prior
real log_sigma; // implied reference prior
}
transformed parameters {
real<lower=0> sigma = exp(log_sigma);
}
model {
Y ~ normal(mu,sigma);
}
Here is the Stan code used for estimation based on sufficient statistics (note the transformed parameter statement has been commented out)
// Estimation of mu and sigma based on Ybar and S2
data {
int<lower=0> nY;
real Ybar;
real<lower=0> S2;
}
parameters {
real mu;
real<lower=0> sigma;
}
transformed parameters {
// real<lower=0> sigma = exp(log_sigma);
}
model {
Ybar ~ normal(mu, sigma / sqrt(nY));
S2 ~ gamma((nY - 1) / 2, (nY - 1) / 2 / sigma^2);
target += -2 * log(sigma);
}
Below is the R code used to execute both of the above models
require(rstan) # Needed for Bayesian nonlinear model
options(mc.cores = parallel::detectCores()) # recommended by Stan
rstan_options(auto_write = TRUE) # avoids recompilation
# Raw data
Y <- c(13.98251, 14.62196, 13.88309, 13.77684, 14.70978,
14.11597, 14.70279, 13.83011, 14.61396, 15.32342)
###### Estimation of mu and sigma based on raw data
###### Assume Jeffreys reference priors
# Define data for Stan
dataList = list(
Y=Y, nY =length(Y)
)
# Define the parameters for Stan
parameters = c("mu","sigma")
# Translate model to C++ and compile to DSO
StabDSO <- stan_model( 'Bob C example.stan' )
# Obtain sample of n.chains*(iter - warmup)/thin posterior draws
n.chains <- 4
options (mc.cores = n.chains)
stanFit <- sampling( object=StabDSO ,
pars=parameters,
data = dataList ,
chains = n.chains ,
iter = 100000 ,
warmup = 20000 ,
thin = 1,verbose=TRUE)
stan.sum<-summary(stanFit,probs=c(0.025,0.5,0.975))
round(stan.sum$summary,digits=3)
###### Estimation of mu and sigma based on sufficient statistics
###### Assume Jeffreys reference priors
# Define data for Stan
dataList = list(
nY =length(Y), Ybar=mean(Y), S2=var(Y)
)
# Define the parameters for Stan
parameters = c("mu","sigma")
# Translate model to C++ and compile to DSO
StabDSO <- stan_model( 'Bob C example SS.stan' )
# Obtain sample of n.chains*(iter - warmup)/thin posterior draws
n.chains <- 4
options (mc.cores = n.chains)
stanFit <- sampling( object=StabDSO ,
pars=parameters,
data = dataList ,
chains = n.chains ,
iter = 100000 ,
warmup = 20000 ,
thin = 1,verbose=TRUE)
stan.sum<-summary(stanFit,probs=c(0.025,0.5,0.975))
round(stan.sum$summary,digits=3)
Below are the estimates obtained based on the raw data vector Y (note these estimates agree exactly with the GPQ algorithm I am trying to duplicate in Stan. This GPQ algorithm however is based on sampling distributions of the sufficient statistics not the Y vector)
mean se_mean sd 2.5% 50% 97.5% n_eff Rhat
mu 14.356 0.000 0.183 13.992 14.356 14.719 154587.6 1
sigma **0.559** 0.000 0.152 0.351 0.531 0.935 135685.9 1
lp__ 1.135 0.003 1.101 -1.842 1.472 2.210 100395.1 1
Below are the estimates based on the sufficient statistics Ybar and S2. Notice that the posterior mean of sigma (0.527) is lower than the estimate based on the Y vector (0.559). I expect these estimates to be identical.
mean se_mean sd 2.5% 50% 97.5% n_eff Rhat
mu 14.356 0.000 0.173 14.012 14.356 14.702 159543.8 1
sigma **0.527** 0.000 0.143 0.331 0.500 0.881 136521.0 1
lp__ 8.422 0.003 1.105 5.462 8.761 9.496 103504.7 1
