Add Kumaraswamy distribution (#2994)

[paulkernfeld]

Please let me know what additional testing you'd like to see for this. It looks like many of the distribution tests compare the PyMC3 implementation to a scipy implementation, but I don't think that scipy.stats implements the Kumaraswamy distribution.

close #2994

[fonnesbeck]

Thanks for the contribution. This looks really good. You have a test for the random method, but not for the logp (which belongs in test_distributions.py).

[paulkernfeld]

Great! I have added a test for the log PDF.

I noticed that the Travis build for this branch was failing due to TestNUTSLKJCholeskyCov, which looks to be fixed in #2997.

[junpenglao]

LGTM - could you please add a line to the RELEASE-NOTES?

[ColCarroll]

Actually, this looks like it fails for A Good Reason, but I can't figure out why. I would guess that it fails whenever a=0 , since then the mean uses 1 / a. Frustratingly, that just convinces me that test_beta should fail in the exact same case...

    def test_kumaraswamy(self):
        # Scipy does not have a built-in Kumaraswamy pdf
        def scipy_log_pdf(value, a, b):
            return np.log(a) + np.log(b) + (a - 1) * np.log(value) + (b - 1) * np.log(1 - value ** a)
>       self.pymc3_matches_scipy(Kumaraswamy, Unit, {'a': Rplus, 'b': Rplus}, scipy_log_pdf)
pymc3/tests/test_distributions.py:588: 
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 
pymc3/tests/test_distributions.py:430: in pymc3_matches_scipy
    model = build_model(pymc3_dist, domain, paramdomains, extra_args)
pymc3/tests/test_distributions.py:157: in build_model
    distfam('value', shape=valuedomain.shape, transform=None, **vals)
pymc3/distributions/distribution.py:37: in __new__
    return model.Var(name, dist, data, total_size)
pymc3/model.py:801: in Var
    total_size=total_size, model=self)
pymc3/model.py:1200: in __init__
    distribution.shape, distribution.dtype) * distribution.default()
pymc3/distributions/distribution.py:62: in default
    return np.asarray(self.get_test_val(self.testval, self.defaults), self.dtype)
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 
self = <pymc3.distributions.continuous.Kumaraswamy object at 0x7f3b5ba367d0>
val = None, defaults = ('median', 'mean', 'mode')
    def get_test_val(self, val, defaults):
        if val is None:
            for v in defaults:
                if hasattr(self, v) and np.all(np.isfinite(self.getattr_value(v))):
                    return self.getattr_value(v)
        else:
            return self.getattr_value(val)
    
        if val is None:
            raise AttributeError("%s has no finite default value to use, "
                                 "checked: %s. Pass testval argument or "
                                 "adjust so value is finite."
>                                % (self, str(defaults)))
E           AttributeError: <pymc3.distributions.continuous.Kumaraswamy object at 0x7f3b5ba367d0> has no finite default value to use, checked: ('median', 'mean', 'mode'). Pass testval argument or adjust so value is finite.

[paulkernfeld]

To fix the tests, I tried to change the domain of the parameters to exclude 0. Unfortunately, running py.test -k umaraswamy did not cause any test failures on my local OSX system. I don't know why that would be. Is this a valid way to run the relevant tests?

[ColCarroll]

I'm not familiar with the -k flag (it looks like it accepts a python expression?), but py.test -v pymc3/tests/test_distributions.py::TestMatchesScipy::test_kumaraswamy should work

[gBokiau]

Maybe move this to where other domains are defined?
I would consider using e-5 as first value, and also use this domain for Beta tests.

[gBokiau]

… actually, the same holds for parameters of gamma, exponential, etc. so changing these tests to use this domain is perhaps better left for a subsequent merge.

[gBokiau]

@paulkernfeld actually the first and last values of the domain are excluded by default, so it's really the 0.01 that is causing trouble in the tests. And indeed, with a=.01 and b=.01, mean=inf.

Hence, discard my review.

[ColCarroll]

Wow! is that why test_beta passes?

[gBokiau]

Further inspection… this is actually a numerical overflow issue. Computing exp(logmean) works, ie:

logmean = tt.log(b) + tt.gammaln(1. + 1./b) + tt.gammaln(b) - tt.gammaln(1 + 1 / a + b)
mean = tt.exp(logmean)

[paulkernfeld]

@gBokiau great insight! I have pushed a patch that computes the mean and variance using more logs. How did you figure out that this was the problem?

@ColCarroll I have reverted my change that added Rpositive. Perhaps it would be worth it to make a separate PR that runs tests for the entire domain of each variable, although that would require separate positive and nonnegative domains, which would probably have ripple effects throughout the code base.

[gBokiau]

@paulkernfeld it didn't make sense for it to have infinite mean since its a unit measure, then noticed the gamma functions, esp. w/ reciprocals. Gamma functions are very prone to n.o., Weibull has the same affliction.

Domains absolutely need a rewrite (#2966)

[junpenglao]

Thanks @paulkernfeld! And congrats on the first PR ;-)