Defining a hierarchical GEV with a modified Beta prior

Modeling
,
1

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

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

3

Thanks for the reply.

Now the problem is that I get:

Chain 1: Rejecting initial value:
Chain 1:   Log probability evaluates to log(0), i.e. negative infinity.
Chain 1:   Stan can't start sampling from this initial value.