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.