如何以数值稳定方式计算t分布的对数CDF和对数概率分位数函数？

我正在寻找数值稳定的方法来计算以下与t分布相关的函数：

对数CDF（Log CDF） #

简易实现 #

import scipy
import numpy as np

def t_log_cdf(x, df):
    p = scipy.stats.t.cdf(x, df=df)
    return np.log(p)

现有实现的数值稳定性问题 #

SciPy提供了t.logcdf方法，但数值稳定性不足。对比测试结果：

SciPy输出：

>>> import scipy
>>> import numpy as np

>>> scipy.stats.t.logcdf(np.linspace(0, 60, 5), df=15)
array([-6.93147181e-01, -9.70638014e-11, -4.21884749e-15,  0.00000000e+00,
        0.00000000e+00])

R语言等效实现输出：

> pt(seq(0, 60, length.out = 5), df = 15, log.p = TRUE)
[1] -6.931472e-01 -9.706380e-11 -4.157404e-15 -1.012752e-17 -1.384523e-19

可见当x较大时，SciPy直接返回0，无法给出精确的极小对数值，而R的实现能保持精度。Tensorflow Probability的实现表现与SciPy类似，PyTorch则未提供t分布的CDF及对数CDF实现。

基于对数概率的分位数函数 #

简易实现 #

import scipy
import numpy as np

def t_log_ppf(log_p, loc, scale):
    p = np.exp(log_p)
    return scipy.stats.t.ppf(p, df=df)

现有库的缺失情况 #

目前SciPy、PyTorch和Tensorflow Probability中，均未找到直接处理对数概率输入的分位数函数，只能先将对数概率转换为普通概率再调用分位数函数，但这种方式在log_p极小时会因数值下溢导致精度损失。

内容的提问来源于stack exchange，提问作者jobrachem