refactor kumaraswamy

[farhanreynaldo]

This PR closes a subset of #4686, which is Kumaraswamy distribution. But, I am not quite sure why the check_pymc_draws_match_reference test failed.

[ricardoV94]

Completely optional for this PR but would you be interested in implementing the logcdf method for this distribution?

Gauging from wikipedia it has a very simple form: https://en.m.wikipedia.org/wiki/Kumaraswamy_distribution

[ricardoV94]

Failing test is unrelated to this PR. Documented in #4661

[farhanreynaldo]

Completely optional for this PR but would you be interested in implementing the logcdf method for this distribution?

Gauging from wikipedia it has a very simple form: https://en.m.wikipedia.org/wiki/Kumaraswamy_distribution

My math is a little bit rusty, but does the following logcdf correct?

-b * at.log(1 - value ** a)

Because when I check against the following code, it doesn't match.

at.log(1 - (1 - value ** a) ** b)

[ricardoV94]

You are missing the first 1 - .... If you want to remove the exponent you would need log(1 - exp(b * log( 1 - value ** a)) but perhaps the outer log will suffice numerically. We can replace log(1 - x) by log1p(-x) which is a bit more stable.

That is log1p(- (1 - value ** a) ** b)

[farhanreynaldo]

You are missing the first 1 - .... If you want to remove the exponent you would need log(1 - exp(b * log( 1 - value ** a)) but perhaps the outer log will suffice numerically. We can replace log(1 - x) by log1p(-x) which is a bit more stable.

That is log1p(- (1 - value ** a) ** b)

I totally missed that! 😩 Thank you for clarifying the equation.

I should also add the following test to test_kumaraswamy, right?

def test_kumaraswamy(self):
        ...
        def scipy_log_cdf(value, a, b):
            return np.log1p(-(1 - value ** a) ** b)

        self.check_logp(Kumaraswamy, Unit, {"a": Rplus, "b": Rplus}, scipy_log_cdf)

[ricardoV94]

It should be check_logcdf, but other than that yes :)

[ricardoV94]

Value above one should give a logcdf of 0, so we need a switch: at.switch(value < 1, normal expression, 0), and to remove value <= 1 as a bound condition. We should still keep the a > 0 and b > 0 conditions as in the logp

[farhanreynaldo]

does tt equal to at? Because I assume tt is Theano tensor and it is not available.

[ricardoV94]

Yes of course. It was a mistake, I am still used to writing with theano's abbreviation :p

Fixed the comment

[ricardoV94]

Looks to be on the right track. I left you some comments. In the meantime you can add a release note mentioning the logcdf method was added to this distribution in the "New features": https://github.com/pymc-devs/pymc3/blob/v4/RELEASE-NOTES.md

[ricardoV94]

The failing test seems to be a precision issue. This means we may want to break the exponent as in your first suggestion. I'll push some changes in a second.

[ricardoV94]

I think it's ready! Thanks for the PR

[farhanreynaldo]

I think it's ready! Thanks for the PR

Thank you so much!