Hi,

I am trying to use pseudo-BMA as described in “Using Stacking to Average Bayesian Predictive Distributions” by Yao et al 2018.

Following R code is directly from “loo_model_weights” function help page

```
library(rstan)
# generate fake data from N(0,1).
N <- 100
y <- rnorm(N, 0, 1)
# Suppose we have three models: N(-1, sigma), N(0.5, sigma) and N(0.6,sigma).
stan_code <- "
data {
int N;
vector[N] y;
real mu_fixed;
}
parameters {
real<lower=0> sigma;
}
model {
sigma ~ exponential(1);
y ~ normal(mu_fixed, sigma);
}
generated quantities {
vector[N] log_lik;
for (n in 1:N) log_lik[n] = normal_lpdf(y[n]| mu_fixed, sigma);
}"
mod <- stan_model(model_code = stan_code)
fit1 <- sampling(mod, data=list(N=N, y=y, mu_fixed=-1))
fit2 <- sampling(mod, data=list(N=N, y=y, mu_fixed=0.5))
fit3 <- sampling(mod, data=list(N=N, y=y, mu_fixed=0.6))
model_list <- list(fit1, fit2, fit3)
log_lik_list <- lapply(model_list, extract_log_lik)
# optional but recommended
r_eff_list <- lapply(model_list, function(x) {
ll_array <- extract_log_lik(x, merge_chains = FALSE)
relative_eff(exp(ll_array))
})
# stacking method:
wts1 <- loo_model_weights(
log_lik_list,
method = "stacking",
r_eff_list = r_eff_list,
optim_control = list(reltol=1e-10)
)
print(wts1)
```

From this, I obtain the weights for each model.

My question is: How can I calculate the model averaged posterior estimates, in this example “model averaged” estimate of the sigma parameter?

(Although this entry (Implementing Model Averaging) looks similar, unfortunately still I did not understand it.)

Thank you very much for your help,

Best regards,

Burak Kürsad Günhan