[avehtari]

• Closes poisson_log cdfs #90: poisson_log cdfs, so that a poisson_log likelihood can be truncated (is this the record for the oldest issue to be closed?).
• Partly addresses Implement missing cdf/lcdf/lccdf functions for distributions #2657 (the poisson_log entry of its list).

This PR depends on #3408, which adds the internal function that it uses.

An earlier attempt, #471 (2016), computed poisson_cdf(n, exp(alpha)) and was closed by its author before review.

This PR was Claude assisted

Summary

poisson_log had only lpmf and rng. Without lcdf and lccdf, Stan cannot truncate it, and users who model with the log rate had to write poisson_lcdf(n | exp(alpha)), which loses digits and gradients when exp(alpha) is large or small. This PR adds the three functions:

• poisson_log_lcdf(n | alpha) = log Q(n + 1, exp(alpha))
• poisson_log_lccdf(n | alpha) = log P(n + 1, exp(alpha))
• poisson_log_cdf(n | alpha) = Q(n + 1, exp(alpha))

Here P and Q are the regularized lower and upper incomplete gamma functions. Each function evaluates log P or log Q with internal::log_gamma_tail from #3408, which works in log space, so the value and the gradient stay finite where the probability underflows. The derivative with respect to alpha is the derivative with respect to the rate times exp(alpha). The probability-scale cdf sums the logs and multiplies by exp(sum) at the end, as the other cdfs in #3408 do.

The argument checks follow poisson_log_lpmf: alpha must not be NaN, and alpha = -inf and alpha = +inf are allowed.

Values at the edges of the support, element by element (an element with probability 0 makes the whole result -inf or 0; an element with probability 1 contributes nothing):

case	lcdf	lccdf	cdf
n < 0	-inf	0	0
n = INT_MAX (an infinite count, as in neg_binomial)	0	-inf	1
exp(alpha) = 0 (alpha = -inf, or below about -745)	0	-inf	1
exp(alpha) = inf (alpha = +inf, or above about 709.8)	-inf	0	0

Accuracy

Against mpmath at 60 digits at 54 points (the three functions at 18 points (n + 1, exp(alpha)) of the incomplete gamma function, from the deep lower tail to the deep upper tail and with counts up to 999), computed by two routes that agree to 1e-25: the worst relative error is 1.1e-12 for values and 1.5e-12 for gradients. Most of the value error comes from the rounding of exp(alpha): near the mean of a large count the value changes by thousands of times the relative change of the rate. The functions agree with poisson_*(n, exp(alpha)) to all printed digits, also in fvar<fvar<double>>.

Testing

• test/unit/math/rev/prob/poisson_log_cdfs_reference.hpp and poisson_log_cdfs_test.cpp: the 54 reference rows with var (value to 1e-11, gradient to 1e-9 relative).
• test/unit/math/prim/prob/poisson_log_cdfs_test.cpp: agreement with poisson_lcdf, poisson_lccdf and poisson_cdf at exp(alpha) for vectors, the edge values of the table above, and the NaN check.
• test/unit/math/mix/prob/poisson_log_cdfs_test.cpp: expect_ad for the three functions at two points each.
• test/prob/poisson/poisson_log_{cdf,cdf_log,ccdf_log}_test.hpp: the generated distribution tests, with the points of the poisson cdf tests at alpha = log(lambda). In the ccdf_log test the point n = 82, lambda = 42 is replaced by n = 40: the test's reference function computes log(1 - cdf), which keeps only about 8 digits where the ccdf is 1.5e-8, and its third derivative then differed from the exact value by more than the test tolerance (this PR agrees with mpmath there to 11 digits).

The new and existing poisson unit suites pass, and all 18 generated test/prob/poisson cdf programs pass. Header checks, make test-math-dependencies, cpplint and clang-format 10 pass, and make doxygen with doxygen 1.9.7 gives 0 warnings. Tested on develop at e50e79f188 with #3408 applied, g++ 11.5.

The stanc3 PR stan-dev/stanc3#TBD adds the signatures.

Release notes

Add poisson_log_cdf, poisson_log_lcdf and poisson_log_lccdf.

Checklist

• Copyright holder: Aalto University

  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

  • new unit tests pass (to run, use: ./runTests.py test/unit)
  • header checks pass, (make test-headers)
  • dependencies checks pass, (make test-math-dependencies)
  • docs build, (make doxygen)
  • code passes the built in C++ standards checks (make cpplint)

• the code is written in idiomatic C++ and changes are documented in the doxygen

• the new changes are tested