Skip to content

Commit ed40f3c

Browse files
committed
MAINT: stats.dpareto_lognorm._cdf: treat special case at x=0
1 parent d821271 commit ed40f3c

File tree

1 file changed

+18
-14
lines changed

1 file changed

+18
-14
lines changed

scipy/stats/_continuous_distns.py

Lines changed: 18 additions & 14 deletions
Original file line numberDiff line numberDiff line change
@@ -1923,20 +1923,24 @@ def _logpdf(self, x, u, s, a, b):
19231923
return out[()]
19241924

19251925
def _logcdf(self, x, u, s, a, b):
1926-
log_y, m = np.log(x), u # compare against [1] Eq. 2
1927-
z = (log_y - m) / s
1928-
x1 = a * s - z
1929-
x2 = b * s + z
1930-
t1 = self._logPhi(z)
1931-
t2 = self._logphi(z)
1932-
t3 = (np.log(b) + self._logR(x1))
1933-
t4 = (np.log(a) + self._logR(x2))
1934-
t1, t2, t3, t4, one = np.broadcast_arrays(t1, t2, t3, t4, 1)
1935-
# t3 can be smaller than t4, so we have to consider log of negative number
1936-
# This would be much simpler, but `return_sign` is available, so use it?
1937-
# t5 = sc.logsumexp([t3, t4 + np.pi*1j])
1938-
t5, sign = sc.logsumexp([t3, t4], b=[one, -one], axis=0, return_sign=True)
1939-
return sc.logsumexp([t1, t2 + t5 - np.log(a + b)], b=[one, -one*sign], axis=0)
1926+
with np.errstate(invalid='ignore', divide='ignore'):
1927+
log_y, m = np.log(x), u # compare against [1] Eq. 2
1928+
z = (log_y - m) / s
1929+
x1 = a * s - z
1930+
x2 = b * s + z
1931+
t1 = self._logPhi(z)
1932+
t2 = self._logphi(z)
1933+
t3 = (np.log(b) + self._logR(x1))
1934+
t4 = (np.log(a) + self._logR(x2))
1935+
t1, t2, t3, t4, one = np.broadcast_arrays(t1, t2, t3, t4, 1)
1936+
# t3 can be smaller than t4, so we have to consider log of negative number
1937+
# This would be much simpler, but `return_sign` is available, so use it?
1938+
# t5 = sc.logsumexp([t3, t4 + np.pi*1j])
1939+
t5, sign = sc.logsumexp([t3, t4], b=[one, -one], axis=0, return_sign=True)
1940+
temp = [t1, t2 + t5 - np.log(a + b)]
1941+
out = np.asarray(sc.logsumexp(temp, b=[one, -one*sign], axis=0))
1942+
out[x == 0] = -np.inf
1943+
return out[()]
19401944

19411945
def _logsf(self, x, u, s, a, b):
19421946
return _log1mexp(self._logcdf(x, u, s, a, b))

0 commit comments

Comments
 (0)