I want to fit a multi-level generalised pareto distribution (GPD) to negated data. I have used stan code from here. While running the model, some chains are not initializing at all but other few are running fine and converging properly. I get the following error for those chains which are not initializing

More specifically, in my study, I have 9 subgroups, for which TG depends upon SG. As single threshold cannot classify extreme in all groups, I am using 9 different thresholds obtained from threshold stability plots. As I am interested in minimum (Prob (TG<=0)), I have used negated value of TG to fit a GPD. Therefore, I am trying to compute the probability (neg_TG>=0).

I’m wondering whether it’s possible to fit a multilevel GPD in STAN to negated data with varying threshold. If yes, why all the chains are not initializing?

**Also, is it possible to introduce the threshold as a parameter in parameter block along with scale and shape parameter instead of as data?**Data has a ragged structure. I am attaching the stan code, r script and data.

data_SG.csv (14.4 KB)

multilevel_model.R (1.9 KB)

stan_code.stan (1.6 KB)

```
functions {
real gpareto_lpdf(vector y, real ymin, real xi, real sigma, real s) {
// generalised Pareto log pdf
real N = s;
real inv_xi = inv(xi);
if (xi<0 && max(y-ymin)/sigma > -inv_xi)
reject("xi<0 and max(y-ymin)/sigma > -1/xi; found xi, sigma =", xi, sigma);
if (sigma<=0)
reject("sigma<=0; found sigma =", sigma);
if (fabs(xi) > 1e-15)
return -(1+inv_xi)*sum(log1p((y-ymin) * (xi/sigma))) -N*log(sigma);
else
return -sum(y-ymin)/sigma -N*log(sigma); // limit xi->0
}
real gpareto_rng(real ymin, real xi, real sigma) {
// generalised Pareto rng
if (sigma<=0)
reject("sigma<=0; found sigma =", sigma);
if (fabs(xi) > 1e-15)
return ymin + (uniform_rng(0,1)^-xi -1) * sigma / xi;
else
return ymin - sigma*log(uniform_rng(0,1)); // limit xi->0
}
}
data {
// number of groups
int<lower=1> K;
//threshold
vector<lower=-10>[K] ymin;
// number of observations
int<lower=1> N;
//observations
vector [N] y;
// group sizes
int<lower=0> s[K];
// ymax
vector<lower=ymin>[K] ymax;
}
parameters {
//sigma mu
real<lower=0>sigma_mu;
//sigma sd
real<lower=0>sigma_sigma;
//y-specific
vector<lower=0>[K] sigma;
vector<lower=0>[K] xi_offset;
}
transformed parameters{
vector[K] xi = -sigma ./ (ymax-ymin) + xi_offset;
}
model {
int pos;
pos = 1;
//priors
sigma_mu ~ normal(0,3);
sigma_sigma ~ normal(0,3);
xi ~ uniform(-1,1);
sigma ~ normal(sigma_mu, sigma_sigma);
for(i in 1:K){
segment(y, pos, s[i]) ~ gpareto(ymin[i], xi[i], sigma[i], s[i]);
pos = pos + s[i];
}
}
****
```