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)`

.