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_poswhich 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;
}
The
student_t_cdf(andstudent_t_lcdfbecause it's just the logged version of the former) is quite poor for large negative and positive values.which gives
and
This function written in the Stan language provides much better precision