Variable precision causing bug?

General
1

I am converting some SICHEL formulas to stan.

A formula is particularly sensible to variable precision, when R calculates the input values I get good results, when stan calculates them I get error.

To replicate, no need data

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

  real[] tofySICHEL2(int[] y, real mu, real sigma, real nu, real lbes, real cvec){

  int maxyp1 = max(y) + 1;
  vector[maxyp1] tofY;
  int iy;
  real alpha;
  real sumT;
  real ans[size(y)];

  for (i in 1:size(y))  {
    iy = y[i]+1;
    tofY[1] = (mu/cvec)*pow((1.0+2.0*sigma*mu/cvec),(-0.5))*exp(lbes);
    alpha = sqrt(1.0 + 2.0*sigma*mu/cvec)/sigma;
    sumT = 0.0;

    for (j in 2:iy)  {
      tofY[j] = ( cvec*sigma*(2.0*((j-1.0)+nu)/mu) + (1.0/tofY[j-1])) * pow(mu/(sigma*alpha*cvec),2);
      sumT += log(tofY[j-1]);

      print(cvec, " ",sigma, " ",2.0, " ",j-1.0, " ",nu, " ",mu, " ",tofY[j-1], " ",mu, " ",sigma, " ",alpha, " ",cvec, " | ", tofY[j-1]);

    }
    ans[i] = sumT;
  }
  return ans;
}

  real sichel(int[] y, real nu, real mu, real sigma){


      if (sigma > 10000 && nu > 0) return neg_binomial_2_log_lpmf(y | mu, 1.0/nu);
      else{

        real nu_1 = nu + 1.0;
        real one_over_sigma = 1.0/sigma;
        real log_bes_nu_one_over_sigma = log_besselk_frac(nu, one_over_sigma);
        real cvec = exp(log_besselk_frac(nu_1, one_over_sigma) - log_bes_nu_one_over_sigma);
        real alpha = sqrt(1.0 + 2.0 * sigma * (mu/cvec))/sigma;
        real log_bes_nu_alpha = log_besselk_frac(nu, alpha);
        real lbes = log_besselk_frac(nu_1, alpha) - log_bes_nu_alpha;
        real sumlty[size(y)] = tofySICHEL2(y, mu,sigma,nu,lbes,cvec);

        real lp[size(y)];

      print(sumlty);

       for(i in 1:size(y))
          lp[i] =  -lgamma(y[i] + 1) - nu * log(sigma * alpha) + sumlty[i] + log_bes_nu_alpha - log_bes_nu_one_over_sigma;

        return sum(lp);

      }
  }

}

parameters{
  real<upper=50> nu;
  real<lower=0> mu;
  real<lower=0> sigma;
}
model {
int x[1];
real rr;
x[1] = 100;

   print(tofySICHEL2(x, 10, 0.1, -30, -0.741379920307050355177125311456620693206787109375, 0.1674189575753287917425637942869798280298709869384765625)); // WORKS
   rr = sichel(x, -30, 10 , 0.1); // CALLED FROM INSIE DOES NT

}

Result:

direct call with R precision

print(tofySICHEL2(x, 10, 0.1, -30, -0.741379920307050355177125311456620693206787109375, 0.1674189575753287917425637942869798280298709869384765625));
[325.391]

Call of the function withing another stan function

real sumlty[size(y)] = tofySICHEL2(y, mu,sigma,nu,lbes,cvec);
[-nan]

P.S. I should add that I think is a precision problem because I can replicate the bug using the proper c++ function called from R

  .C(
        "tofySICHEL2",
        as.double(100),
        as.double(10),
    as.double(0.1),
        as.double(-30),
        as.double(-0.741379920307050355177125311456620693206787109375),
        as.double( 0.1674189575753287917425637942869798280298709869384765625),
        ans = double(length(100)),
       as.integer(length(100)),
        as.integer(max(100) + 1),
     PACKAGE = "gamlss.dist"
      )$ans

[1] 325.3915

  .C(
    "tofySICHEL2",
    as.double(100),
    as.double(10),
    as.double(0.1),
    as.double(-30),
    as.double(-0.7413799),
    as.double( 0.167419),
    ans = double(length(100)),
    as.integer(length(100)),
    as.integer(max(100) + 1),
    PACKAGE = "gamlss.dist"
  )$ans

-nan