I am fitting a GEV distribution with a time trend in the location parameter \mu, such that
where t_c = \lfloor N / 2 \rfloor, i.e. the index of the observation in the middle of the data set. Shown below is the model implementation in code:
functions {
real gevt_lpdf(vector y, real mu_0, real sigma, real xi, real delta) {
int n = rows(y);
int c = n / 2;
vector[n] mu;
for (i in 1:n) { mu[i] = mu_0 * (1 + delta * (i - c)); }
vector[n] eta = (abs(xi) < 1e-12)
? exp(-(y - mu) / sigma)
: fmax(1e-12, pow(1 + xi * (y - mu) / sigma, -1 / xi));
real lp = -n * log(sigma) + sum((1 + xi) * log(eta)) - sum(eta);
return lp;
}
}
data {
int<lower = 0> n;
vector[n] y;
}
parameters {
real mu_0;
real<lower = 0> sigma;
real<lower = 0, upper = 1> xi_0;
real delta;
}
transformed parameters {
real xi = xi_0 - 0.5;
}
model {
mu_0 ~ normal(0, 10);
sigma ~ exponential(3);
xi_0 ~ beta(4, 4);
delta ~ normal(0, 10);
y ~ gevt(mu_0, sigma, xi, delta);
}
generated quantities {
vector[n] log_lik;
for (i in 1:n)
log_lik[i] = gevt_lpdf([y[i]]' | mu_0, sigma, xi, delta);
}
Model convergence is good, \widehat R \approx 1 and the effective sample size is larger than the number of draws. But I am having this confusing issue where the loo
functions psis
and loo
seemingly disagree on the Pareto k values.
This is what happens when I run psis
on the samples from the above model (timetrend_fit
):
> timetrend_loglik <- loo::extract_log_lik(timetrend_fit)
> loo::psis(timetrend_loglik)
Computed from 40000 by 132 log-weights matrix.
MCSE and ESS estimates assume independent draws (r_eff=1).
All Pareto k estimates are good (k < 0.7).
See help('pareto-k-diagnostic') for details.
Whereas for the loo
function this is the produced output:
> loo::loo(timetrend_loglik)
Computed from 40000 by 132 log-likelihood matrix.
Estimate SE
elpd_loo -494.3 20.7
p_loo 15.6 13.1
looic 988.6 41.4
------
MCSE of elpd_loo is NA.
MCSE and ESS estimates assume independent draws (r_eff=1).
Pareto k diagnostic values:
Count Pct. Min. ESS
(-Inf, 0.7] (good) 131 99.2% 23774
(0.7, 1] (bad) 0 0.0% <NA>
(1, Inf) (very bad) 1 0.8% <NA>
See help('pareto-k-diagnostic') for details.
Warning message:
Some Pareto k diagnostic values are too high. See help('pareto-k-diagnostic') for details.
There is one problematic observation w.r.t. to the Pareto k value. Is this a bug in the implementation, or is there a difference in how psis
and loo
estimate the Pareto k diagnostic values? What can I do to fix this?