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

[nicky van foreest]

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