Initial value error (log(0)) for user defined beta prior distribution in hierarchical generalized extreme value (GEV) distribution

I have a user defined Beta distribution.

Beta(\theta+0.5|p,q) = (0.5 + \theta)^{p-1} (0.5 - \theta)^{q-1} / B(p,q) between [-0.5,0.5]

where B(p,q) = \Gamma(p) \Gamma(q) / \Gamma(p+q)

real mbeta_lpdf(real theta, real p, real q) {
      real lpp;
      real b;
      lpp = pow(0.5 + theta,p-1) * pow(0.5 - theta,q-1);
      b = (tgamma(p-1)*tgamma(q-1))/tgamma(p+q-1);
      return log(lpp) - log(b);

This Beta distribution is implemented into a hierarchical generalized extreme value (GEV) distribution for the shape parameter (xi). For the location (mu) and scale (sigma) parameters I use a normal and gamma prior distributions respectively with uniform hyperpriors (alpha = uniform(-infinity,+infinity) and delta = uniform(0,+infinity).

This is the whole model:

functions {
    real gev_lpdf(vector y, real mu, real sigma, real xi) {
        int N = rows(y);
        
        vector[rows(y)] z;
        vector[rows(y)] lp;
        
        
        for(n in 1:N) {
            z[n] = 1 + (y[n] - mu) * xi / sigma;
            lp[n] = (1 + (1 / xi)) * log(z[n]) + pow(z[n], -1/xi);
        }        
        return -sum(lp) -  N*log(sigma);
    }
    
    real mbeta_lpdf(real theta, real p, real q) {
      real lpp;
      real b;
      lpp = pow(0.5 + theta,p-1) * pow(0.5 - theta,q-1);
      b = (tgamma(p-1)*tgamma(q-1))/tgamma(p+q-1);
      return log(lpp) - log(b);
    }
}

data {
  int<lower=0> N;    //observation and duration vector length
  int<lower=1> J;    //number of weather stations
  vector[N] y[J];      //data matrix
}

parameters {
real<lower=-0.5,upper=0.5> xi[J];
real mu[J];
real<lower=0> sigma[J];
real alpha;
real<lower=0> delta;
}

model {
  for(j in 1:J){
    xi[j] ~ mbeta(6.0,9.0);
    mu[j] ~ normal(alpha,25);
    sigma[j] ~ gamma(delta,delta);
    y[j] ~ gev(mu[j],sigma[j],xi[j]);
    }
}

When I run this code unfortunately it throws an error:

...
Chain 1: Rejecting initial value:
Chain 1:   Log probability evaluates to log(0), i.e. negative infinity.
Chain 1:   Stan can't start sampling from this initial value.
Chain 1: 
Chain 1: Initialization between (-2, 2) failed after 100 attempts. 
Chain 1:  Try specifying initial values, reducing ranges of constrained values, or reparameterizing the model.

Where have I gone wrong here?

Are you able to post some example code that runs the model? That would help with debugging. Could you also post the mathematical expression you’re using for the generalized GEV distribution?

At first glance, I’m noticing the 1/xi terms in gev_lpdf, which will blow up when xi is close to zero (it looks like xi can be from -0.5 to 0.5). I’m not sure if that’s the issue, but it should probably be addressed.

Thanks for the reply.

Yes xi should not be equal to zero otherwise the GEV becomes the Gumbel distribution.

Maybe I should constrain xi to be != 0, but where and how to do that I don’t know (I’m fairly new to Stan).

Here is my R code:

library(rstan)

vec =  readRDS("..../file.rds")

col_len = length(vector[,,1])

y = matrix(NA,nrow = 20, ncol = 41)

for (i in 1:20) {
  y[i,] = IDFlist[,1,i]
}

N = 41
J = 20

dir = ".../stackoverflow/"
model = stan_model(paste0(dir,"GEV_hier.stan"))

fit = sampling(model,list(N=N,D=D,J=J,y=y,d=d),iter=500,chains=1,seed=10)

Thanks! I can’t run that code though because it depends on some data files that I don’t have. Maybe you could simulate some sample data as an example? (Fitting to simulated data is a good thing to do anyway, as you can see if the estimated parameters match the simulation parameters.)

As good general debugging strategy when you have an error message but you’re not sure what exactly is causing it is to start removing model components until you get to a simpler model that works. Then, build slowly back up to the original model. For example, you could remove the mbeta prior on xi. That is, comment out

xi[j] ~ mbeta(6.0,9.0);

Then run the code and see if you still have an error. If you do, keep removing model components – maybe assign a simpler response distribution to y instead of gev, for example. If we can localize the problem it will be much easier to troubleshoot.

Thanks for the advice.

I’ve rewritten the code and added simulated GEV samples.

Please see below:

library(rstan)

y = matrix(NA,nrow = 20, ncol = 41)

mu = runif(20,15,20)
sigma = runif(20,3,5)
xi = runif(20,-0.5,0.5)

library(extRemes)

for (i in 1:20) {
  y[i,] = revd(41,loc=mu[i],scale=sigma[i],shape=xi[i],type="GEV")
}

N = 41
J = 20


dir = ".../stackoverflow/"
model = stan_model(paste0(dir,"GEV_hier.stan"))

fit = sampling(model,list(N=N,D=D,J=J,y=y),iter=500,chains=1,seed=10)

I was able to reproduce the error with this code.

Thanks for the example code.

It looks like the issue is that z[n] can be negative depending on the value of mu, sigma, and xi, which then causes log(z[n]) to be undefined in the gev_lpdf calculation.

I think this comes down to how the support of the GEV distribution depends on the values of mu, sigma, and xi. (See “Support” in the sidebar here: Generalized extreme value distribution - Wikipedia).

One way this could be fixed by constraining the values of mu, sigma, and xi so that the observations are always within the support of the GEV. You might look to this post: Bayesian data analysis of extreme values with Stan and rstan. – Georgios Boumis, which appears to take the approach of constraining xi.

1 Like

Excellent! That solved the problem.