# Do we need more unit tests for distributions?

**URL:** <https://discourse.mc-stan.org/t/do-we-need-more-unit-tests-for-distributions/12482>\
**Category:** Developers\
**Tags:** math\
**Created:** [January 2, 2020, 4:26pm UTC](https://discourse.mc-stan.org/t/do-we-need-more-unit-tests-for-distributions/12482 "2020-01-02T16:26:20Z")\
**Posts on this page:** 9\
**Page:** 1

<div class="post-metadata">

**Author:** ![martinmodrak](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/martinmodrak/32/133_2.png) [@martinmodrak](https://discourse.mc-stan.org/u/martinmodrak)\
**Post date:** [January 2, 2020, 4:26pm UTC](https://discourse.mc-stan.org/t/do-we-need-more-unit-tests-for-distributions/12482/1 "2020-01-02T16:26:20Z")

</div>

Recently a number of both subtle and not-so-subtle issues were reported (partly by me) for multiple probability distributions:

[https://github.com/stan-dev/math/issues/1495](https://github.com/stan-dev/math/issues/1495), [https://github.com/stan-dev/math/issues/1496](https://github.com/stan-dev/math/issues/1496), [https://github.com/stan-dev/math/issues/1498](https://github.com/stan-dev/math/issues/1498), [https://github.com/stan-dev/math/issues/1531](https://github.com/stan-dev/math/issues/1531)

This involves even the `neg_binomial_2`, which I thought is very “basic” and is presumably used in thousands of deployed models. Now most of the issues I’ve seen would rarely completely invalidate a model, but I can see some of them causing divergences here and there (many of the issues could manifest as large discontinuities in density for extreme parameters).

This might indicate we need more/better unit tests for distributions. As an example, when fixing the initial issue I found in `neg_binomial_2`, I discovered 4 additional issues and wrote or overhauled 7 tests to achieve some confidence that all code paths work (see [https://github.com/stan-dev/math/pull/1497](https://github.com/stan-dev/math/pull/1497) for more details).

We rarely test the distribution for anything but “normal” values of parameters - extremes (large, close to 0, …) are tested at most for not throwing anything or not returning nonsense, but not for accuracy. Many tests actually involve only one parameter combination. Some distributions (e.g. inv. gamma) have no tests for derivatives.

Some relatively simple ways to strengthen the tests are IMHO:

- Iterate over a grid of parameter values spanning several order of magnitudes for all extant tests
- Compute values for both the result and the derivatives in Mathematica and test against those
- Whenever a discrete choice is made (as when delegating from Neg. binomial to Poisson), test continuity at the boundary
- Test derivatives from autodiff against derivatives computed with complex step.

I’ve tried to build a few tests for some distributions that are less “obvious”, and all of them pass, which is nice (see [https://github.com/stan-dev/math/pull/1575](https://github.com/stan-dev/math/pull/1575)). So maybe `neg_binomial_2` was an outlier with unexpectedly tricky implementation while everything else is mostly OK?

I also understand that adding the tests might be a big effort and not necessarily worth it, hence this is very much up to discussion.

---

<div class="post-metadata">

**Author:** ![martinmodrak](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/martinmodrak/32/133_2.png) [@martinmodrak](https://discourse.mc-stan.org/u/martinmodrak)\
**Post date:** [January 3, 2020, 9:17am UTC](https://discourse.mc-stan.org/t/do-we-need-more-unit-tests-for-distributions/12482/2 "2020-01-03T09:17:16Z")

</div>

OK, so I tweaked a bit the tests for `neg_binomial_lpmf` and I am quite convinced it too has issues with derivatives for extreme values: [https://github.com/stan-dev/math/issues/1578](https://github.com/stan-dev/math/issues/1578)

---

<div class="post-metadata">

**Author:** ![ahartikainen](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/ahartikainen/32/321_2.png) [@ahartikainen](https://discourse.mc-stan.org/u/ahartikainen)\
**Post date:** [January 3, 2020, 9:09pm UTC](https://discourse.mc-stan.org/t/do-we-need-more-unit-tests-for-distributions/12482/3 "2020-01-03T21:09:14Z")

</div>

> [@martinmodrak](#):
>
> I also understand that adding the tests might be a big effort and not necessarily worth it,

Maybe this could be one topic that GSoC (Google Summer of Code) student could do, at least part of the work.

---

<div class="post-metadata">

**Author:** ![Bob\_Carpenter](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bob_carpenter/32/9230_2.png) [@Bob\_Carpenter](https://discourse.mc-stan.org/u/Bob_Carpenter)\
**Post date:** [January 4, 2020, 11:36pm UTC](https://discourse.mc-stan.org/t/do-we-need-more-unit-tests-for-distributions/12482/4 "2020-01-04T23:36:08Z")

</div>

I agree that we need to test the distributions better, including for derivatives.

I think the best course forward is to expand the finite differences test package. We can test against known derivatives in some cases, but it’s combinatorially prohibitive to do it by hand for first, second, and third-order derivatives using all of our autodiff options.

We need to think more about the possible accuracies we can get. It’s challenging to test over several orders of magnitude.

> [@martinmodrak](#):
>
> maybe `neg_binomial_2` was an outlier with unexpectedly tricky implementation while everything else is mostly OK?

We won’t know until we test :-)

> [@martinmodrak](#):
>
> I also understand that adding the tests might be a big effort and not necessarily worth it, hence this is very much up to discussion.

Very much worth it. I’d make testing a higher priority, but it’s hard to get funding for testing and hard to find people that want to do it.

> [@ahartikainen](#):
>
> Maybe this could be one topic that GSoC (Google Summer of Code) student could do,

Sure—anyone could do this. Is anyone applying for a GSoC student?

---

<div class="post-metadata">

**Author:** ![ahartikainen](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/ahartikainen/32/321_2.png) [@ahartikainen](https://discourse.mc-stan.org/u/ahartikainen)\
**Post date:** [January 5, 2020, 12:09am UTC](https://discourse.mc-stan.org/t/do-we-need-more-unit-tests-for-distributions/12482/5 "2020-01-05T00:09:32Z")

</div>

> [@Bob\_Carpenter](#):
>
> Sure—anyone could do this. Is anyone applying for a GSoC student?

Stan project could apply one under NumFocus umbrella.

> **[Google Summer of Code](https://summerofcode.withgoogle.com/)**
>
> Google Summer of Code is a global program focused on bringing more developers into open source software development.

I know Stan tried this once before, but didn’t get any students, but maybe this year it could work.

Last summer ArviZ had 2 students (I was mentoring one) and the experience was positive.

---

<div class="post-metadata">

**Author:** ![martinmodrak](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/martinmodrak/32/133_2.png) [@martinmodrak](https://discourse.mc-stan.org/u/martinmodrak)\
**Post date:** [January 6, 2020, 12:13pm UTC](https://discourse.mc-stan.org/t/do-we-need-more-unit-tests-for-distributions/12482/6 "2020-01-06T12:13:41Z")

</div>

> [@Bob\_Carpenter](#):
>
> I think the best course forward is to expand the finite differences test package.

In my work with the tests for neg. binomial, testing against finite differences was very hard to setup to both provide strong test while not producing false positives. Using `boost::math::differentiation::complex_step_derivative` allowed me to keep the relative acuraccy for the test at `1e-8` (except for few extreme cases) so I would suggest moving to complex step whenever possible. The downside is that this requires the densities to work with `std::complex<double>` substituted for all differentiable arguments. This may not work for some functions - notably `lgamma(std::complex<double>)` is not available, but workarounds are relatively easy to write.

---

<div class="post-metadata">

**Author:** ![Bob\_Carpenter](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bob_carpenter/32/9230_2.png) [@Bob\_Carpenter](https://discourse.mc-stan.org/u/Bob_Carpenter)\
**Post date:** [January 15, 2020, 2:01am UTC](https://discourse.mc-stan.org/t/do-we-need-more-unit-tests-for-distributions/12482/7 "2020-01-15T02:01:10Z")

</div>

> [@ahartikainen](#):
>
> Stan project could apply one under NumFocus umbrella.

We can. Someone told me that they put in a project to NumFOCUS but NumFOCUS dropped the ball and never listed it. I’m not 100% sure what the story was. So if you do go this way, make double sure it gets listed!

> [@martinmodrak](#):
>
> I would suggest moving to complex step whenever possible.

How do we know when it’s possible? A lot of our distributions use `lgamma` or even more complicated functions.

> [@martinmodrak](#):
>
> workarounds are relatively easy to write.

Can you explain more?

---

<div class="post-metadata">

**Author:** ![martinmodrak](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/martinmodrak/32/133_2.png) [@martinmodrak](https://discourse.mc-stan.org/u/martinmodrak)\
**Post date:** [January 15, 2020, 1:16pm UTC](https://discourse.mc-stan.org/t/do-we-need-more-unit-tests-for-distributions/12482/8 "2020-01-15T13:16:27Z")

</div>

> [@Bob\_Carpenter](#):
>
> How do we know when it’s possible? A lot of our distributions use `lgamma` or even more complicated functions.  
> Can you explain more?

Sorry, I was unnecessarily mysterious.

There is a bunch of functions that don’t admit `std::complex<double>` as argument by default, but I’ve succesfully made them work enough to compute complex step derivative of functions involving them:

- `log1p`, `expm1` etc. those have quite short code and usually fail due to lines like `if(x < 1e-4) { ... }`, providing a complex specialization often means just copying the code and replacing comparison of values with comparison of absolute value (e.g. `if(abs(x) < 1e-4) { ... }`) - the underlying approximations motivating many of the implementations are usually valid in complex plane as well.
- `lgamma()` and other complex functions can be approximated via Taylor expansion around `x.imag() == 0`. This might not work when the argument to the function is already heavily transformed (e.g. `lgamma(exp(15 * x))`, but when the argument is “close” to the original (e.g. `lgamma(x + 1)`), the imaginary part will be the small “complex step” and the approximation holds. For `lgamma` this gives:

```no-highlight
auto lgamma_c_approx = [](const std::complex<double>& x) {
  return std::complex<double>(lgamma(x.real()),
                              x.imag() * boost::math::digamma(x.real()));
  };

```

There are two possible approaches to computing complex step derivative with such custom implementations:

- build specializations of necessary library functions admitting `std::complex`, making the current implementation of the function work with complex step. This could be a test-only library as we probably don’t want to invest in vouching for 100% correctness of the implementations. On the other hand, as there is work on [`std::complex<var>`](https://discourse.mc-stan.org/t/help-testing-release-candidate-for-std-complex/12314), such functions might be necessary anyway. For `lgamma` (and possibly others) there are complex implementations we could easily use, e.g. the [one on Wikipedia](https://en.wikipedia.org/wiki/Lanczos_approximation#Simple_implementation)
- write a specialized test-only version of the function to be tested (or it’s simplest formula) using specialized functions - this is what I currently do for the [`neg_binomial_2_lpmf` test against complex step](https://github.com/stan-dev/math/blob/fbe958d633c0701b91edc51405626bc253371866/test/unit/math/rev/prob/neg_binomial_2_test.cpp#L440)

Hope that is more specific and actionable.

---

<div class="post-metadata">

**Author:** ![Bob\_Carpenter](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bob_carpenter/32/9230_2.png) [@Bob\_Carpenter](https://discourse.mc-stan.org/u/Bob_Carpenter)\
**Post date:** [January 23, 2020, 7:00pm UTC](https://discourse.mc-stan.org/t/do-we-need-more-unit-tests-for-distributions/12482/9 "2020-01-23T19:00:48Z")

</div>

> [@martinmodrak](#):
>
> Hope that is more specific and actionable.

Thanks! Yes, that’s clear.

I have [a branch in stan-math](https://github.com/stan-dev/math/tree/feature/0123-complex-var) that’s going to fold complex number support into Stan math for autodiff vars and matrices.

It’d be great if we could start expanding our library to handle complex inputs, even if we just have to define everything templated to start and autodiff through it. The branch above will provide support for autodiff for every function defined for `std::complex` in the standard library.