# Defining a hierarchical GEV with a modified Beta prior

I’ve defined a hierarchical generalized extreme value distribution and I’m trying to implement a beta distribution with a slight modification as a prior for this.

The beta distribution is:

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

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

I have the following Stan code for this:

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 betamod_lpdf(vector y, real p, real q) {
int N = rows(y);
vector[rows(y)] lp;
for(n in 1:N) {
lp[n] = pow(0.5+y[n],p-1) * pow(0.5-y[n],q-1) / lbeta(p,q);
}
return -sum( lp );
}
}

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

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

real gamma_1;
real<lower=0> lambda_1;
real gamma_2;
real<lower=0> lambda_2;
}

model {
//hierarchical hyperprior
gamma_1 ~ normal(0,25);
lambda_1 ~ exponential(1);
gamma_2 ~ normal(0,25);
lambda_2 ~ exponential(1);

for(j in 1:J){
xi[j] ~ betamod(6.0,9.0);
mu[j] ~ normal(gamma_1,lambda_1);
sigma[j] ~ lognormal(gamma_2,lambda_2);
y[j] ~ gev(mu[j],sigma[j],xi[j]);
}
}


So, first of all, when I check this it throws an error saying:

SYNTAX ERROR, MESSAGE(S) FROM PARSER:
No matches for:

real ~ betamod(real, real)

Available argument signatures for betamod:

vector ~ betamod(real, real)

Real return type required for probability function.
error in 'model23b411265b72_hiergev' at line 51, column 29
-------------------------------------------------
49:
50:   for(j in 1:J){
51:     xi[j] ~ betamod(6.0,9.0);
^
52:     mu[j] ~ normal(gamma_1,lambda_1);
-------------------------------------------------

Error in stanc(filename, allow_undefined = TRUE) :
failed to parse Stan model 'hiergev' due to the above error.


And also I wanted to check if the defined function for… Beta(\theta+0.5|p,q) = (0.5 + \theta ) ^ {p-1} (0.5 - \theta)^{q-1} / B(p,q) between [-0.5,0.5]
…is also correct (see below)

    real betamod_lpdf(vector y, real p, real q) {
int N = rows(y);
vector[rows(y)] lp;
for(n in 1:N) {
lp[n] = pow(0.5+y[n],p-1) * pow(0.5-y[n],q-1) / lbeta(p,q);
}
return -sum( lp );
}


Hi :)

In your function definition, you defined y to be a vector, that is why you are getting the error message I think.

You could either change betamod_lpdf(vector y for betamod_lpdf(real y, or avoid the loop by defining xi ~ betamod(6.0, 9.0).

Hope that helps!
Lucas

Chain 1: Rejecting initial value: