Logpdf support

[mborland]

No description provided.

[jzmaddock]

Thanks Matt, that's more or less along the lines I was thinking too.

I think we can safely add a default implementation:

template <class Distribution>
inline typename Distribution::value_type logpdf(const Distribution& dist, const typename Distribution::value_type& x)
{
  using std::log;
   return log(pdf(dist, x));
}

For those distributions that do not have exponential-like pdf's. The Arcsine is actually a good case in point, do we need an actual implementation here, or is the default good enough? I can see some arguments either way: if this->hi is super large then (x-lo)(hi-x) could overflow (but no one has ever complained), on the other hand, evaluating by logs as you do, may lead to cancellation error when the pdf is close to 1. And I doubt there is much performance difference between them.

With regard to testing... we already have decent tests for the pdf's, so I suspect that a basic sanity check against the pdf may be good enough in thee cases?

[mborland]

@jzmaddock This is clean and good for review. Logpdf has been specialized when there are exponentials in the pdf and there rest will default to the naive log(pdf). The only place where I had to adjust tolerances was in the test of long doubles on the normal distribution. The naive, specialized, nor hard coding the result from WolframAlpha got me into the existing tolerance band.

[jzmaddock]

There's a few comments above, I haven't checked the details of the formulas, and am replying on your tests.

We should also add logpdf to the DistributionConcept checks in test/compile_test/test_compile_result.hpp

Otherwise, looks good, much thanks for taking this on!

[NAThompson]

One thing I worry about these tests is that if they fail at some later date, it's not clear why (say) 5.7630258931329868780772138043668005779060097243996L is the correct value. Is there a property that can be tested instead?

[mborland]

@jzmaddock Your review comments have been addressed; I apologize for the delay. I also ran the sonar lint and clang-tidy combs through this as suggested by @ckormanyos.

[jzmaddock]

Hi Matt, thanks again for this, I just spotted a couple of location where you're still using log(tgamma(x)) rather than lgamma(x) but other than that it all looks good to go to me.

[NAThompson]

@HDembinski: Looks like this is basically finished. Let us know if you hit any bugs in the usage.

[mborland]

CI is clean so merging and closing linked issue.

[HDembinski]

Thanks for this. It looked fine, but perhaps you can remove some of the now duplicated code between pdf and logpdf. For some distributions (poisson, normal, exponential), the pdf can be computed from the logpdf.

[jzmaddock]

Thanks for this. It looked fine, but perhaps you can remove some of the now duplicated code between pdf and logpdf. For some distributions (poisson, normal, exponential), the pdf can be computed from the logpdf.

Technically yes, but would tend to loose precision due to cancellation error in the logpdf calculation?

[NAThompson]

Technically yes, but would tend to loose precision due to cancellation error in the logpdf calculation?

Yeah I'd like to see an ulps plot before going for that . . . also log and exp are not particularly cheap.

[HDembinski]

Fair enough, but you can do it at least for the exponential distribution, no? In case of the normal, you may be right. In case of Poisson, I don't know.

[NAThompson]

Fair enough, but you can do it at least for the exponential distribution, no?

For the rate λ, we'd have log(pdf(x)) = log(λ) - λx, and then exp(log(pdf)) = exp(log(λ) - λx). That gives two transcendental function calls over 1 using the code-duplicated λexp(-λx). Losing (or gaining!) precision on the exp(log(λ)) also seems possible, although generally you have half-ulp guarantees with those functions. Maybe I'm overthinking it but I just feel like an ulps plot would help convince me.