Sampling statement does not work at all!

General
1

Hello Stan professors.

I would like you to help me.
A strange error occurred, and I’m in trouble with it.
I wrote the following Stan program:

functions{
  //Gamma-GPD Model from Cibele(2004)
  real ggpd_lpdf(real x, real alpha, real beta, real u, real sigma, real xi){
    real result;
    real ccdfu = gamma_lccdf(u | alpha, beta);
    if(x <= u){
      result = gamma_lpdf(x | alpha, beta);
    }else{//x>u
      if(xi>0){
        result = ccdfu + pareto_type_2_lpdf(x-u | 0, sigma/xi, 1/xi);
      }else if(xi==0){
        result = ccdfu + exponential_lpdf(x-u | 1/sigma);
      }else{//xi<0
        result = ccdfu - xi*beta_lpdf(-xi*(x-u)/sigma | 1, -1/xi)/sigma;
      }
    }
    return result;
  }

  
  //solve function
  vector algebra_system(vector y,//xi
                      vector theta,//alpha, beta, u, d
                      real[] x_r,
                      int[] x_i){
    real xi = y[1];
    
    real alpha = theta[1];
    real beta = theta[2];
    real u = theta[3];
    real d = theta[4];
    
    real cJ = 1 - gamma_cdf(u, alpha, beta);
    real sigma = cJ/exp(gamma_lpdf(u | alpha, beta));

    vector[1] result;//xi
    
    result[1] = d-u-(sigma/xi)*((cJ/0.1)^xi-1);
    return result;
  }
}
data{
  int n;
  real Data[n];
  real u_guess;
}
transformed data{
  vector[1] y_guess = [0.5]';//guess value of xi
  real x_r[0];
  int x_i[0];
}
parameters{
  real<lower = 0> alpha;
  real<lower = 0> beta;
  real<lower = 0, upper = max(Data)> u;
  real<lower = u> d;
}
transformed parameters{
  vector[4] theta = [alpha, beta, u, d]';
  real sigma = (1-gamma_cdf(u, alpha, beta))/exp(gamma_lpdf(u | alpha, beta));
  real xi = algebra_solver(algebra_system, y_guess, theta, x_r, x_i, 1e-10, 1e-6, 1e5)[1];
}
model{
  print("u = ", u);
  for(i in 1:n){
    Data[i] ~ ggpd(alpha, beta, u, sigma, xi);
  }
  
  alpha ~ gamma(1,1);
  beta ~ gamma(1,1);
  u ~ uniform(19.9, 20.1);//Here is a cause of trouble.
  d ~ gamma(12, 0.6);
}
generated quantities{
  
}

What I try to do is Implementation of Gamma-GPD Model.
Gamma-GPD Model is application of Extreme Value Theory, and this is a mixture model of Gamma distribution and Generalized Pareto distribution.
Gamma-GPD Model is defined by the following equation:

F(x) = \begin{cases} 1-(1-J(u\;|\;\alpha,\beta))(1-H(x-u\;|\sigma,\xi)) & \text{where}\;x>u \\ J(x\;|\;\alpha,\beta) & \text{otherwise} \end{cases}

In here, J and H denote CDF of Gamma distribution and Generalized Pareto distribution respectively.
And clearly the log likelihood of the distribution is:

\log f(x) = \begin{cases} \log(1-J(u))+\log h(x-u) & \text{where}\;x>u \\ \log j(x) & \text{otherwise} \end{cases}

In here, j and h denote PDF of …
in addition, h changes shape depending on sign of \xi.

h(x) = \begin{cases} f_{\text{ParetoType2}(0,\sigma/\xi,1/\xi)}(x) & \text{where}\;\xi>0\\ f_{\text{Exponential}(1/\sigma)}(x) & \text{where}\;\xi=0\\ (-\xi/\sigma)f_{\text{Beta}(1,-1/\xi)}(-\xi x/\sigma) & \text{where}\;\xi<0 \end{cases}

I expect to have written correct definition in functions block.

And now, probably I write correct code in all blocks. (I have an anxiety, though.)
Moreover, I succeeded compile this code.

But, I always get following error message:

Unrecoverable error evaluating the log probability at the initial value.
Exception: algebra_solver: the norm of the algebraic function is: 0.31713 but should be lower than the function tolerance: 1e-006. Consider increasing the relative tolerance and the max_num_steps. (in ‘model1ac2bfc7fa3_if’ at line 64)

When I put a print statement in algebra_system to find a bug, I found strange behavior.

  1. theta=[0.25219,0.980642,245.788,246.106], result=[0.31713]

  2. theta=[3.5211,0.141341,127.991,128.164], result=[0.1727]

  3. theta=[2.33726,0.169544,104.544,106.857], result=[2.31328]

  4. …

Although theta = (alpha, beta, u, xi), u is obviously large, because prior distribution of u is defined as uniform(19.9, 20.1) (In actual modeling, I’m planning to use appropriate prior distribution for example Normal distribution having a large variance)

Please tell me the solution to this (or any other) problem.
I will show you the version of R, rstan, etc when I was requested.

2

Hi,
Stan expects that all parameters have support (non-zero posterior density) everywhere within their constraints. However, u has zero density outside [19.9,20.1]. Uniform prior should be expressed with constraints, e.g.:

real<lower = 19.9, upper = 20.1> u;

while the uniform sampling statement is used basically as a comment. In general, uniform priors don’t work well with Stan, you might be better of with something like u ~ normal(20,0.1).

3

Thank you your reply, @martinmodrak!!

But No, I have already tried u ~ normal(20,0.1).
However print statement outputs strange large values such about u = 200, 100 and, same error message occurred.
Frankly speaking, sampling statement of u is not work.
In addition, from the model assumption, (as u means threshold of Data, )

real<lower = 0, upper = max(Data)> u;

should be correct.

The most strange thing is that sampling statement of parameter other than u works correctly.
I cannot find cause of error…

4

I think you are misunderstanding some fundamentals of the way Stan works. This is OK - Stan is not exactly intuitive. However, sampling statements do not “work” in isolation, they simply add a value to the log posterior density. It is this value and its derivatives (and hence combination of all sampling statements together) that drive sampling in Stan. What most likely happens is that either

  1. The algorithm initializes u far from the plausible region around 20 (by default, inits are uniform in [-2,2] on the unconstrained scale - see the manual for details). Since the plausible region is narrow and the posterior geometry difficult, the sampler is unable to move towards this region; or
  2. Other parts of the model force u to grow beyond expected limits.

If 1) is the case you can change the default init (see manual for RStan or other interface you are using how to specify inits). Or - provided you use normal prior for u, you could reparametrize with something like:

parameters {
  real u_raw; 
}

transformed parameters {
  real u = u_raw * 0.1 + 20;
}

model {
  u ~ normal(0,1);
}

The reparemtrization could also help if 2) is the case. To check for 2) you should try recovering parameters from simulated data.

In general, you should start with a very simple model, check that it recovers parameters from simulated data and then add complexity step by step, always checking with simulated data. This is the only way to make sure you didn’t make any mistakes in the coding. Starting with such a complex model straight away is a recipe for problems. The model may actually be difficult for Stan thanks to the discrete switches in \xi (unless the posterior is smooth on these boundaries).

I also have to insist, that whenever you use uniform(a,b) sampling statement, you need to add constraints <lower=a,upper=b on the respective variable, otherwise sampling may not work. Unfortunately, this is not explicitly stated in the Stan manual, but it is true. In particular, constraints are taken into account when initializing the values of parameters, so u is unlikely to be initialized to a value between a and b unless the constraint is given.

5

I see…
I surely don’t understand Stan works well.
But as you say, today, I replaced uniform prior with normal prior, and in order to simplify the model, I daringly removed algebra_solver from my program code.

Then fortunately (or unfortunately?), Stan worked and I got posterior distribution!!!
For simulated data, the estimation result was good.

Although I don’t understand the cause of this issue perfectly. I accept this nice result.
After I try various in this model, I retry to use algebra_solver again.

Thank you for your suggestion.

6

I would also suggest trying some small model that contains only the algebra_solver before adding the solver to the large model.
I am glad to have helped you and I wish you best of luck with this model!

2 Likes