[SciPy-Dev] scipy.stats: algorithm to for ticket 1493

nicky van foreest vanforeest@gmail....
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 [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

In [26]: print  invnorm.cdf(sol, 7.24000019602, scale=2.51913630166)

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.


More information about the SciPy-Dev mailing list