# Student_t quantile function

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

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

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