Best way to check implementation of density, distribution function and random generation What is the best way to check if implementation of density, distribution function, quantile function and random generation for some distribution are correct? For example, base R lacks Laplace distribution. Let's say we implement it. We can use random generation function and check if generated values are "similar" to density and distribution functions, but this is quite idem per idem kind of check. We know that if $X,Y \sim \mathrm{Exp}(\lambda)$ then $X-Y \sim \mathrm{Laplace}(0, \lambda^{-1})$, so we can generate exponential distributed values and use them to imitate Laplace distribution, however how do we know that we are close enough? We can check also if $X = F_X^{-1}(F_X(X))$, but what more? What other issues need to be checked (e.g. underflow or overflow) and what are the best practices?
 A: After creating my own R package implementing a number of probability distributions, I have some thoughts about verifying correctness of the functions. For a nice starting point one could check the tests implemented in base R for testing the default distribution, that can be found in tests/d-p-q-r-tests.R and tests/p-r-random-tests.R files.
There is a number of formal properties that need to be met and checks that should be made:

*

*It is a good idea to plot everything: functions, the random draws, empirical cumulative distributions against theoretical cumulative distribution functions etc. There is a number of plots that can and should be considered. This tremendously helps in finding bugs.


*Check if $0 \le f(x) \le 1$ for discrete random variables and $0 \le f(x) < \infty$ for continuous random variables.


*For discrete distributions $\left|1 - (\sum_x f(x))\right| \le \varepsilon$ for $x$'s in the support of $f$ ranging from some arbitrary small, up to arbitrary large value (for some arbitrary small $\varepsilon$) and $F(k) = \sum_k f(k)$.


*Check if $0 \le F(x) \le 1$.


*Check if $F(-\infty) = 0$ and $F(\infty) = 1$.


*Check if $x = F^{-1}(F(x))$ or $p = F(F^{-1}(p))$


*As noted by Xi'an, $F(X)$ should be uniformly distributed. Remember not to be very strict about uniformity of random draws from discrete distributions.


*Moreover, in tests/p-r-random-tests.R R implements test based on an inequality of
Massart:
$$
\Pr\left( \sup |\hat F_n(x) - F(x)| > \lambda \right) \le 2\exp(-2n\lambda^2)
  $$
where $\hat F_n(x)$ is the empirical distribution function, that can be used to compare the randomly generated samples to the cumulative distribution function. When using R, the code can be simply copy-and-pasted from R's source and re-used for testing.


*Check the boundary conditions, i.e. $x \in \{-\infty, 0, \infty\}$ especially since some distributions  have discontinuities on $0$ that may be easily overseen (e.g. slash distribution). For bounded distribution it is important to check what happens on bounds and if $f(x) = 0$ for $x$ outside the support of $f$.


*Run checks on a wide range of parameter values (valid and invalid, e.g. against negative scale parameters).


*It is good to have a deeper thought about handling invalid parameter values, missing data, NaN's etc. For example, base R propagates missing values and NaN's, it returns NaN's and throws warnings NaNs produced for invalid parameter values in the d/p/q functions and returns NA's and throws NAs produced warnings in r functions, etc.


*Checking smoothness of the functions for very small (1e-13, 1e-14, 1e-15, 1e-16, ...), or very big values, may help you to diagnose the problems with numerical precision.
Some helpful hints are also given on slides  Software for Distributions in R by David Scott, Diethelm Wurtz and Christine Dong.
Moreover, it is important not to make equality checks when dealing with non-integers, since due to numerical precision, they will never be passed. Recall the advice given in Writing R Extensions document:

Only test the accuracy of results if you have done a formal error
analysis. Things such as checking that probabilities numerically sum
to one are silly: numerical tests should always have a tolerance. That
the tests on your platform achieve a particular tolerance says little
about other platforms. R is configured by default to make use of long
doubles where available, but they may not be available or be too slow
for routine use. Most R platforms use ‘ix86’ or ‘x86_64’ CPUs: these
use extended precision registers on some but not all of their FPU
instructions. Thus the achieved precision can depend on the compiler
version and optimization flags—our experience is that 32-bit builds
tend to be less precise than 64-bit ones. But not all platforms use
those CPUs, and not all81 which use them configure them to allow the
use of extended precision. In particular, ARM CPUs do not (currently)
have extended precision nor long doubles, and long double was 64-bit
on HP/PA Linux.
If you must try to establish a tolerance empirically, configure and
build R with --disable-long-double and use appropriate compiler flags
(such as -ffloat-store and -fexcess-precision=standard for gcc,
depending on the CPU type82) to mitigate the effects of
extended-precision calculations.
Tests which involve random inputs or non-deterministic algorithms
should normally set a seed or be tested for many seeds.

