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 good 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.RR 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.
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). For bounded distribution it is important to check what happens on bounds and if $f(x) = 0$ for $x$ outside the support of $f$.
It is important to 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 andNaN's, it returnsNaN's and throws warningsNaNs producedfor invalid parameter values in thed/p/qfunctions and returnsNA's and throwsNAs producedwarnings inrfunctions, etc.
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.