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