I am looking for numerically stable implementations of the functions below. Since my application concerns the t distribution, I am using the t distribution here as an example.

Log CDF

# Naive Python implementation of the function I need

import scipy
import numpy as np

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

In SciPy, there is a logcdf, but it is not sufficiently numerically stable. Compare the behavior of the scipy function:

>>> 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])

to the equivalent in R:

# R code
> 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

Tensorflow Probability offers an implementation, but its results resemble the scipy results.

I could not find an implementation in PyTorch (which does not offer the t distribution’s CDF at all).

Quantile function based on log probability:

# Naive Python implementation of the function I need

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))

I could not find quantile functions that operate on log probabilities in SciPy, PyTorch, or Tensorflow Probability.