Student_t quantile function

Modeling
1

Here’s the student_t_qf that handles dof for all cases. Inspired by

functions {
 real student_t_qf(real p, real ndf) {
    real eps = 1e-12;
    real d_epsilon = 2.220446e-16;
    real d_max = 1.797693e+308;
    real d_min = 2.225074e-308;
    int d_mant_dig = 53;
    real q;
    
    if (is_nan(p) || is_nan(ndf)) {
      return p + ndf;
    }
    
    if (ndf <= 0) {
      reject("Invalid value for ndf");
    }
    
    if (ndf < 1) {
      // Find the upper and lower bounds
      real accu = 1e-13;
      real Eps = 1e-11;
      if (p > 1 - d_epsilon) {
        return positive_infinity();
      }
      real pp = min({1 - d_epsilon, p * (1 + Eps)});
      real ux = 1;
      while (student_t_cdf(ux | ndf, 0., 1.) < pp) {
        ux = ux * 2;
      }
      pp = p * (1 - Eps);
      real lx = -1;
      while (student_t_cdf(lx | ndf, 0., 1.) > pp) {
        lx = lx * 2;
      }
      
      // Find the quantile using interval halving
      real nx = 0.5 * (lx + ux);
      int iter = 0;
      while ((ux - lx) / abs(nx) > accu && iter < 1000) {
        iter += 1;
        if (student_t_cdf(nx | ndf, 0., 1.) > p) {
          ux = nx;
        } else {
          lx = nx;
        }
        nx = 0.5 * (lx + ux);
      }
      return 0.5 * (lx + ux);
    }
    
    if (ndf > 1e20) {
      return inv_Phi(p);
    }
    
    int neg = p < 0.5 ? 1 : 0;
    int is_neg_lower = neg;
    real P = neg == 1 ? 2 * p : 2 * (0.5 - p + 0.5);
    
    P = min({max({P, 0}), 1});
    
    if (abs(ndf - 2) < eps) {
      if (P > d_min) {
        if (3 * P < d_epsilon) {
          q = 1 / sqrt(P);
        } else if (P > 0.9) {
          q = (1 - P) * sqrt(2 / (P * (2 - P)));
        } else {
          q = sqrt(2 / (P * (2 - P)) - 2);
        }
      } else {
        q = positive_infinity();
      }
    } else if (ndf < 1. + eps) {
      if (P == 1.) {
        q = 0;
      } else if (P > 0) {
        q = 1 / tan(pi() * p / 2);
      } else {
        q = negative_infinity();
      }
    } else {
      real x = 0;
      real y;
      real log_P2 = 0;
      real a = 1 / (ndf - 0.5);
      real b = 48 / (a * a);
      real c = ((20700 * a / b - 98) * a - 16) * a + 96.36;
      real d = ((94.5 / (b + c) - 3) / b + 1) * sqrt(a * pi() / 2) * ndf;
      
      y = pow((d * P), (2.0 / ndf));
      int P_ok = y >= d_epsilon ? 1 : 0;
      
      if (P_ok != 1) {
        log_P2 = is_neg_lower == 1 ? log(p) : log1m_exp(p);
        x = (log(d) + log2() + log_P2) / ndf;
        y = exp(2 * x);
      }
      
      if ((ndf < 2.1 && P > 0.5) || y > 0.05 + a) {
        if (P_ok == 1) {
          x = inv_Phi(0.5 * P);
        } else {
          x = inv_Phi(log_P2);
        }
        
        y = square(x);
        
        if (ndf < 5) {
          c += 0.3 * (ndf - 4.5) * (x + 0.6);
        }
        
        c = (((0.05 * d * x - 5) * x - 7) * x - 2) * x + b + c;
        y = (((((0.4 * y + 6.3) * y + 36) * y + 94.5) / c - y - 3) / b + 1)
            * x;
        y = expm1(a * square(y));
        q = sqrt(ndf * y);
      } else if (P_ok != 1 && x < -log2() * d_mant_dig) {
        q = sqrt(ndf) * exp(-x);
      } else {
        y = ((1
              / (((ndf + 6) / (ndf * y) - 0.089 * d - 0.822) * (ndf + 2) * 3)
              + 0.5 / (ndf + 4))
             * y - 1)
            * (ndf + 1) / (ndf + 2) + 1 / y;
        
        q = sqrt(ndf * y);
      }
      
      if (P_ok == 1) {
      int it = 0;
      while (it < 10) {
        y = exp(student_t_lpdf(q | ndf, 0., 1.));
        if (y <= 0 || is_inf(y)) {
          break;
        }
        
        real t = (exp(student_t_lccdf(q | ndf, 0., 1.)) - P / 2) / y;
        if (abs(t) <= 1e-14 * abs(q)) {
          break;
        }
        
        q = q + t * (1. + t * q * (ndf + 1) / (2 * (q * q + ndf)));
        
        it += 1;
      }
      }
    }
    
    if (neg) {
      q = -q;
    }
    
    return q;
  }
  
  vector student_t_qf(vector u, real nu) {
    int N = num_elements(u);
    vector[N] x;
    
    for (i in 1 : N) {
      x[i] = student_t_qf(u[i], nu);
    }
    
    return x;
  }

}
3 Likes
2

Is there any chance of code like this getting into stan proper?

3

The problem with the above is that it follows R’s implementation and we can’t add that into Stan because of licensing restrictions. That’s not a death knell though because someone can just rewrite the below parts with their own code:

  1. for ndf < 1 write your own halving algorithm
  2. Hill’s algorithm can be used verbatim and starts at
} else { 
 real x = 0;
 real y = 0; ...
  1. the while loop at the end needs to be re-written that starts at
int it = 0;
while (it < 10) { ...

The other issue for any devs is to just let autodiff differentiate through this and don’t write your own because the derivative wrt to nu is nasty and slow.

2 Likes