Generalized normal distribution #3157
Conversation
|
Thanks for this! Can you remove the changes to the |
Jenkins Console Log Machine informationNo LSB modules are available. Distributor ID: Ubuntu Description: Ubuntu 20.04.3 LTS Release: 20.04 Codename: focalCPU: G++: Clang: |
|
I discarded any changes to the existing distributions |
Jenkins Console Log Machine informationNo LSB modules are available. Distributor ID: Ubuntu Description: Ubuntu 20.04.3 LTS Release: 20.04 Codename: focalCPU: G++: Clang: |
|
Sorry for a broken PR, I believed the tests pass, but apparently I made some mistake when running the tests locally. Would appreciate any help deciphering the error messages, it doesn't seem very understandable what the problem is exactly to me. But apparently it somehow has to do with the distribution being possibly non-differentiable at mu=y. The suspect line is:
In case when mu=y, we get:
which I think should evaluate to 0, and also when beta=1
should also evaluate to 0 I would expect? Is pow(0,0) actually undefined in STAN math? And what would be the recommended technique for special-case handling? |
No worries! It's often quite a miracle when these things pass first try. You can run the exact test that is given issues locally by running
It's defined (#3031), but it's not twice differentiable: #2993 (comment), which I believe is the problem here. This second comment is by @andrjohns who can probably also give you suggestions about what to do in this case. It's possible that what you're doing in the code is fine and you just need to have the test call it somewhere else where it will be better behaved |
|
OK, so it may be just a testing problem. But this relates to something else: The case when beta=1 where the tests are apparently breaking is maybe just a tip of the iceberg: when beta<1, I write in the code comments the partials are undefined when y=mu, because we then get |
| auto ops_partials | ||
| = make_partials_propagator(y_ref, mu_ref, alpha_ref, beta_ref); | ||
|
|
||
| if (!is_constant_all<T_y, T_loc>::value) { |
There was a problem hiding this comment.
For any if statements with known compile time values we can use if constexpr
| if (!is_constant_all<T_y, T_loc>::value) { | |
| if constexpr (!is_constant_all<T_y, T_loc>::value) { |
| const auto& rep_deriv | ||
| = to_ref_if < !is_constant<T_y>::value | ||
| && !is_constant<T_loc>::value | ||
| > (sign(diff) * beta_val * pow(scaled_abs_diff, beta_val - 1) | ||
| * inv_alpha); |
There was a problem hiding this comment.
Once you are in this block you know that T_y and T_loc are not constant
| const auto& rep_deriv | |
| = to_ref_if < !is_constant<T_y>::value | |
| && !is_constant<T_loc>::value | |
| > (sign(diff) * beta_val * pow(scaled_abs_diff, beta_val - 1) | |
| * inv_alpha); | |
| auto rep_deriv = eval(sign(diff) * beta_val * pow(scaled_abs_diff, beta_val - 1) * inv_alpha); |
| // note: The partial derivatives for y, mu are undefined when y == mu && | ||
| // beta < 1. The derivative limit as mu -> y goes: | ||
| // to 0 from both sides if beta > 1 (defined as 0) | ||
| // to +1/alpha from right but -1/alpha from left if beta == 1 (defined as | ||
| // 0, consistent with double_exponential_lpdf) to +∞ from right but -∞ | ||
| // from left as y -> mu if beta < 1 (undefined) |
There was a problem hiding this comment.
I had a hard time understanding your comment here. Is this the right interpretation? Either way can you clean it up a bit so the cut points and conditions are laid out a bit nicer
| // note: The partial derivatives for y, mu are undefined when y == mu && | |
| // beta < 1. The derivative limit as mu -> y goes: | |
| // to 0 from both sides if beta > 1 (defined as 0) | |
| // to +1/alpha from right but -1/alpha from left if beta == 1 (defined as | |
| // 0, consistent with double_exponential_lpdf) to +∞ from right but -∞ | |
| // from left as y -> mu if beta < 1 (undefined) | |
| // note: The partial derivatives for y, mu are undefined when | |
| // y == mu && beta < 1. | |
| // The derivative limit as mu -> y has the following cases: | |
| // beta > 1: 0 from both sides(defined as 0) | |
| // beta == 1: +1/alpha from right, but -1/alpha from left (defined as | |
| // 0, consistent with double_exponential_lpdf) to +∞ from right | |
| // beta < 1: -∞ from left as y -> mu (undefined) |
|
Some notes on the distribution (I made some mistakes in the #3133 comments but here should be corrected): We have the density: The parameters The log-likelihood is then: Wolfram Mathematica outputs the following partials: Using algebraic manipulations I change this to the form: This is where I believe the test errors are comming from when testing the case The this one has a also potential boundary problem when The I tend to believe that @WardBrian is right that this is perhaps a problem with testing the distribution at these boundaries, but idk what a good approach to this is. PDF plot from wikipedia for |
|
Incorporated changes requested by @SteveBronder although that doesn't solve the test errors. I think those are more like just performance+stylistic improvements. See the notes and discussion above. |
Jenkins Console Log Machine informationNo LSB modules are available. Distributor ID: Ubuntu Description: Ubuntu 20.04.3 LTS Release: 20.04 Codename: focalCPU: G++: Clang: |
Summary
This is my take on the request #3133 I created a while ago.
The Generalized normal distribution a.k.a. exponential power distribution generalizes the normal, Laplace (double exponential) and in limit also uniform distribution. This makes it an interesting prior choice as it generalizes the corresponding MAP regularization coefficient to L_p norm and can also be interesting to model noise in regression problems for basically the same reasons.
This is my first attempt at implementing a distribution to STAN, I tried to follow the guide https://mc-stan.org/math/md_doxygen_2contributor__help__pages_2adding__new__distributions.html but there are likely some problems, so I welcome any comments.
So far I implemented only the lpdf, if other files like lcdf etc are strictly required, please let me know.
Tests
I added one test file analogical to the one for normal distribution.
Side Effects
I included some very minor code polishing to the related normal and double exponential distributions (I used those as a starting point) which shouldn't do anything and can be discarded if undesired.
Release notes
Added generalized_normal_lpdf
Checklist
Copyright holder: Me
The copyright holder is typically you or your assignee, such as a university or company. By submitting this pull request, the copyright holder is agreeing to the license the submitted work under the following licenses:
- Code: BSD 3-clause (https://opensource.org/licenses/BSD-3-Clause)
- Documentation: CC-BY 4.0 (https://creativecommons.org/licenses/by/4.0/)
the basic tests are passing
./runTests.py test/unit)make test-headers)make test-math-dependencies)make doxygen)make cpplint)the code is written in idiomatic C++ and changes are documented in the doxygen
the new changes are tested