After [creating my own R package][1] 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`][2] and [`tests/p-r-random-tests.R`][3] files.

There is a number of formal properties that need to be met and checks that should be made:

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

 2. Check if $0 \le f(x) \le 1$ for discrete random variables and $f(x) \ge 0$ for continuous random variables.

 3. 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)$.

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

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

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

 7. 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.

 8. Moreover, [in `tests/p-r-random-tests.R` R implements][3] 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.

 9. It is important to 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][5]). For bounded distribution it is important to check what happens on bounds and if $f(x) = 0$ for $x$ outside the support of $f$.

 10. It is important to run checks on a wide range of parameter values (valid and invalid, e.g. against negative scale parameters).

 11. 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.

Some helpful hints are also given on slides  [*Software for Distributions in R*
by David Scott, Diethelm Wurtz and Christine Dong][7].

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*][6] 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.



  [1]: http://%20https://CRAN.R-project.org/package=extraDistr
  [2]: https://github.com/wch/r-source/blob/a84a1ba90a6db8c999d352d253b9088f4177e8ea/tests/d-p-q-r-tests.R
  [3]: https://github.com/wch/r-source/blob/a84a1ba90a6db8c999d352d253b9088f4177e8ea/tests/p-r-random-tests.R
  [4]: http://www.itl.nist.gov/div898/handbook/eda/section3/eda33.htm
  [5]: https://en.wikipedia.org/wiki/Slash_distribution
  [6]: http://%20https://cran.r-project.org/doc/manuals/R-exts.html
  [7]: https://www.rmetrics.org/files/Meielisalp2009/Presentations/Scott.pdf