Arc length as a statistical functional

The idea

Most statistical functionals summarise a curve by its height somewhere: a moment, a quantile, a density value. Arc length summarises it by how far you travel along it. For a quantile function \(Q\) on \((0,1)\) the arc length is

\[\int_0^1 \sqrt{1 + Q'(u)^2}\, du,\]

which is large when the curve is steep somewhere and small when it is flat. It is a measure of curve complexity rather than of location or spread, and it responds to structure that the usual functionals average away.

This package collects the arc-length constructions: a goodness-of-fit test, a family of distributions indexed by the shape of the quantile density, an equivalence family, the characteristic-function version, and a Bayesian test.

Arc length of a quantile family

The arcq family has quantile density \(Q'(u) = \sigma[1 + \sum_k c_k P_k(u)]_+\) with \(P_k\) the shifted Legendre polynomials.

o <- arcq(coef = c(0.5, -0.3), mu = 0, sigma = 1)
arclength(o)
#> [1] 1.43581

Arc length is a shape functional, so shifting the distribution cannot change it:

c(at_zero = arclength(arcq(0.5, mu = 0)), at_seventeen = arclength(arcq(0.5, mu = 17)))
#>      at_zero at_seventeen 
#>     1.429349     1.429349

The integrand is analytic on \([0,1]\), so Gauss-Legendre converges geometrically. Sixteen nodes already reach machine precision, where the equally spaced grid it replaced converges linearly:

vapply(c(8L, 16L, 32L, 64L), function(k) arclength(o, nodes = k), 0)
#> [1] 1.43581 1.43581 1.43581 1.43581

The integral is generally not elementary. With \(Q'\) of degree one the antiderivative is an inverse hyperbolic sine; with \(Q'\) of degree two it is an elliptic integral, which by Liouville’s theorem has no elementary antiderivative; beyond that it is hyperelliptic.

Goodness of fit

Under the null the probability integral transform is uniform, and the arc length of the resulting probability plot has a known distribution with an analytic saddlepoint tail. The test is powerful against local density structure — multimodality, clustering, heaping — and weak against smooth location and scale departures, which is the opposite of the empirical-distribution tests.

set.seed(1)
al_test(runif(200))$p.value                       # null holds
#> [1] 0.3974529
al_test(c(runif(100), rbeta(100, 8, 8)))$p.value  # a clustered middle
#> [1] 0.0002620575

Characteristic-function arc length

The same functional applied to the characteristic-function curve is scale free, and it has closed forms for several families. Two of them anchor the whole construction:

c(normal = cf_arclength_family("normal"),
  exponential = cf_arclength_family("exponential", lambda = 1),
  pi = pi)
#>      normal exponential          pi 
#>    2.000000    3.141593    3.141593

A symmetric monotone (Polya) law carries the value two, and the exponential carries \(\pi\) exactly. Because total arc length is scale free, the exponential rate cannot matter:

c(rate_1 = cf_arclength_family("exponential", lambda = 1),
  rate_2 = cf_arclength_family("exponential", lambda = 2))
#>   rate_1   rate_2 
#> 3.141593 3.141593

The empirical version approaches the theoretical one on a large sample:

set.seed(3)
cf_arclength(rnorm(50000))
#> [1] 2.013334

The equivalence family

A separate strand asks when two readings of a curve agree. The shoulder of a quantile density solves \(3q'^2 = q q''\), transcendental in general. On the two-exponent family \(q(u) = u^{\alpha}(1-u)^{\beta}\) it collapses to a quadratic, so the shoulder is available in closed form, and the returned root satisfies the original equation:

al <- -0.60; be <- -0.35
ub <- eq_ub_quad(al, be)[1]
gp  <- al / ub - be / (1 - ub)
gpp <- -al / ub^2 - be / (1 - ub)^2
c(shoulder = ub, residual_of_original_equation = abs(2 * gp^2 - gpp))
#>                      shoulder residual_of_original_equation 
#>                     0.1231096                     0.0000000

Equivalence itself is the vanishing of a closed-form discrepancy, and the solution curve is increasing in \(\alpha\):

als <- seq(-0.62, -0.54, by = 0.02)
bss <- vapply(als, eq_bstar, 0)
rbind(alpha = als, beta_star = round(bss, 6))
#>                [,1]      [,2]      [,3]      [,4]      [,5]
#> alpha     -0.620000 -0.600000 -0.580000 -0.560000 -0.540000
#> beta_star -0.578076 -0.378118 -0.243343 -0.148972 -0.081551
c(discrepancy_on_the_curve = eq_E(-0.60, eq_bstar(-0.60)))
#> discrepancy_on_the_curve 
#>             6.661338e-16

L-moments

Population and sample L-moments are included, since the families above are estimated by matching them. A uniform grid has first L-moment \(1/2\) and L-scale \(1/6\) exactly:

u <- (seq_len(20000) - 0.5) / 20000
sample_lmoments(u, 4L)[1:2]
#> [1] 0.500000 0.166675

Implementation

The numerical primitives run in a shared C back-end that is also bound from Python, and the two front ends are checked against each other value by value.