Log(0) error, when defining shape parameter with normal distribution containing regression parameters in non-stationary GEV

Modeling
,
1

I have define a hierarchical GEV with a non-stationary shape parameter.

I have defined some samples y with a changing shape parameter shvec[i] and the covariates COV for the shape parameter are simply 1+(3*shvec).

library(rstan)
library(extRemes)

set.seed(10)

shvec  = seq(-0.3,0.2,length.out=20)
y = array(0,dim=c(40,20))

for (i in 1:20) {
  y[,i] = revd(40,loc=20,scale=0.5,shape=shvec[i])
}

# Bayesian hierarchical model
library(rstan)
N   = 40
J   = 20
COV = 1 + (3*shvec)

model = stan_model("GEV.stan")
fit = sampling(model,list(N=N,J=J,y=y,COV=COV),iter=2000,chains=1,seed=10)

The model GEV.stan is below:

functions {
    real gev_lpdf(vector y, real mu, real sigma, real xi) {
        vector[rows(y)] z;
        vector[rows(y)] lp;
        real loglike;
        int N;
         
        N = rows(y);
         
        if (fabs(xi) > 1e-10) {
            for(i in 1:N){
                z[i] = 1 + xi * (y[i] - mu) / sigma;
                lp[i] = (1 + 1 / xi) * log(z[i]) + pow(z[i], -1 / xi);
            }
        loglike = - sum(lp) - N * log(sigma);
        } else {
            for (i in 1:N){
                z[i] = (y[i] - mu) / sigma;
                lp[i] = z[i] + exp(-z[i]);
            }
            loglike = - sum(lp) - N * log(sigma);
        }
        return(loglike);
    }
}


data {
  int<lower=1> N;      
  int<lower=1> J;      
  vector[J] y[N];
  vector[J] COV;
}

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

model {
  for (jj in 1:J) {
    xi[jj] ~ normal(alpha + delta*COV[jj],10);
    y[jj] ~ gev(mu[jj], sigma[jj], xi[jj]);
  }
}

Unfortunaly I get Log probability evaluates to log(0), i.e. negative infinity for this. Anybody know what I’m doing wrong here?

I think it may have something to do with shape parameter xi[jj] ~ normal(alpha + delta*COV[jj],10).

2

Doesn’t it tell you exactly which parameter the error refers to? It can be something that is wrongly (intentionally or unintentionally) set to zero than used to compute the likelihood, or something that is randomly initialized to a very low probability that ends up underflowing to zero. The latter is common when the parameter space is large, but can be easily solved by setting the initial parameter values to values known to have nonzero likelihood.

3

Do not force the sigma parameter to be positive but rather model log(sigma), i.e., replace sigma with exp(sigma) in your likelihood function. See if that fixes the issue…

4

I found my problem.

Obviously, the shape parameter also needs to be constrained:

1 + \xi \left( \frac{y-\mu}{\sigma} \right) > 0

So the code for that:

transformed data {
  vector[J] ymin;
  for (j in 1:J) {
    ymin[j] = min(y[j]);
  }
}

parameters {
  real mu[J];
  real<lower=0> sigma[J];
  real<lower=-sigma[J]/(ymin[J] - mu[J])> xi[J];
  
  real alpha;
  real delta;
}
5

Nice! I wonder why you need to model shape as a function of covariates… Usually, this is not a good idea (see Coles et al. 2001).

1 Like
6

Yes, it’s probably better to keep the covariates in the location or, if needed, the scale parameter.

1 Like