Pearson Update: Advanced Dependence Measures

• Pearson Update is a modern reworking of Pearson’s ideas, refining normalization and extending correlation measures to nonlinear, multivariate, and online settings.
• It introduces rearrangement correlation (r♯), which expands the measurable range for monotone dependence and outperforms traditional measures in accuracy.
• The update also improves statistical inference using Edgeworth calibration and corrects Pearson-based distances to ensure true metric properties in applied contexts.

“Pearson Update” (Editor's term) denotes a contemporary reworking of Pearsonian ideas across several distinct but connected domains. In this usage, the term covers revisions to Pearson’s $r$ as a dependence measure, refinements of its sampling-theoretic calibration, exact online update and sensitivity formulas, corrections to Pearson-based dissimilarities that fail metric axioms, extensions from bivariate Euclidean settings to multi-way and Riemannian geometries, and parallel developments in the Pearson family of densities and Pearson diffusions (Ai, 2022, Vrbik, 2022, Harary, 2024, Solo, 2019, Taylor, 2020, Michl, 2020, Afendras et al., 2012, Beghin et al., 11 May 2025). Taken together, these results do not replace Pearson’s framework so much as sharpen its normalization, delimit its proper scope, and transplant it into settings where the original formulation is either conservative or formally inadequate.

1. Reassessing what Pearson’s $r$ can measure

Pearson’s $r$, the most widely-used correlation coefficient, is traditionally regarded as exclusively capturing linear dependence. The central revision advanced by Ai is that this exclusion of Pearson’s $r$ from nonlinear settings is too coarse: the real issue is not covariance itself, but the bound used to normalize it. Since Pearson’s $r$ is a scaled covariance, different scaling bounds yield coefficients with different capture ranges, and tighter bounds can expand these ranges rather than shrink them (Ai, 2022).

The classical normalization comes from the Cauchy–Schwarz inequality,

$$|\operatorname{Cov}(X,Y)| \le \sqrt{\operatorname{Var}(X)\operatorname{Var}(Y)}.$$

Ai replaces this with a rearrangement-based bound. For real random variables $X,Y$ with finite second moments, define the increasing rearrangements

$$X^\uparrow = F^{-1}(U), \qquad Y^\uparrow = G^{-1}(U),$$

and define

$$Y^\updownarrow = \begin{cases} Y^\uparrow, & E[XY]\ge 0,\ G^{-1}(1-U), & E[XY]<0. \end{cases}$$

Then

$$|E[XY]| \le |E[X^\uparrow Y^\updownarrow]| \le \sqrt{E[X^2]E[Y^2]}.$$

This refinement changes the interpretive status of Pearson-type normalization. Under the rearrangement bound, equality in the first inequality holds if and only if $r$0 are almost surely monotone, whereas equality in the Cauchy–Schwarz step retains the usual linear-scaling condition. A plausible implication is that the conventional “Pearson equals linear only” dictum confounds covariance with a particular denominator.

2. Rearrangement correlation and arbitrary monotone dependence

Ai’s new coefficient, called rearrangement correlation and written $r$1, is defined at the population level by

$r$2

with sample analogue

$r$3

where $r$4 are the sorted data vectors and $r$5 is $r$6 or $r$7 according to the sign of $r$8 (Ai, 2022).

The coefficient satisfies $r$9, and $r$0 if and only if the variables, or the sample vectors, are perfectly monotone increasing or decreasing. In linear settings it reduces exactly to Pearson’s $r$1: if $r$2, or in the sample case $r$3 up to permutation, then $r$4, so $r$5. It is sign-consistent, since $r$6, and it is monotone-invariant under strictly increasing or decreasing transformations. It also dominates classical Pearson in absolute value on monotone data: $r$7, with equality if and only if the relation is linear in distribution.

The connection to rank statistics is exact on ranked data. When applied to ranks $r$8, one has $r$9, which explains why Spearman’s $r$0 also attains the full range on monotone relationships, but only after discarding magnitude information. Rearrangement correlation preserves the covariance-based magnitude scale.

Its computation is simple. The sample algorithm first computes $r$1, determines its sign, sorts $r$2, sorts $r$3 increasingly or decreasingly according to that sign, computes $r$4, and returns $r$5. The dominating cost is sorting, so the complexity is $r$6, with all remaining steps $r$7.

The empirical comparison in Ai’s study is sharply favorable on monotone data. Across 50 canonical monotone functions with Gaussian noise, repeated 10 times at each prescribed $r$8-level and evaluated by mean absolute error against the true $r$9, $r$0 had the lowest MAE, approximately $r$1, ahead of Spearman’s $r$2 $r$3, distance correlation $r$4, Pearson $r$5 $r$6, Kendall’s $r$7 $r$8, MIC $r$9, Chatterjee’s $$|\operatorname{Cov}(X,Y)| \le \sqrt{\operatorname{Var}(X)\operatorname{Var}(Y)}.$$0 $$|\operatorname{Cov}(X,Y)| \le \sqrt{\operatorname{Var}(X)\operatorname{Var}(Y)}.$$1, additivity $$|\operatorname{Cov}(X,Y)| \le \sqrt{\operatorname{Var}(X)\operatorname{Var}(Y)}.$$2 $$|\operatorname{Cov}(X,Y)| \le \sqrt{\operatorname{Var}(X)\operatorname{Var}(Y)}.$$3, and HSIC $$|\operatorname{Cov}(X,Y)| \le \sqrt{\operatorname{Var}(X)\operatorname{Var}(Y)}.$$4. On five NIST monotone data sets—Chwirut1, Hahn1, Rat43, Roszman1, and Thurber—$$|\operatorname{Cov}(X,Y)| \le \sqrt{\operatorname{Var}(X)\operatorname{Var}(Y)}.$$5 again achieved the smallest MAE, $$|\operatorname{Cov}(X,Y)| \le \sqrt{\operatorname{Var}(X)\operatorname{Var}(Y)}.$$6, compared with $$|\operatorname{Cov}(X,Y)| \le \sqrt{\operatorname{Var}(X)\operatorname{Var}(Y)}.$$7 $$|\operatorname{Cov}(X,Y)| \le \sqrt{\operatorname{Var}(X)\operatorname{Var}(Y)}.$$8, MIC $$|\operatorname{Cov}(X,Y)| \le \sqrt{\operatorname{Var}(X)\operatorname{Var}(Y)}.$$9, distance correlation $X,Y$0, $X,Y$1 $X,Y$2, $X,Y$3 $X,Y$4, $X,Y$5 $X,Y$6, HSIC $X,Y$7, and $X,Y$8 $X,Y$9.

The principal limitation is equally clear. In 16 classic non-monotone scenarios, including functions with peaks, cycles, and circles, $$X^\uparrow = F^{-1}(U), \qquad Y^\uparrow = G^{-1}(U),$$0 performs poorly in absolute terms, with MAE approximately $$X^\uparrow = F^{-1}(U), \qquad Y^\uparrow = G^{-1}(U),$$1, even though it still outperforms $$X^\uparrow = F^{-1}(U), \qquad Y^\uparrow = G^{-1}(U),$$2. The coefficient is therefore a monotone-dependence measure, not a general-purpose dependence measure. The paper explicitly notes open problems on robustness to outliers, non-monotone dependence, multivariate or conditional extensions, and statistical inference for $$X^\uparrow = F^{-1}(U), \qquad Y^\uparrow = G^{-1}(U),$$3.

3. Sampling distribution, Fisher calibration, and online sensitivity

A separate update concerns the inferential treatment of the sample Pearson correlation. For i.i.d. samples from a bivariate Normal distribution with true correlation $$X^\uparrow = F^{-1}(U), \qquad Y^\uparrow = G^{-1}(U),$$4, the Edgeworth expansion can be used to approximate the distribution of a smooth transform $$X^\uparrow = F^{-1}(U), \qquad Y^\uparrow = G^{-1}(U),$$5. The leading $$X^\uparrow = F^{-1}(U), \qquad Y^\uparrow = G^{-1}(U),$$6 skewness term is

$$X^\uparrow = F^{-1}(U), \qquad Y^\uparrow = G^{-1}(U),$$7

Setting this to zero for all $$X^\uparrow = F^{-1}(U), \qquad Y^\uparrow = G^{-1}(U),$$8 forces

$$X^\uparrow = F^{-1}(U), \qquad Y^\uparrow = G^{-1}(U),$$9

whose nontrivial solution is $$Y^\updownarrow = \begin{cases} Y^\uparrow, & E[XY]\ge 0,\ G^{-1}(1-U), & E[XY]<0. \end{cases}$$0, recovering the classical Fisher transform. The same Edgeworth analysis then yields higher-order corrections: $$Y^\updownarrow = \begin{cases} Y^\uparrow, & E[XY]\ge 0,\ G^{-1}(1-U), & E[XY]<0. \end{cases}$$1

$$Y^\updownarrow = \begin{cases} Y^\uparrow, & E[XY]\ge 0,\ G^{-1}(1-U), & E[XY]<0. \end{cases}$$2

with $$Y^\updownarrow = \begin{cases} Y^\uparrow, & E[XY]\ge 0,\ G^{-1}(1-U), & E[XY]<0. \end{cases}$$3 and $$Y^\updownarrow = \begin{cases} Y^\uparrow, & E[XY]\ge 0,\ G^{-1}(1-U), & E[XY]<0. \end{cases}$$4. Replacing $$Y^\updownarrow = \begin{cases} Y^\uparrow, & E[XY]\ge 0,\ G^{-1}(1-U), & E[XY]<0. \end{cases}$$5 by $$Y^\updownarrow = \begin{cases} Y^\uparrow, & E[XY]\ge 0,\ G^{-1}(1-U), & E[XY]<0. \end{cases}$$6 in the correction terms gives the improved transform

$$Y^\updownarrow = \begin{cases} Y^\uparrow, & E[XY]\ge 0,\ G^{-1}(1-U), & E[XY]<0. \end{cases}$$7

whose Normal approximation error is $$Y^\updownarrow = \begin{cases} Y^\uparrow, & E[XY]\ge 0,\ G^{-1}(1-U), & E[XY]<0. \end{cases}$$8, compared with $$Y^\updownarrow = \begin{cases} Y^\uparrow, & E[XY]\ge 0,\ G^{-1}(1-U), & E[XY]<0. \end{cases}$$9 for the classical Fisher choice $$|E[XY]| \le |E[X^\uparrow Y^\updownarrow]| \le \sqrt{E[X^2]E[Y^2]}.$$0 (Vrbik, 2022).

The reported numerical error constants are materially smaller. At $$|E[XY]| \le |E[X^\uparrow Y^\updownarrow]| \le \sqrt{E[X^2]E[Y^2]}.$$1, the maximal CDF error for the classical Fisher approximation is approximately $$|E[XY]| \le |E[X^\uparrow Y^\updownarrow]| \le \sqrt{E[X^2]E[Y^2]}.$$2 at $$|E[XY]| \le |E[X^\uparrow Y^\updownarrow]| \le \sqrt{E[X^2]E[Y^2]}.$$3, $$|E[XY]| \le |E[X^\uparrow Y^\updownarrow]| \le \sqrt{E[X^2]E[Y^2]}.$$4 at $$|E[XY]| \le |E[X^\uparrow Y^\updownarrow]| \le \sqrt{E[X^2]E[Y^2]}.$$5, and $$|E[XY]| \le |E[X^\uparrow Y^\updownarrow]| \le \sqrt{E[X^2]E[Y^2]}.$$6 at $$|E[XY]| \le |E[X^\uparrow Y^\updownarrow]| \le \sqrt{E[X^2]E[Y^2]}.$$7; for the improved $$|E[XY]| \le |E[X^\uparrow Y^\updownarrow]| \le \sqrt{E[X^2]E[Y^2]}.$$8, the corresponding values are approximately $$|E[XY]| \le |E[X^\uparrow Y^\updownarrow]| \le \sqrt{E[X^2]E[Y^2]}.$$9, $r$00, and $r$01. In a Monte Carlo with $r$02 replicates from a bivariate Normal with $r$03 and $r$04, 95\% confidence intervals achieved empirical coverage $r$05 under the classical Fisher transform and $r$06 under the improved transform, with average interval length reduced by about $r$07.

Operationally, Pearson correlation has also been updated for online and streaming settings. Harary and collaborators develop Welford-based identities for maintaining running means $r$08, sums of squares $r$09, covariance sum $r$10, and hence

$r$11

When a new point $r$12 arrives, the updated correlation has a closed-form $r$13 expression: $r$14 Under a rectangular feasibility set $r$15, the maximal possible change

$r$16

is attained either at the four corners or at up to four edge-intersection points determined by the regression lines, leaving at most eight candidates. The same geometry governs worst-case p-value change $r$17, except for the additional zero-crossing case where $r$18 becomes attainable. If Welford summaries are already maintained, both sensitivities can be computed in $r$19; otherwise a single $r$20 pass suffices (Harary, 2024).

These two lines of work update Pearson in complementary senses. One recalibrates finite-sample inference by removing skewness at the $r$21 level; the other makes correlation and significance responsive to adversarial or merely newly arriving data.

4. Pearson-based dissimilarity and the metric correction

Another corrective update is negative rather than expansive: some Pearson-derived quantities that are commonly treated as distances are not metrics. Solo studies the Pearson distance

$r$22

and the sign-invariant variant

$r$23

showing that both can violate the triangle inequality (Solo, 2019).

The failures are explicit. For a valid correlation structure with $r$24, $r$25, and $r$26, one gets

$r$27

so $r$28. Likewise, with $r$29, $r$30, and $r$31,

$r$32

again violating the triangle inequality.

The repair is classical but often neglected in applications. On the space of zero-mean, unit-variance variables,

$r$33

is a metric because

$r$34

so $r$35. Likewise,

$r$36

is a metric on projective space, since it is the sine of the acute angle $r$37. The first repair preserves sign, whereas the second is sign-invariant.

This correction matters in precisely those settings in which Pearson distance has been used as if it were a metric, including clustering problems in gene expression analysis, brain imaging, and cyber security. A common misconception is that any bounded, symmetric function derived from correlation is a legitimate distance. Solo’s result shows that metric validity must be established rather than presumed.

5. Higher-dimensional and nonlinear geometric extensions

Pearson’s original coefficient is bivariate, Euclidean, and signed. Taylor’s multi-way correlation coefficient extends the idea of a single-number summary of dependence to $r$38 variables without designating any one variable as an outcome. Given observations $r$39, form the empirical correlation matrix

$r$40

let $r$41 be its eigenvalues, and define

$r$42

Since $r$43, the eigenvalue mean is $r$44. The coefficient lies in $r$45, equals $r$46 when $r$47, equals $r$48 under perfect linear dependence, is symmetric in its arguments, and is invariant under centering and rescaling. In the two-dimensional case it reduces exactly to $r$49, because the eigenvalues are $r$50 and $r$51. Its computational cost is $r$52: $r$53 to form $r$54 and $r$55 for the symmetric eigendecomposition (Taylor, 2020).

Taylor’s construction remains purely linear. Michl’s Riemann–Pearson correlation addresses the nonlinear case by replacing straight-line geometry with a smooth principal manifold $r$56. If $r$57 denotes orthogonal projection onto $r$58, define the reliability of coordinate $r$59 by

$r$60

and define the local sensitivity

$r$61

The squared Riemann–Pearson correlation is then

$r$62

If $r$63 is elliptically $r$64-distributed and $r$65 is in fact a linear subspace $r$66, then $r$67 and the construction reduces to the linear $r$68-correlation; in the bivariate linear case it recovers the square of the ordinary Pearson coefficient. The half-circle example demonstrates the intended gain: when data lie exactly on the upper half-circle $r$69, classical Pearson correlation vanishes by symmetry, yet the Riemann–Pearson correlation equals $r$70, reflecting perfect one-to-one nonlinear dependence (Michl, 2020).

These generalizations show two distinct routes beyond bivariate Pearson. One route aggregates linear interrelatedness by spectral dispersion of the correlation matrix; the other redefines correlation relative to a curved latent geometry. Neither should be confused with rearrangement correlation: $r$71 is joint and linear, Riemann–Pearson is nonlinear and manifold-based, and $r$72 is covariance-based but restricted to monotone dependence.

6. Pearson families beyond correlation: integrated densities and stretched diffusions

The name “Pearson” also denotes a family of probability laws defined by differential structure rather than by correlation. A density $r$73 on $r$74 belongs to the ordinary Pearson family if

$r$75

where $r$76 has degree at most $r$77, and $r$78 has degree at most $r$79. The Integrated Pearson family is the subclass of absolutely continuous laws with finite mean $r$80 for which there exists a quadratic

$r$81

such that

$r$82

This implies

$r$83

so every Integrated Pearson density solves a Pearson differential equation, but the converse need not hold; the integrated family is a strict subset. The review by Afendras and Papadatos identifies orthogonality of the first three Rodrigues-type polynomials with the Integrated Pearson property, gives a central-moment recurrence,

$r$84

establishes a covariance identity

$r$85

and organizes the system, up to affine transformations, into six types: Normal-type, Gamma-type, Beta-type, Student-type, Reciprocal-Gamma, and Snedecor-type. It also provides an endpoint-adjustment algorithm for deciding whether a given Pearson density belongs to the integrated subclass (Afendras et al., 2012).

Pearson diffusions extend the same tradition to stochastic processes. A classical Pearson diffusion solves

$r$86

with generator

$r$87

and Fokker–Planck operator

$r$88

Beghin, Leonenko, Papić, and Vaz define stretched non-local Pearson diffusions by time-changing $r$89 with a non-Markovian process $r$90 whose Laplace transform is the three-parameter Kilbas–Saigo function,

$r$91

The associated stretched Caputo derivative is

$r$92

and the Kilbas–Saigo function is its eigenfunction. The fractional Cauchy problem

$r$93

admits both analytic spectral solutions and stochastic representations through the time-changed process $r$94. The forward equation

$r$95

has an analogous expansion, and the limiting distribution remains the same invariant density $r$96 as in the untimed-changed Pearson diffusion. The framework also extends to stretched fractional hyperbolic, telegraph-type problems and yields a Mellin–Barnes representation and asymptotic formula for the Kilbas–Saigo function with complex argument (Beghin et al., 11 May 2025).

A persistent misconception is that “Pearson” names only the bivariate correlation coefficient. The broader literature shows instead a family of related constructions: covariance normalization, orthogonal polynomial systems, diffusion generators, and fractional time changes. The update is therefore both local and structural. Locally, it sharpens formulas such as Pearson’s $r$97 or Fisher’s transform. Structurally, it enlarges the Pearsonian program into domains where geometry, online computation, nonlocal dynamics, or orthogonality are the primary organizing principles.