hello
E_loo uses the weights generated by loo to estimate the leave-one-out mean (and quantiles). However, “A small Pareto k estimate is necessary, but not sufficient, for E_loo() to give reliable estimates.” Mixis is more robust for estimating elpd. Can the log-weights generated by mixis be used to improve the estimate of the leave-one-out mean?
Below is a simple gaussian example with synthetic data including one clear outlier. Putting eloov as the vector of mean values obtained from E_loo, I compare the true (brute force) expected value at the outlier, eloov at the outlier, and a proposed way to use mixis to give an another estimate, which is better at the outlier in this example. Is this a valid method for this purpose? Is there a better way? Does it perform better than eloov in general?
thanks
Greg
R
library(brms)
library(rstan)
rstan_options(auto_write = TRUE)
library(brms)
library(cmdstanr)
options(mc.cores=2, brms.normalize=TRUE,brms.backend="cmdstanr")
library(ggplot2)
library(loo)
library("matrixStats")
set.seed(1)
N=400
x<-runif(N)
f<-function(t) $1.8+sin(12*t)*exp(-3*t)$
y<-$f(x)+0.3*rnorm(N)$
y[1]=-1
dat<-data.frame(x,y)
bf1<-bf(y~s(x))
mod1<-brm(bf1,data=dat,family="gaussian",adapt_delta=0.9)
LL<-log_lik(mod1)
eLL<-exp(LL)
cid<-c(rep(1,1000),rep(2,1000),rep(3,1000),rep(4,1000))
reff<-relative_eff(eLL,chain_id=cid)
loo1<-loo(LL,r_eff=reff,save_psis=TRUE)
loo1
loopt1<-loo1$pointwise
head(loopt1)
pospred<-posterior_epred(mod1)
psis_object<-loo1$psis_object
log_ratios<- -1*LL
eloo<-E_loo(pospred,psis_object,type="mean",log_ratios=log_ratios)
eloov<-eloo$value
#brute force estimate
datx1<-dat[-1,]
modx1<-brm(bf1,data=datx1,family="gaussian",adapt_delta=0.9)
newd<-data.frame(x=x[1])
ppd<-posterior_predict(modx1,newdata=newd)
mean(ppd)
#1.787043
eloov[1]
#1.795551
#mixis
cat(mod1$model,file="mod1.stan")
#save as mod1mix.stan
#amend the model block with mixis code lines
#amend the generated quantities block to generate mu and log_lik
#save as mod1mix.stan
the two amended blocks now read
model {
// likelihood including constants
if (!prior_only) {
// initialize linear predictor term
vector[N] log_lik = rep_vector(0.0, N);
vector[N] mu = rep_vector(0.0, N);
mu += Intercept + Xs * bs + Zs_1_1 * s_1_1;
for (n in 1:N) {
log_lik[n] = normal_lpdf(Y[n] | mu[n], sigma);
}
target += sum(log_lik);
target += log_sum_exp(-log_lik);
}
// priors including constants
target += lprior;
target += std_normal_lpdf(zs_1_1);
}
generated quantities {
// actual population-level intercept
real b_Intercept = Intercept;
vector[N] log_lik = rep_vector(0.0, N);
vector[N] mu = rep_vector(0.0, N);
mu += Intercept + Xs * bs + Zs_1_1 * s_1_1;
for (n in 1:N) {
log_lik[n] = normal_lpdf(Y[n] | mu[n], sigma);
}
#make Cmdstanr model and run it
mod1.mix<-cmdstan_model("mod1mix.stan")
class(mod1.mix)
mod1.mix$print()
stand<-standata(mod1)
mix1<-mod1.mix$sample(data=stand,chains=4,parallel_chains=2,refresh=100,
iter_warmup=1000,iter_sampling=1000,save_warmup=TRUE,max_treedepth=13,adapt_delta=0.9)
mix1$save_object(file="mix1.RDS")
#mix1<-readRDS("mix1.RDS")
#elpd from mixis
draws<-mix1$draws()
drawmat<-as_draws_matrix(draws)
dim(drawmat)
grep("log_lik",colnames(drawmat))[c(1,N)]
dwLL<-drawmat[,24:423]
dim(dwLL)
colnames(dwLL)[c(1,N)]
l_common_mix <- rowLogSumExps(-dwLL)
log_weights <- -dwLL - l_common_mix
elpd_mixis <- logSumExp(-l_common_mix) - rowLogSumExps(t(log_weights))
save(elpd_mixis,file="elpd_mixis.rds")
#compare elpd_loo
elpd_loo<-loopt1[,1]
plot(elpd_mixis,elpd_loo)
#proposed method for E_mix
lw<-log_weights
dim(lw)
#make normalised weights
lse<-apply(lw,2,logSumExp)
lseo<-outer(rep(1,4000),lse)
lwn<-lw-lseo
w<-exp(lwn)
dim(w)
range(w)
apply(w,2,sum)[1:10]
dwmu<-drawmat[,424:823]
dim(dwmu)
colnames(dwmu)[c(1,N)]
emix<-apply(dwmu*w,2,sum)
plot(emix,eloov)
mean(ppd)
#1.787043
eloov[1]
#1.795551
emix[1]
#1.783195