Writing Statistical Software as Executable Mathematics

Most software tests have the same basic shape.

Give a function an input.

Check that the output equals the expected value.

For example:

1
expect_equal(mean(c(1, 2, 3)), 2)

There is nothing wrong with that style of test. Statistical software needs ordinary unit tests too.

But many of the strongest tests available to us are not examples of the form

They are mathematical relations between different parts of the software.

A test produces a p-value and a rejection decision. A confidence set is constructed by inverting those decisions. A projection aggregates pointwise evidence over nuisance values. A scalar transformation should reduce to a coordinate projection when the transformation is the identity. A fast implementation should agree with a simple reference definition. Reordering an arbitrary grid should not change the inferential object.

Those are not isolated outputs.

They are equations.

My preferred approach to statistical software is therefore:

The test suite then becomes more than a collection of remembered answers.

It becomes executable mathematics.

An Example Output Is Weaker Than an Identity

Suppose a test statistic is (T), a Monte Carlo p-value is (p), and the test rejects at level (\alpha).

A conventional test might use one fixture and require

1
2
3
expect_equal(result<span class="math-expression-inline math-expression--source" data-math-alt="p_value, 0.04) expect_true(result" data-math-source="p_value, 0.04)
expect_true(result" aria-label="p_value, 0.04) expect_true(result" tabindex="0">$p_value, 0.04)
expect_true(result$</span>reject)

That checks one output.

But if the procedure is defined so that

then we possess a stronger statement.

For every valid fixture,

That is the contract I want the test suite to protect.

Why is it stronger?

Because a refactor can change both the p-value and the rejection decision in a mutually inconsistent way while still accidentally preserving the expected result for a handful of stored examples.

A relational test checks whether the pieces still fit together mathematically.

Decision Duality Should Survive Every Layer

Suppose a confidence set is obtained by inverting pointwise tests.

At parameter value (\theta), define

If rejection is defined by

then set membership is not an independent design choice.

It follows immediately that

Notice the strict inequality.

If the test rejects when

then a point with exactly

is rejected and therefore excluded from the inverted set.

This is a tiny boundary detail that can easily be lost when two different functions are written months apart.

The code for the test may use

1
reject <- p_value <= alpha

while the code for inversion later uses

1
accepted <- p_value >= alpha

Every ordinary example away from the boundary can pass.

The mathematics is still wrong.

A relational regression test should state the whole chain:

1
2
3
expect_identical(result<span class="math-expression-inline math-expression--source" data-math-alt="reject, result" data-math-source="reject, result" aria-label="reject, result" tabindex="0">$reject, result$</span>p_value <= alpha)
expect_identical(result<span class="math-expression-inline math-expression--source" data-math-alt="accepted, !result" data-math-source="accepted, !result" aria-label="accepted, !result" tabindex="0">$accepted, !result$</span>reject)
expect_identical(result<span class="math-expression-inline math-expression--source" data-math-alt="accepted, result" data-math-source="accepted, result" aria-label="accepted, result" tabindex="0">$accepted, result$</span>p_value > alpha)

This is not duplicated testing.

It is testing a theorem connecting three representations of the same decision.

Finite Monte Carlo Makes Boundary Logic More Important

Monte Carlo p-values often use the plus-one form

where (K) is the number of simulated statistics at least as extreme as the observed statistic and (B) is the number of simulated replicates.

Then the attainable p-values lie on the grid

The rejection rule

may also be represented by a simulated critical value (c_\alpha).

If the software reports both quantities, they should satisfy

The exact convention around ties matters.

A generic empirical quantile does not automatically produce a critical value decision-equivalent to the plus-one p-value for every ((B,\alpha)).

This is a perfect candidate for executable mathematics because the identity defines what the two outputs mean together.

A test can generate tie-heavy simulated statistics deliberately and verify the relation across many valid replication counts.

Some Significance Levels May Be Unattainable

Finite simulation produces another useful contract.

The smallest possible plus-one p-value is

Therefore if

rejection is impossible no matter how extreme the observed statistic is.

A well-designed implementation should not hide that fact behind an arbitrary empirical quantile.

The mathematical implication is

That should become a test.

The interesting cases in statistical software are often exactly these edge cases where the mathematics gives a discrete boundary and generic numerical utilities give something merely plausible.

Confidence-Set Construction Is an Equation

Suppose we evaluate pointwise tests on a finite grid

The represented inverted set is

Then the returned accepted values should equal exactly

Again, this sounds obvious.

That is precisely why it should be executable.

The set-building layer can be tested independently of any particular expected set:

1
2
expected <- grid[p_values > alpha]
expect_identical(accepted_values, expected)

This is more robust than storing one expected vector from one fixture.

If the pointwise p-values legitimately change after a numerical improvement, the relational test remains correct.

It protects the inferential construction rather than yesterday's floating-point output.

Projection Gives Us Another Contract

Now suppose the parameter is multidimensional,

where (\psi) is the target and (\lambda) is nuisance.

On a represented finite grid, define a profile p-value

The projected set is then

There are two identities to test:

and

Those equations are a much better specification than a prose sentence saying that the function "projects a confidence set."

They tell us exactly what every returned value must mean.

Reduction Identities Are Extremely Powerful

A complicated API often contains special cases that should reduce exactly to simpler APIs.

Those reductions are some of the best tests we can write because they compare independent code paths.

Suppose a general coordinate-projection function operates on a multidimensional parameter grid.

If the model has only one parameter, there is no nuisance dimension.

Therefore projection should reduce to ordinary inversion:

The software should satisfy the same identity.

Likewise, suppose a scalar projection accepts a transformation

If we choose

the arbitrary scalar projection should reduce to coordinate projection onto parameter (j).

So we obtain another contract:

A regression test can compare:

• represented target values,
• profile p-values,
• decisions,
• accepted values,
• connected components,
• witness semantics.

No hard-coded numerical answer is necessary.

Two independently implemented APIs become mutual checks on each other.

Think in Commutative Diagrams

There is a useful mathematical way to think about this.

Suppose two paths through the software are supposed to represent the same operation.

Then we want the diagram to commute.

For example,

The two paths should lead to the same result.

Software engineering usually talks about unit tests and integration tests.

For mathematical software, I also want commutativity tests.

If two mathematically equivalent routes through the API disagree, at least one route is wrong.

Invariance Tests Encode What Must Not Matter

Some mathematical statements say that a transformation should leave the answer unchanged.

These are equally valuable.

Suppose the rows of a finite parameter grid are merely a representation of a set.

Then permuting those rows should not change the inferential object.

If (\pi) is a row permutation,

and

we should have

after canonical ordering of the result.

The same applies to projected p-values and accepted target values.

This tests something fundamentally different from an ordinary fixture.

It asks whether the implementation depends accidentally on representation order.

But Not Everything Must Be Invariant

A good invariant test also requires knowing what is allowed to change.

Suppose the projection object reports a witness row attaining

If the input grid is permuted, the numeric row index of that witness can change even though the inferential result does not.

So the correct contract is not

It is

This distinction matters.

Overly rigid tests can be just as misleading as weak tests. They freeze incidental representation choices instead of protecting the mathematics.

Deterministic Tie-Breaking Can Be Part of the API

Sometimes the mathematical solution is not unique.

Suppose two nuisance rows satisfy

The profile p-value is unambiguous.

The witness is not.

If the software exposes a single witness row, it needs a deterministic convention.

For example:

That rule is not a theorem about statistics.

It is a theorem about the software contract.

Once documented, it should be tested with an exact tie fixture rather than hoping two floating-point calculations happen to tie.

A good fixture can deliberately make the nuisance coordinate irrelevant so several rows produce mathematically identical tests.

Then the tie is structural, not accidental.

Construct Fixtures From the Mathematics

This is a general principle.

If we want to test equality, build a fixture where equality holds exactly by construction.

If we want to test sign symmetry, generate one dataset and its negative.

If we want to test translation invariance, add a common constant.

If we want to test permutation invariance, permute the same represented set.

If we want to test a reduction identity, choose the transformation that makes the general case algebraically identical to the simple case.

If we want to test an empty selected set, construct moments that are all comfortably slack.

The mathematics should design the test data.

Random fixtures are useful for stress testing, but exact structural fixtures tell us why a failure occurred.

Numerical Refinement Needs Its Own Mathematics

Even numerical routines have useful executable contracts.

Suppose a confidence-set boundary is refined by bisection inside a bracket

whose endpoints have opposite decisions.

If convergence is defined by

then the software should report convergence only when that condition is actually true.

Now consider adjacent representable floating-point numbers.

There may be no machine number strictly between them.

Then the mathematical midpoint exists in (\mathbb R), but the numerical midpoint satisfies

The algorithm cannot subdivide the bracket further.

That is not convergence to tolerance.

It is a numerical stall.

The contract should distinguish

This is another example where precise vocabulary becomes executable behaviour.

A Reference Implementation Can Be a Mathematical Specification

Optimized statistical code is often much harder to inspect than the formula it implements.

Suppose a direct method computes a quantity in (O(n^2)) memory while an optimized derivation computes it in (O(n)).

The direct version may be too slow for production but excellent as a reference.

Then test

across carefully chosen inputs.

I like retaining such simple reference implementations when they are small enough.

They are executable definitions.

The production code can be aggressively optimized while the test suite continues to ask whether the optimization represents the same mathematics.

Exact Equality and Numerical Equivalence Are Different Contracts

Not every identity should use the same comparison operator.

Some relations are discrete and should be exact:

• rejection decisions,
• accepted grid membership,
• names,
• dimensions,
• deterministic ordering,
• RNG-state restoration.

Others involve reassociated floating-point arithmetic and should be tested numerically:

The tolerance (\tau) should come from numerical reasoning, not convenience.

A test that uses all.equal() everywhere can hide meaningful errors.

A test that demands bit identity after a stable algebraic optimization can reject perfectly valid code.

The mathematical contract tells us which one is appropriate.

Documentation and Tests Should Share the Same Equation

One practice I find especially useful is to write the defining identity in the documentation and then mirror that identity in the regression test.

For example, documentation says

The test says

1
expect_identical(accepted, p_value > alpha)

Documentation says

The test recomputes that maximum independently.

Documentation says identity scalar projection reduces to coordinate projection.

The test calls both APIs and compares the inferential fields.

This creates a useful triangle:

If one corner changes, the mismatch becomes visible.

Do Not Test Only the Happy Interior

Mathematical contracts often become most informative at boundaries.

Test

Test a tied simulated critical value.

Test one represented grid point.

Test an empty set.

Test all points accepted.

Test no nuisance coordinates.

Test exact duplicate maximizers.

Test adjacent floating-point numbers.

Test a significance level unattainable at the chosen Monte Carlo replication count.

Those cases are not pathological distractions.

They are where conventions become observable.

A Statistical API Is a Network of Identities

As a package grows, it is tempting to think of it as a list of exported functions.

I find another mental model more useful.

Think of it as a graph.

The nodes are procedures:

• statistic construction,
• pointwise testing,
• p-value calculation,
• critical values,
• inversion,
• refinement,
• projection,
• transformations,
• summaries.

The edges are mathematical relationships:

Bugs often appear on the edges rather than inside the nodes.

Each individual function may look locally reasonable while two functions disagree about what the same mathematical quantity means.

Relational tests target those edges directly.

Why This Matters More in Statistical Software

Ordinary application software can often be specified by examples and business rules.

Statistical software makes claims about mathematical objects.

A function may claim to return a level-(1-\alpha) inverted set. Another may claim to project it. Another may report the p-value supporting the projected decision.

If those objects do not satisfy their defining relations, the problem is not merely an implementation bug.

The software is representing a different statistical procedure from the one its documentation describes.

That is why I want equations in the tests.

The Workflow I Prefer

When implementing a statistical method, I now like to proceed in roughly this order.

First, write the mathematical object.

For example,

Second, identify relations to objects already implemented.

For example,

Third, identify reductions and invariances.

For example,

in one dimension.

Fourth, decide which conventions are needed where mathematics is non-unique.

For example, first-max witness selection under ties.

Fifth, implement the algorithm.

Sixth, write regression tests that reconstruct the relations independently.

That order matters because it makes the test derive from the method rather than from whatever the first implementation happened to do.

Tests Should Not Merely Ratify the Code

A weak testing workflow is:

1. write the implementation,
2. run it once,
3. copy the output into the test,
4. declare regression coverage.

That can be useful for detecting changes, but it does not tell us whether the original output was correct.

A stronger workflow is:

1. derive a relation independently,
2. construct a fixture where the relation is informative,
3. compute each side through different code paths,
4. require them to agree.

The test then has an epistemic source outside the implementation itself.

The Best Regression Test Often Looks Like a Theorem

Some of my favorite statistical tests can be written almost directly as mathematics:

for every permutation (\pi) when row order is irrelevant.

These equations are concise because the mathematics has already done most of the specification work.

The Deeper Point

Statistical software is not just code that happens to perform mathematics.

At its best, the architecture can reflect the mathematics itself.

Definitions become functions.

Theorems become relations between functions.

Invariances become metamorphic tests.

Special cases become reduction tests.

Non-uniqueness becomes an explicit deterministic convention.

Numerical limits become separate states rather than silently changing inferential meaning.

Then a refactor does not merely have to reproduce a collection of old numbers.

It has to preserve the mathematical structure of the method.

That is the standard I want from statistical software:

$$ \boxed{ \text{the implementation should not merely run the mathematics; the test suite should enforce it.} }