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{
//GammaGPD 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(xu  0, sigma/xi, 1/xi);
}else if(xi==0){
result = ccdfu + exponential_lpdf(xu  1/sigma);
}else{//xi<0
result = ccdfu  xi*beta_lpdf(xi*(xu)/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] = du(sigma/xi)*((cJ/0.1)^xi1);
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 = (1gamma_cdf(u, alpha, beta))/exp(gamma_lpdf(u  alpha, beta));
real xi = algebra_solver(algebra_system, y_guess, theta, x_r, x_i, 1e10, 1e6, 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 GammaGPD Model.
GammaGPD Model is application of Extreme Value Theory, and this is a mixture model of Gamma distribution and Generalized Pareto distribution.
GammaGPD Model is defined by the following equation:
F(x) = \begin{cases} 1(1J(u\;\;\alpha,\beta))(1H(xu\;\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(1J(u))+\log h(xu) & \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: 1e006. 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.

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

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

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

…
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.