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