[SciPy-Dev] scipy.stats: algorithm to for ticket 1493
nicky van foreest
vanforeest@gmail....
Sun Apr 22 13:29:08 CDT 2012
Hi,
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 [21]: from scipy.stats import invnorm
In [22]: from scipy import optimize
In [23]: invnorm.xb = 10
In [24]: sol = invnorm.ppf(0.8455, 7.24000019602, scale= 2.51913630166)
In [25]: print sol
25.1878345282
In [26]: print invnorm.cdf(sol, 7.24000019602, scale=2.51913630166)
0.8455
In [27]:
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.
def findppf(q):
# 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: \
invnorm.cdf(x, 7.24000019602,
scale=2.51913630166) - q,\
left, right)
# 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 sol
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.
Nicky
More information about the SciPy-Dev
mailing list