From 38fa415c15b0c0469fbd9cad3a3b9ba974fc8733 Mon Sep 17 00:00:00 2001 From: Dominik Date: Thu, 5 Sep 2019 15:33:57 +0200 Subject: [PATCH 1/4] remove eps from numpy StickBreaking #3577 --- pymc3/distributions/transforms.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pymc3/distributions/transforms.py b/pymc3/distributions/transforms.py index 17c400d2af..cdc7219b6c 100644 --- a/pymc3/distributions/transforms.py +++ b/pymc3/distributions/transforms.py @@ -485,7 +485,7 @@ def backward_val(self, y_): Km1 = y.shape[0] k = np.arange(Km1)[(slice(None),) + (None,) * (y.ndim - 1)] eq_share = nplogit(1.0 / (Km1 + 1 - k).astype(str(y_.dtype))) - z = expit(y + eq_share, self.eps) + z = expit(y + eq_share) yl = np.concatenate([z, np.ones(y[:1].shape)]) yu = np.concatenate([np.ones(y[:1].shape), 1 - z]) S = np.cumprod(yu, 0) From 434068a8c6e21b30613000f5c37582d5b8f86fe5 Mon Sep 17 00:00:00 2001 From: Dominik Date: Fri, 6 Sep 2019 16:48:53 +0200 Subject: [PATCH 2/4] simplified StickBreaking --- pymc3/distributions/transforms.py | 82 +++++++++++-------------------- 1 file changed, 30 insertions(+), 52 deletions(-) diff --git a/pymc3/distributions/transforms.py b/pymc3/distributions/transforms.py index cdc7219b6c..5095449ed5 100644 --- a/pymc3/distributions/transforms.py +++ b/pymc3/distributions/transforms.py @@ -1,3 +1,5 @@ +import warnings + import theano import theano.tensor as tt @@ -89,7 +91,8 @@ def backward(self, z): raise NotImplementedError def jacobian_det(self, x): - """Calculates logarithm of the absolute value of the Jacobian determinant for input `x`. + """Calculates logarithm of the absolute value of the Jacobian determinant + of the backward transformation for input `x`. Parameters ---------- @@ -431,77 +434,52 @@ def jacobian_det(self, x): class StickBreaking(Transform): """ Transforms K - 1 dimensional simplex space (k values in [0,1] and that sum to 1) to a K - 1 vector of real values. - Primarily borrowed from the STAN implementation. - - Parameters - ---------- - eps : float, positive value - A small value for numerical stability in invlogit. """ name = "stickbreaking" - def __init__(self, eps=floatX(np.finfo(theano.config.floatX).eps)): - self.eps = eps + def __init__(self, eps=None): + if eps is not None: + warnings.warn("The argument `eps` is depricated and will not be used.", + DeprecationWarning) def forward(self, x_): x = x_.T - # reverse cumsum - x0 = x[:-1] - s = tt.extra_ops.cumsum(x0[::-1], 0)[::-1] + x[-1] - z = x0 / s - Km1 = x.shape[0] - 1 - k = tt.arange(Km1)[(slice(None),) + (None,) * (x.ndim - 1)] - eq_share = logit(1.0 / (Km1 + 1 - k).astype(str(x_.dtype))) - y = logit(z) - eq_share + n = x.shape[0] + lx = tt.log(x) + shift = tt.sum(lx, 0, keepdims=True) / n + y = lx[:-1] - shift return floatX(y.T) - def forward_val(self, x_, point=None): + def forward_val(self, x_): x = x_.T - # reverse cumsum - x0 = x[:-1] - s = np.cumsum(x0[::-1], 0)[::-1] + x[-1] - z = x0 / s - Km1 = x.shape[0] - 1 - k = np.arange(Km1)[(slice(None),) + (None,) * (x.ndim - 1)] - eq_share = nplogit(1.0 / (Km1 + 1 - k).astype(str(x_.dtype))) - y = nplogit(z) - eq_share + n = x.shape[0] + lx = np.log(x) + shift = np.sum(lx, 0, keepdims=True) / n + y = lx[:-1] - shift return floatX(y.T) def backward(self, y_): y = y_.T - Km1 = y.shape[0] - k = tt.arange(Km1)[(slice(None),) + (None,) * (y.ndim - 1)] - eq_share = logit(1.0 / (Km1 + 1 - k).astype(str(y_.dtype))) - z = invlogit(y + eq_share, self.eps) - yl = tt.concatenate([z, tt.ones(y[:1].shape)]) - yu = tt.concatenate([tt.ones(y[:1].shape), 1 - z]) - S = tt.extra_ops.cumprod(yu, 0) - x = S * yl + y = tt.concatenate([y, -tt.sum(y, 0, keepdims=True)]) + # "softmax" with vector support and no deprication warning: + e_y = tt.exp(y - tt.max(y, 0, keepdims=True)) + x = e_y / tt.sum(e_y, 0, keepdims=True) return floatX(x.T) def backward_val(self, y_): y = y_.T - Km1 = y.shape[0] - k = np.arange(Km1)[(slice(None),) + (None,) * (y.ndim - 1)] - eq_share = nplogit(1.0 / (Km1 + 1 - k).astype(str(y_.dtype))) - z = expit(y + eq_share) - yl = np.concatenate([z, np.ones(y[:1].shape)]) - yu = np.concatenate([np.ones(y[:1].shape), 1 - z]) - S = np.cumprod(yu, 0) - x = S * yl + y = np.append(y, -np.sum(y, 0, keepdims=True)) + x = np.exp(y)/np.sum(np.exp(y), 0, keepdims=True) return floatX(x.T) - def jacobian_det(self, y_): - y = y_.T - Km1 = y.shape[0] - k = tt.arange(Km1)[(slice(None),) + (None,) * (y.ndim - 1)] - eq_share = logit(1.0 / (Km1 + 1 - k).astype(str(y_.dtype))) - yl = y + eq_share - yu = tt.concatenate([tt.ones(y[:1].shape), 1 - invlogit(yl, self.eps)]) - S = tt.extra_ops.cumprod(yu, 0) - return tt.sum(tt.log(S[:-1]) - tt.log1p(tt.exp(yl)) - tt.log1p(tt.exp(-yl)), 0).T - + def jacobian_det(self, x_): + x = x_.T + n = x.shape[0] + # stable according to: http://deeplearning.net/software/theano_versions/0.9.X/NEWS.html + lse = tt.log(tt.sum(tt.exp(x), 0, keepdims=True)) + d = tt.log(n) - (n*lse) + return d.T stick_breaking = StickBreaking() From 67ece13de83a62325a7027870c24bd5ce1b9f3c1 Mon Sep 17 00:00:00 2001 From: Dominik Date: Fri, 6 Sep 2019 18:30:38 +0200 Subject: [PATCH 3/4] fix StickBreaking jacobian --- pymc3/distributions/transforms.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/pymc3/distributions/transforms.py b/pymc3/distributions/transforms.py index 5095449ed5..663a022846 100644 --- a/pymc3/distributions/transforms.py +++ b/pymc3/distributions/transforms.py @@ -476,9 +476,9 @@ def backward_val(self, y_): def jacobian_det(self, x_): x = x_.T n = x.shape[0] - # stable according to: http://deeplearning.net/software/theano_versions/0.9.X/NEWS.html - lse = tt.log(tt.sum(tt.exp(x), 0, keepdims=True)) - d = tt.log(n) - (n*lse) + sx = tt.sum(x, 0, keepdims=True) + r = tt.log(floatX(1) + tt.sum(tt.exp(x + sx), 0,keepdims=True)) + d = tt.log(n) + (n*sx) - (n*r) return d.T stick_breaking = StickBreaking() From d82ded70339075173f39d956b87fd9037e756d2d Mon Sep 17 00:00:00 2001 From: Dominik Date: Fri, 6 Sep 2019 18:45:21 +0200 Subject: [PATCH 4/4] use stable logsumexp in StickBreaking --- pymc3/distributions/transforms.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/pymc3/distributions/transforms.py b/pymc3/distributions/transforms.py index 663a022846..ce6da7964d 100644 --- a/pymc3/distributions/transforms.py +++ b/pymc3/distributions/transforms.py @@ -469,7 +469,7 @@ def backward(self, y_): def backward_val(self, y_): y = y_.T - y = np.append(y, -np.sum(y, 0, keepdims=True)) + y = np.concatenate([y, -np.sum(y, 0, keepdims=True)]) x = np.exp(y)/np.sum(np.exp(y), 0, keepdims=True) return floatX(x.T) @@ -477,8 +477,10 @@ def jacobian_det(self, x_): x = x_.T n = x.shape[0] sx = tt.sum(x, 0, keepdims=True) - r = tt.log(floatX(1) + tt.sum(tt.exp(x + sx), 0,keepdims=True)) - d = tt.log(n) + (n*sx) - (n*r) + r = tt.concatenate([x+sx, tt.zeros(sx.shape)]) + # stable according to: http://deeplearning.net/software/theano_versions/0.9.X/NEWS.html + sr = tt.log(tt.sum(tt.exp(r), 0, keepdims=True)) + d = tt.log(n) + (n*sx) - (n*sr) return d.T stick_breaking = StickBreaking()