import numpy as np from scipy import stats, optimize def findppf(dist, q, *args, **kwds): if q <= np.finfo(float).eps: return -np.inf # or should we return xa? if q >= 1. - np.finfo(float).eps: return np.inf # or should we return xb? # in case the user has been fiddling with xa and xb, what to choose? # choice 1, at the expense of being broken for some distributions? left, right = -10., 10. # choice 2, at the risk that the user has been fiddling with xa and # xb. if xb < 0 or xa > 0 the algo below will not work, unless we # include an if test on the sign of xa and xb. left, right = dist.xa, dist.xb # search until cdf(left) < q counter = 0 while dist.cdf(left, *args, **kwds) > q: left *= 2 counter += 1 print counter, left # search until cdf(right) > q while dist.cdf(right, *args, **kwds) < q: right *= 2 counter += 1 print counter, right return optimize.brentq(lambda x: dist.cdf(x, *args, **kwds) - q, left, right) print 'exponential' mu = 3. sol = findppf(stats.expon, 0.8455, scale = mu) print sol sol = findppf(stats.expon, 1., scale = mu) print sol sol = findppf(stats.expon, 1. - 1e-15, scale = mu) print sol print print 'invgauss' s = 7.24000019602 sol = findppf(stats.invgauss, 0.8455, s) print sol sol = findppf(stats.invgauss, 1-1e-8, s) print 'roundtrip', 1-1e-8, sol, stats.invgauss.cdf(sol, s) print 1e-30, stats.invgauss.cdf(findppf(stats.invgauss, 1e-30, s), s) print '\nt' print findppf(stats.t, 1-1e-8, s), stats.t.ppf(1-1e-8, s) print findppf(stats.t, 1e-8, s), stats.t.ppf(1e-8, s) print '\ncauchy' print findppf(stats.cauchy, 1e-8), stats.cauchy.ppf(1e-8) print '\nf' print findppf(stats.f, 1-1e-8, 2, 10), stats.f.ppf(1-1e-8, 2, 10)