Student T CDF precision for large values

[spinkney]

The student_t_cdf (and student_t_lcdf because it's just the logged version of the former) is quite poor for large negative and positive values.

out <- matrix(nr = 30, nc = 3)
a <- -1
for (i in 1:30) {
  out[i, 1] <- a
  out[i, 2] <- stan_student_t_cdf(a, 0.1)
  out[i, 3] <- my_pt_stan(a, 0.1)
  a <- a * 10
}

out

out_pos <- matrix(nr = 30, nc = 3)
a <- 1
for (i in 1:30) {
  out_pos[i, 1] <- a
  out_pos[i, 2] <- stan_student_t_cdf(a, 0.1)
  out_pos[i, 3] <- my_pt_stan(a, 0.1)
  a <- a * 10
}

out_pos

which gives

      [,1]       [,2]         [,3]
 [1,] -1e+00 0.41632824 0.4163282425
 [2,] -1e+01 0.33152829 0.3315282881
 [3,] -1e+02 0.26334911 0.2633491054
 [4,] -1e+03 0.20918568 0.2091856843
 [5,] -1e+04 0.16616210 0.1661620957
 [6,] -1e+05 0.13198724 0.1319872442
 [7,] -1e+06 0.10483700 0.1048411947
 [8,] -1e+07 0.08371484 0.0832783211
 [9,] -1e+08 0.00000000 0.0661503218
[10,] -1e+09 0.00000000 0.0525450683
[11,] -1e+10 0.00000000 0.0417380314
[12,] -1e+11 0.00000000 0.0331536968
[13,] -1e+12 0.00000000 0.0263349174
[14,] -1e+13 0.00000000 0.0209185685
[15,] -1e+14 0.00000000 0.0166162096
[16,] -1e+15 0.00000000 0.0131987244
[17,] -1e+16 0.00000000 0.0104841195
[18,] -1e+17 0.00000000 0.0083278321
[19,] -1e+18 0.00000000 0.0066150322
[20,] -1e+19 0.00000000 0.0052545068
[21,] -1e+20 0.00000000 0.0041738031
[22,] -1e+21 0.00000000 0.0033153697
[23,] -1e+22 0.00000000 0.0026334917
[24,] -1e+23 0.00000000 0.0020918568
[25,] -1e+24 0.00000000 0.0016616210
[26,] -1e+25 0.00000000 0.0013198724
[27,] -1e+26 0.00000000 0.0010484119
[28,] -1e+27 0.00000000 0.0008327832
[29,] -1e+28 0.00000000 0.0006615032
[30,] -1e+29 0.00000000 0.0005254507

and

  [,1]      [,2]      [,3]
 [1,] 1e+00 0.5836718 0.5836718
 [2,] 1e+01 0.6684717 0.6684717
 [3,] 1e+02 0.7366509 0.7366509
 [4,] 1e+03 0.7908143 0.7908143
 [5,] 1e+04 0.8338379 0.8338379
 [6,] 1e+05 0.8680128 0.8680128
 [7,] 1e+06 0.8951630 0.8951588
 [8,] 1e+07 0.9162852 0.9167217
 [9,] 1e+08 1.0000000 0.9338497
[10,] 1e+09 1.0000000 0.9474549
[11,] 1e+10 1.0000000 0.9582620
[12,] 1e+11 1.0000000 0.9668463
[13,] 1e+12 1.0000000 0.9736651
[14,] 1e+13 1.0000000 0.9790814
[15,] 1e+14 1.0000000 0.9833838
[16,] 1e+15 1.0000000 0.9868013
[17,] 1e+16 1.0000000 0.9895159
[18,] 1e+17 1.0000000 0.9916722
[19,] 1e+18 1.0000000 0.9933850
[20,] 1e+19 1.0000000 0.9947455
[21,] 1e+20 1.0000000 0.9958262
[22,] 1e+21 1.0000000 0.9966846
[23,] 1e+22 1.0000000 0.9973665
[24,] 1e+23 1.0000000 0.9979081
[25,] 1e+24 1.0000000 0.9983384
[26,] 1e+25 1.0000000 0.9986801
[27,] 1e+26 1.0000000 0.9989516
[28,] 1e+27 1.0000000 0.9991672
[29,] 1e+28 1.0000000 0.9993385
[30,] 1e+29 1.0000000 0.9994745

This function written in the Stan language provides much better precision

 real my_pt_stan(real x, real n) {
  int lower_tail = 1;
  real val;
  
  if (n <= 0) {
    return not_a_number();
  }
  if (is_inf(x)) {
    return not_a_number();
  }
  if (is_inf(n)) {
    if(lower_tail == 1)
      return Phi(x);
    else
      return exp(normal_lccdf(x | 0, 1));
  }

  real nx = 1 + (x / n) * x;
  if (nx > 1e100) {
    real lval = -0.5 * n * (2 * log(abs(x)) - log(n)) - lbeta(0.5 * n, 0.5) - log(0.5 * n);
    val = exp(lval);
  } else {
    val = n > x * x ? exp(beta_lccdf(x * x / (n + x * x) | 0.5, n / 2)) : beta_cdf(1 / nx | n / 2, 0.5);
  }

  if (x <= 0) {
    lower_tail = 0;
  }

    val /= 2;
    if (lower_tail == 1) {
      return 0.5 - val + 0.5;
    } else {
      return val;
    }
    
  return val;
}