Skip to content
This repository

BUG: solved issue with ppf #216

Closed
wants to merge 2 commits into from

4 participants

nokfi Denis Laxalde Ralf Gommers Josef Perktold
nokfi
nokfi commented

I issued this first as a branch called ticket1493. This got messed up with the branch repair-typo. Therefore I removed ticket1493, and pushed this new branch. There are also some comments related to the old, but now removed, branch. I include this as one comment here.

nokfi
nokfi commented

I replaced the ppf function. Once this works, I'll remove xa and xb.

josef-pkt commented 8 hours ago
the new ppf function looks good, as discussed on the mailing list, the existing test should cover this case.

possible extra optimization
_ppf_to_solve uses cdf which needs to go through the generic wrapping code each time.
we could use _cdf instead, but that would require additional checking/tests to see whether check_args are properly set (I think they are)

nokfi commented 2 hours ago
As a test I naively replaced self.cdf by self._pdf in _ppf_to_solve, and ran test_continuous_basic.py This fails, so a simple replacement will not suffice.

josef-pkt commented 2 hours ago
"replaced self.cdf by self._pdf" did you actually use _pdf or _cdf ?

my guess was, that, since you limit the range to the interval [.a, .b], _cdf might work, but maybe it would require open intervall (.a, .b).
I tried some cases like this in the past, and sometimes it worked often it broke on a few distributions.

nokfi commented 2 hours ago
On 17 May 2012 21:30, Josef Perktold
reply@reply.github.com
wrote:
"replaced self.cdf by self._pdf" did you actually use _pdf or _cdf ?
To be explicit, I did this:

def _ppf_to_solve(self, x, q,*args):
    return apply(self._cdf, (x, )+args)-q

this self._cdf rather than self.cdf

my guess was, that, since you limit the range to the interval [.a, .b], _cdf might work, but maybe it would require open intervall (.a, .b).
I tried some cases like this in the past, and sometimes it worked often it broke on a few distributions.
I like to include .a and .b in the search interval to cover the cases
q = 0 and q =1. Sure, mathematically speaking the return value of
cdf(x) = 1 can be set to np.inf. However, this is, at least for me,
less informative than self.b. (and likewise for self.a).

So we stick to cdf, rather than _cdf, at least for the moment?
$B!D

josef-pkt commented an hour ago
Yes, stick with .cdf. Worth a try, but if it doesn't work it's unrelated to current pull request.
Your limiting to [a,b] is the right thing to do

(aside: q=0, q=1 is supposed to be handled by generic part, but using left=.a and right=.b if finite, let's brentq search arbitrarily close to the boundary, which is singular in some cases (pdf goes to inf). Although, it might still break at a singularity, or if cdf uses intquad and we want ppf(1e-10))

Denis Laxalde
Collaborator
dlax commented

Does this replace PR #214? If so, you should close the latter @nokfi.

nokfi
nokfi commented
Ralf Gommers
Owner

general_cont_ppf can indeed be removed.

Ralf Gommers
Owner

I didn't see the git question in PR-214 before, but I'm working on disentangling some of these PRs now.

Nicky, I think you are now used to creating different branches for new features, but keep in mind that you should create each branch separately based on master, and not on each other (you did that correctly for this PR now). Otherwise the same commits show up in several PRs.

Ralf Gommers
Owner

The change to test_continuous_basic.py is already in master, that's why this can't be merged automatically. Once this is done, I can rebase and merge it.

nokfi
nokfi commented

Sorry for the mess. Is there something that I should do about this, or is it easy to fix?

Ralf Gommers
Owner

No, it's no problem. This one's not too bad, PR-205 was the one requiring some surgery and that's merged now.

nokfi
nokfi commented
Josef Perktold
Collaborator

docstring line 879 and following still contain explanation for xa xb AFAICS

The only problem left here is the removal of xa and xb from the init. This can possibly break existing code even if it's not necessary anymore.

I'm a bit puzzled that I don't find any xa, xb defined in any of the distributions when the instances are created. A quick look at my own code, also doesn't use it anymore. I was convinced that I changed xa and xb for some distributions.

If we want to be conservative with respect to backwards compatibility, we could add for example a **kwds to the init method, and warn,

if 'xa' in kwds or 'xb' in kwds: 
    import warnings
    <deprecation warning....>

Since I recently advertised xa, xb on stackoverflow, http://stackoverflow.com/questions/10678546/creating-new-distributions-in-scipy it might be better to raise a deprecation warning instead of raising an exception if someone defines xa, xb.

Other than the xa,xb in the init it's a good change. If the test suite passes (including slow tests), then there are no hidden problems with any distribution.

Ralf Gommers
Owner

Better to just not remove xa and xb keywords from the signature in that case, add the deprecation warning for it and open a new ticket to remove them for 0.12. Adding a new **kwds parameter is not so nice.

Ralf Gommers
Owner

I can do that, remove the left over doc and also remove general_cont_ppf and merge this.

Josef Perktold
Collaborator

sounds good,
(I didn't like **kdws much either, but it would have removed xa, xb from the visible signature, although almost no one will look at the init signature)

Josef Perktold
Collaborator

in case it's not obvious: if xa, xb stay as kwd arguments, the defaults should be set to None, so it's easy to check whether someone used a value as argument.

Ralf Gommers
Owner

Opened http://projects.scipy.org/scipy/ticket/1667 to track removal of these params.

Ralf Gommers
Owner

recipinvgauss and foldcauchy were still initialized with an xb parameter by the way.

Josef Perktold
Collaborator

I'm glad you found them, I searched everywhere but only for xa, not for xb. (and I can still trust my memory)

Ralf Gommers
Owner

foldcauchy doesn't seem to have a test for this specific issue though. recipinvgauss apparently wasn't working anyway, according to the FIXME note above it.

Ralf Gommers
Owner

Pushed as ed7f037 and 49ed1d2. Thanks Nicky and Josef.

Ralf Gommers rgommers closed this
jnothman jnothman referenced this pull request from a commit
Commit has since been removed from the repository and is no longer available.
Ralf Gommers rgommers referenced this pull request from a commit in rgommers/scipy
Ralf Gommers DEP: remove deprecated keywords `xa, xb` from continuous distributions.
Closes gh-2192.  See gh-216 for discussion.
388cf0c
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
This page is out of date. Refresh to see the latest.
31  scipy/stats/distributions.py
@@ -434,11 +434,9 @@ def _build_random_array(fun, args, size=None):
434 434
 ##  (needs cdf function) and uses brentq from scipy.optimize
435 435
 ##  to compute ppf from cdf.
436 436
 class general_cont_ppf(object):
437  
-    def __init__(self, dist, xa=-10.0, xb=10.0, xtol=1e-14):
  437
+    def __init__(self, dist, xtol=1e-14):
438 438
         self.dist = dist
439 439
         self.cdf = eval('%scdf'%dist)
440  
-        self.xa = xa
441  
-        self.xb = xb
442 440
         self.xtol = xtol
443 441
         self.vecfunc = sgf(self._single_call,otypes='d')
444 442
     def _tosolve(self, x, q, *args):
@@ -1077,7 +1075,7 @@ def _pdf:
1077 1075
 
1078 1076
     """
1079 1077
 
1080  
-    def __init__(self, momtype=1, a=None, b=None, xa=-10.0, xb=10.0,
  1078
+    def __init__(self, momtype=1, a=None, b=None,
1081 1079
                  xtol=1e-14, badvalue=None, name=None, longname=None,
1082 1080
                  shapes=None, extradoc=None):
1083 1081
 
@@ -1095,8 +1093,6 @@ def __init__(self, momtype=1, a=None, b=None, xa=-10.0, xb=10.0,
1095 1093
             self.a = -inf
1096 1094
         if b is None:
1097 1095
             self.b = inf
1098  
-        self.xa = xa
1099  
-        self.xb = xb
1100 1096
         self.xtol = xtol
1101 1097
         self._size = 1
1102 1098
         self.m = 0.0
@@ -1177,7 +1173,28 @@ def _ppf_to_solve(self, x, q,*args):
1177 1173
         return apply(self.cdf, (x, )+args)-q
1178 1174
 
1179 1175
     def _ppf_single_call(self, q, *args):
1180  
-        return optimize.brentq(self._ppf_to_solve, self.xa, self.xb, args=(q,)+args, xtol=self.xtol)
  1176
+        left = right = None
  1177
+        if self.a > -np.inf:
  1178
+            left = self.a
  1179
+        if self.b < np.inf:
  1180
+            right = self.b
  1181
+
  1182
+        factor = 10.
  1183
+        if  not left: # i.e. self.a = -inf
  1184
+            left = -1.*factor
  1185
+            while self._ppf_to_solve(left, q,*args) > 0.:
  1186
+                right = left
  1187
+                left *= factor
  1188
+            # left is now such that cdf(left) < q
  1189
+        if  not right: # i.e. self.b = inf
  1190
+            right = factor
  1191
+            while self._ppf_to_solve(right, q,*args) < 0.:
  1192
+                left = right
  1193
+                right *= factor
  1194
+            # right is now such that cdf(right) > q
  1195
+
  1196
+        return optimize.brentq(self._ppf_to_solve, \
  1197
+                               left, right, args=(q,)+args, xtol=self.xtol)
1181 1198
 
1182 1199
     # moment from definition
1183 1200
     def _mom_integ0(self, x,m,*args):
5  scipy/stats/tests/test_continuous_basic.py
@@ -305,8 +305,9 @@ def check_sample_meanvar(sm,m,msg):
305 305
 
306 306
 @_silence_fp_errors
307 307
 def check_cdf_ppf(distfn,arg,msg):
308  
-    npt.assert_almost_equal(distfn.cdf(distfn.ppf([0.001,0.5,0.999], *arg), *arg),
309  
-                            [0.001,0.5,0.999], decimal=DECIMAL, err_msg= msg + \
  308
+    values = [0.001,0.5,0.999]
  309
+    npt.assert_almost_equal(distfn.cdf(distfn.ppf(values, *arg), *arg),
  310
+                            values, decimal=DECIMAL, err_msg= msg + \
310 311
                             ' - cdf-ppf roundtrip')
311 312
 
312 313
 @_silence_fp_errors
Commit_comment_tip

Tip: You can add notes to lines in a file. Hover to the left of a line to make a note

Something went wrong with that request. Please try again.