The gradient of inner_integral is nan 0

Modeling
1

I my attempts to implement custom probability densities I have never encountered this error. I am trying to implement a SICHEL lpmf the reparametrisation here.

I could not find much online on this error

Chain 1: 0.378937 0.0977251 0.355928
-5.8373
-6.85597
-6.92309
-6.44776
-5.93319
-5.69185
-5.18989
-6.28615
-5.78904
-7.88485
-5.74056
-5.36186
-5.49455
-5.13792

Chain 1: Exception: Exception: inner_integral_grad_v: The gradient of inner_integral is nan 0  (in 'model3f1f757baa8e_poisson_GIG' at line 10)
  (in 'model3f1f757baa8e_poisson_GIG' at line 36)

[1] "Error in sampler$call_sampler(args_list[[i]]) : "                                                                                       
[2] "  Exception: Exception: inner_integral_grad_v: The gradient of inner_integral is nan 0  (in 'model3f1f757baa8e_poisson_GIG' at line 10)"
[3] "  (in 'model3f1f757baa8e_poisson_GIG' at line 36)"

I get log probabilities as negative number as expected, and I am not even sampling with the sichel function, just printing

The model:

functions{
  real log_besselk_frac(real v, real z);

  real sichel_lpmf(int y, real gamma, real alpha, real epsilon){
    // https://www.researchgate.net/publication/247109928_Parameter_Estimation_for_the_Sichel_Distribution_and_Its_Multivariate_Extension


    real w = hypot(epsilon, alpha) - epsilon;

    return
      ( (gamma * (log(w) - log(alpha))) - log_besselk_frac(gamma, w) ) +
      (y * (log(epsilon) + log(w) - log(alpha))) - lgamma(y+1) +
      log_besselk_frac(y + gamma, alpha);

  }

}
data{
  int<lower=0> N;
  int y[N];

}
parameters{
  real<lower=-0.5, upper=0.5> gamma;
  real<lower=0, upper=0.5> alpha;
  real<lower=0, upper=0.5> epsilon;
}
model {
gamma ~ normal(0,1);
alpha ~ normal(0,1);
epsilon ~ normal(0,1);

print(gamma, " " , alpha, " " , epsilon);
for(i in 1:N) print(sichel_lpmf(y[i] | gamma, alpha , epsilon));

}
2

I had a quick grep in a fresh version of cmdstan and didn’t find an inner_integral.

Separate out and print the values of every one of the functions you use in sichel_lpmf. That’s the next step. You’ll want to identify exactly what is blowing up.

Even though you are only printing the sichel_lpmf stuff, it’ll still end up on the autodiff stack (if that was part of the question). Everything on the autodiff stack gets its grad function called even if it didn’t directly contribute to the target log density.

3
1 Like
4

Oh, neat, it’s alive. This was for overflow stuff, right? The quadrature tries to evaluate stuff like 1e300 and sometimes things just break?

@stemangiola where’s the integrade1d live? Add some prints around there I s’pose.

5

@martinmodrak now I see, I think it is in the BesselK. Strange in the other SICHEL formulation was not complaining, when I have a minute I check the inputs of besselK

6

Yeah, this is totally my custom code (although it has adapted parts from integral_1d). Will look into this, but it 100% is not a Stan issue.