[SciPy-Dev] scipy.stats: algorithm to for ticket 1493
nicky van foreest
Sun Apr 22 13:29:08 CDT 2012
I saw Ralf's suggestions on issues of scipy.stats to work on. In
another project I had to solve a similar problem as the one mentioned
in http://projects.scipy.org/scipy/ticket/1493. I like to propose the
algorithm below to tackle this. I am looking forward to your feedback.
Once the algo seems to work impeccably, I'll try to convert it such
that it can be incorporated in scipy.stats and add test cases.
Lets first test the example of the ticket:
In : from scipy.stats import invnorm
In : from scipy import optimize
In : invnorm.xb = 10
In : sol = invnorm.ppf(0.8455, 7.24000019602, scale= 2.51913630166)
In : print sol
In : print invnorm.cdf(sol, 7.24000019602, scale=2.51913630166)
Weird that this leads to a solution, as the solution is not in the
search interval. Trying xb = 5 leads to the reported error.
Now a solution:
# Proposed solution.
# A quintessential property in the algorithm is that the cdf is an
# non-decreasing function.
# search until xb is large enough
left, right = invnorm.xa, invnorm.xb
while invnorm.cdf(right, 7.24000019602, scale=2.51913630166) < q:
left = right
right *= 2
return optimize.brentq(lambda x: \
scale=2.51913630166) - q,\
# Perhaps increasing with a larger factor than 2 is faster, as brentq
# convergest very fast. Taking `right` too large is problably not a real
# problem. On the other hand, when `right` increases too fast, cdf(right) may
# become numerically equal to 1. So, what would be a good factor?
sol = findppf(0.8455)
print invnorm.cdf(sol, 7.24000019602, scale=2.51913630166)
This code (I ran it in a file, not in ipython) gives the right
results, at least for these values.
More information about the SciPy-Dev