Implementing Generalized Exteme Value distribution (GEV) constraints

Hi all,

I’m very new to Stan and I’d like your help on a simple model I’ve been setting up. I want to fit some data to the GEV distribution y~GEV(loc, scale, shape). I want to make sure though that the quantity z = 1+shape(y-loc)/scale* remains strictly positive for any combination of loc/scale/shape that are sampled/updated by the algorithm. How can I do that? My model is given below:

functions{

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

}

data {
int<lower=0> N;
vector[N] y;
}

transformed data {
real miny;
real maxy;

miny = min(y);
maxy = max(y);

}

parameters {
real<lower=-1.0,upper=1.0> shape;
real<lower=0> scale;
real<lower=(shape > 0 ? miny : negative_infinity()),
upper=(shape > 0 ? positive_infinity() : maxy) > loc;
}

model {
scale ~ normal(0, 10);
shape ~ normal(0, 10);
loc ~ normal(0, 10);
y ~ gev(loc, scale, shape);
}