Variable bounds depending on parameter

Modeling
,
1

I would like to fit some rainfall data to a GEV. I am trying different versions of the model, some considering constant parameters, others where parameters present a linear trend and the final idea would be to fit a hierarchichal model where parameters are month specific and include a month-specific trend.

The problem appears when I try to impose the condition of the GEV parameters that defines the support of the observations z:

1 + shape * (z - location) / scale > 0

where location, scale and shape are the parameters of the GEV

For simple models, I do not need to enforce the condition. I only reject the wrong samples, but Stan is able to converge easily. As more parameters are added (time trend, month-specific parameters), Stan cannot even initialize the iterations.

The condition mentioned above implies that:
if shape > 0 --> -Infinity < location < min(z) + scale / shape
if shape < 0 --> max(z) + scale / shape < location < Infinity

Is it possible to reparameterize the problem in order to impose the condition maintaining a one-to-one relation between the unconstrained and the constrained parameters? I have read that it is bad practice to include conditionals in the transformation of the unconstrained parameter into the constrained one, but I am quite blocked trying to find a solution not involving conditionals.

2 Likes
2

You can (and should) do

transformed data {
  real min_z = min(z);
  real max_z = max(z);
}
parameters {
  real shape;
  real<lower=0> scale;
  real<lower = shape > 0 ? negative_infinity() : max_z + scale / shape,
       upper = shape < 0 ? positive_infinity() : min_z + scale / shape> location;
}

Discontinuities are bad for NUTS but in this case, it is continuous at zero. But it is worse to have part of the parameter space have zero density.

1 Like
3

Thank you very much for your response @bgoodri. I think I have not explained myself clearly, so I will try to develop it a little more to see if I am able to transmit my point.

In any case, your answer has clarified me one point that I also had in mind, and generated an additional question. I understand that the way of coding the model that you propose would work for the stationary case if I use an improper prior but, would it be possible to set a proper prior? That is, may I truncate the distribution in the same way:

location ~ normal(0, 1) [ shape > 0 ? negative_infinity() : max_z + scale / shape, 
                          shape < 0 ? positive_infinity() : min_z + scale / shape ]

Going back now to my original question, my problem is that location is not a parameter, but a transformed parameter. That is, location is not sampled directly but generated as a linear combination of other parameters, specifically a constant value location_base and the product of a coefficient alfa and time (linear trend).

data {
    int<lower = 1> N;
    vector[N] time;
    vector[N] z;
}
transformed data {
    real min_z = min(z);
    real max_z = max(z);
}
parameters {
    real location_base;
    real scale_base_raw;
    real shape_base;
    vector[3] alfa;
}
transformed parameters {
    vector<lower = 0>[N] scale;
    vector[N] shape;
    vector<lower = shape > 0 ? negative_infinity() : max_z + scale / shape,
           upper = shape < 0 ? positive_infinity() : min_z + scale / shape>[N] location;
    scale = exp(scale_base_raw + alfa[2] * time);
    shape = shape_base + alfa[3] * time;
    location = location_base + alfa[1] * time;
}
model {
    z ~ gev(location, scale, shape);
}
generated quantities {
    real scale_base;
    scale_base = exp(scale_base_raw);
}

In some cases I am able to set a reparameterization that allows me to sample only valid values for the parameters. For instance, the scale parameter must be positive and thus I generate it as:

scale = exp(scale_base_raw + alfa[2] * time);

No matter the distribution that I use for scale_base_raw and alfa[2] the sampler is able to only generate valid candidates for scale.

However, in the case of the location parameter, because of the discontinuity introduced by the value of shape, I am not being able to formulate a proper reparameterization. I now understand how to set the lower and upper limits, which is great, but I understand that this is not enough to ensure that only valid parameters will be generated. Am I missing some point here? Is reparameterization not the proper strategy to follow here?

In essence, I believe that my question is: how do I enforce, by construction, that a linear combination of parameters (or any other function for that sake) always satisfies a constraint like the one location has in the GEV case.

Sorry for the long message. I hoped to have been able to explain myself in less lines but I may not have acquired the proper lingo yet. Thank you very much for you help and attention.

4

In this case, I think you need to restrict the magnitude of alfa[3] to ensure that the constraint is satisfied. This would be easier if you separated alfa into three different scalars. I would create two functions called lb and ub that return the necessary lower bound and upper bound on alpha_3 respectively. For sufficiently small alpha_3, it should work, or at least it would with informative priors.

functions {
  real lb(real min_z, real max_z, real shape_base, vector time) {

  }
  real ub(real min_z, real max_z, real shape_base, vector time) {

  }
}
parameters {
  real location_base;
  real scale_base_raw;
  real shape_base;
  real alpha_1;
  real alpha_2;
  real<lower=lb(min_z, max_z, shape_base, time),
       upper=ub(min_z, max_z, shape_base, time)> alpha_3;
}
transformed parameters {
    vector[N] location = location_base + alfa_1 * time;
    vector[N] scale = exp(scale_base_raw + alfa_2 * time);
    vector[N] shape = alfa_3 * (shape_base + time);
}
5

Thank you very much @bgoodri, I am going to try the strategy that you propose.

I don’t know why I had understood that setting the bounds was only half of the work and that you needed to enforce that the bound was respected by constructing an adequate prior distribution or reparameterization.

I will let you know if I am able to make this work.

Thanks again.

6

Hi @ManudJ, I am planning to use GEV to study the behaviour of a bunch of meteorological extremes.
Sorry for the silly question, but where did you find the implementation of the log probability of the GEV distribution?
Thanks in advance,
F.

7

A google search for gev log density stan] yielded several results, the first of which has a simple gev density implementation:

https://groups.google.com/forum/#!topic/stan-users/B1x42n0ecQw

1 Like
8

Thank you @Bob_Carpenter . I’ll use the one at the end of the topic you pointed me at.

9

Hi @fabio

I implemented the GEV myself based on different online resources, among which the one linked by @Bob_Carpenter

Hope it works well for you