Fitting a multilevel GPD to negated data

Modeling
1

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

2023-07-13 23_17_02-04-Paper_3_BGPD - RStudio
2023-07-13 23_17_02-04-Paper_3_BGPD - RStudio763×872 57.5 KB

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];
  }
}

****
2

You should be getting more debug information out of this, like the line of the failure. I’m not sure why you’re not seeing that.

You need to initialize chains at a value of the parameters where the density is positive (i.e., log density is finite). Your model’s too complicated for me to understand easily, but you’re presumably rejecting because of illegal values.

One thing you can do if you can’t get more debug info out of the command (I don’t know how to do that), then you can write print() statements into your code where you can do print("place marker", target()) and see when target() becomes -infinity to try to diagnose where things are going wrong.

1 Like
3

Hi @Bob_Carpenter

Sir, thank you very much for taking the time to read my question.

I am a beginner in stan and do not have much coding knowledge. In my data set, 3 subgroups have fewer observations (less than 10). After removing those, the problem got fixed, and all the chains are running and mixing well. Is that might have been the issue?

Further, I would like your comment on my second question “Is it possible to introduce the threshold as a parameter in the parameter block along with scale and shape parameters instead of as data?”