---
title: "Lecture 4: Fréchet Means — Inference in General Metric Spaces"
subtitle: "Scalar variance asymptotics and Fréchet ANOVA"
number-sections: true
crossref:
chapters: true
format:
html:
code-fold: true
code-tools: true
code-copy: true
pdf:
documentclass: scrartcl
pdf-engine: xelatex
toc: true
number-sections: true
geometry:
- margin=1in
colorlinks: true
include-in-header:
text: |
% scrartcl has sections but no chapter counter. Quarto's chapter-style
% cross-references need this counter to produce section-qualified labels.
\newcounter{chapter}
\renewcommand{\thechapter}{\thesection}
\makeatletter
\AtBeginDocument{%
\@addtoreset{theorem}{section}%
\@addtoreset{proposition}{section}%
\@addtoreset{lemma}{section}%
\@addtoreset{codelisting}{section}%
}
\makeatother
bibliography: source/ref.bib
---
# Learning Goals
By the end of this lecture, learners should be able to:
- Explain why inference for Fréchet means is more delicate than point estimation and why scalar-based approaches are widely used.
- State the scalar CLT for the empirical Fréchet variance under metric-space regularity conditions.
- Formulate the $k$-sample Fréchet ANOVA hypotheses — equality of Fréchet means ($H_0^\mu$), equality of Fréchet variances ($H_0^V$), and their joint null ($H_0$) — and distinguish the mean-sensitive contrast, variance test, and omnibus test.
- State Propositions 3 and 4 of @DubeyMuller2019: the negligibility of $F_n$ under $H_0^\mu$ and the $\chi^2_{k-1}$ limit of the variance contrast $U_n$ under $H_0^V$.
- Define the population contrasts $F$ and $U$ and explain the power conditions $\sqrt n\,a_n\to\infty$ and $n b_n\to\infty$ from Theorem 3 of @DubeyMuller2019.
# Why Inference Is Different from Estimation
The preceding lectures treated existence, consistency, convergence rates, and computation of Fréchet means. **Inference** — hypothesis testing, confidence sets, $p$-values — is more delicate.
In Euclidean space, inference for a mean is straightforward:
1. Subtract the hypothesized mean from the observations.
2. Studentize the resulting vector (divide by an estimated covariance).
3. Appeal to a central limit theorem (CLT) in a **fixed linear space**.
In a general metric space, none of these steps may be available:
- **No subtraction**: There is no vector-space structure — you cannot form $X_i - \mu_0$.
- **No tangent space**: Without a smooth manifold structure, there is no canonical linear approximation.
- **No covariance matrix**: Even a notion of "direction" may be absent.
- **No smooth local chart**: Singularities, non-unique means, stickiness, and curvature can produce non-classical limits.
A broadly applicable strategy is to work through **scalar quantities** built from the Fréchet
function—especially squared distances and Fréchet variances. This lecture develops that
strategy from the scalar variance CLT to the $k$-sample Fréchet ANOVA and its power analysis.
# Preliminaries: Fréchet Variance and Its Scalar CLT
Following @DubeyMuller2019, let $(\Omega,d)$ be a **bounded** metric space and let
$X,X_1,\ldots,X_n$ be independent and identically distributed $\Omega$-valued random objects.
Boundedness is part of the paper's standing framework; in particular, it guarantees that all
powers of the distances used below have finite moments. Define the population Fréchet function,
mean, and variance by
$$
F(\omega) = \mathbb{E}\, d^2(X, \omega),
\qquad
\mu = \arg\min_{\omega \in \Omega} F(\omega),
\qquad
V = F(\mu),
\qquad
\sigma^2 = \operatorname{var}\{d^2(X, \mu)\}.
$$
Here $V$ is the **Fréchet variance**, whereas $\sigma^2$ is the ordinary variance of
the scalar random variable $d^2(X,\mu)$. Their empirical counterparts are
$$
\hat{F}_n(\omega) = \frac{1}{n}\sum_{i=1}^n d^2(X_i, \omega), \qquad
\hat{\mu}_n = \arg\min_{\omega \in \Omega} \hat{F}_n(\omega), \qquad
\hat{V}_n = \hat{F}_n(\hat{\mu}_n).
$$
## Regularity Assumptions
For consistency and asymptotic normality, @DubeyMuller2019 impose the following three
conditions. We label them (A1)–(A3), corresponding respectively to Assumptions 1–3 in the paper.
::: {.callout-note title="Assumptions (A1)–(A3)"}
**(A1) Existence, uniqueness, and separation.** The population and sample Fréchet means
$\mu$ and $\hat{\mu}_n$ exist and are unique, and for every $\varepsilon>0$,
$$
\inf_{d(\omega,\mu)>\varepsilon}F(\omega)>F(\mu).
$$
The separation condition ensures that approximate minimizers cannot remain a positive distance
from $\mu$ and yields $d(\hat{\mu}_n,\mu)=o_P(1)$.
**(A2) Local entropy.** For $\omega\in\Omega$, let
$B_\delta(\omega)=\{\omega'\in\Omega:d(\omega',\omega)<\delta\}$ and define
$$
J(\delta,\omega)
=
\int_0^1
\left[1+\log N\{\varepsilon\delta/2,B_\delta(\omega),d\}\right]^{1/2}
\,d\varepsilon.
$$
For every $\omega\in\Omega$, $\delta J(\delta,\omega)\to0$ as $\delta\downarrow0$.
**(A3) Global entropy.** The entropy integral of the whole space is finite:
$$
\int_0^1
\{1+\log N(\varepsilon,\Omega,d)\}^{1/2}
\,d\varepsilon<\infty.
$$
:::
Assumptions (A1)–(A3) are used for the scalar Fréchet-variance CLT. A useful distinction for
the ANOVA results below is that Propositions 3 and 4 and Theorem 2 of the paper require only
the groupwise versions of (A1)–(A2), together with existence and uniqueness of the pooled
sample Fréchet mean and stable sampling proportions. Assumption (A3) reappears in the paper's
power analysis under shrinking alternatives.
## Spaces Covered by These Entropy Assumptions
@DubeyMuller2019 note that these conditions are satisfied by several metric spaces. A useful sufficient entropy condition is
$$
\log N(\varepsilon, \Omega, d)
\le
K\varepsilon^{-\alpha}
\qquad\text{for some }K > 0\text{ and }0 < \alpha < 2.
$$
Indeed, this bound gives the global entropy integral because
$$
\int_0^1
\sqrt{1 + \log N(\varepsilon, \Omega, d)}\, d\varepsilon
<
\infty,
$$
and, applying the same bound to balls $B_\delta(\mu)$,
$$
J(\delta,\mu)
\lesssim
\int_0^1
\{1 + K(\varepsilon\delta)^{-\alpha}\}^{1/2}\, d\varepsilon
= O(\delta^{-\alpha/2}).
$$
Therefore
$$
\delta J(\delta,\mu) = O(\delta^{1-\alpha/2}) \to 0
\qquad\text{as }\delta \downarrow 0,
$$
which gives the local entropy condition. The remaining part of (A1) is an ordinary uniqueness and separation requirement for the relevant population and sample Fréchet means.
The examples highlighted by @DubeyMuller2019 include the following.
- **Univariate distributions under Wasserstein distance.** If the distributions have common compact support in $\mathbb{R}$, then the $L^2$-Wasserstein distance has the quantile representation
$$
d_W^2(F, G)
=
\int_0^1
\{F^{-1}(t) - G^{-1}(t)\}^2\, dt.
$$
Thus the object space can be represented by a bounded class of monotone quantile functions in $L^2(0,1)$. Entropy bounds for monotone functions yield the required local and global entropy conditions.
- **Correlation matrices of fixed dimension.** With Frobenius distance $d_F(A, B) = \|A - B\|_F$, the space of $p \times p$ correlation matrices is a bounded subset of a finite-dimensional Euclidean space. Hence its covering number grows at most polynomially in $1/\varepsilon$, so the metric entropy is only logarithmic in $1/\varepsilon$ and satisfies the entropy integrals.
- **Graph Laplacians of fixed-size weighted simple graphs.** For connected, undirected, simple graphs with a fixed number of vertices and edge weights bounded in $[0,1]$, the graph Laplacian is a symmetric matrix with entries in a bounded finite-dimensional set. Under Frobenius distance, the same finite-dimensional polynomial-covering argument applies. @DubeyMuller2019 describe this as a minor extension of the correlation-matrix argument.
- **Function classes with polynomial entropy exponent less than two.** Standard entropy bounds imply the sufficient condition above for classes such as distribution functions on $\mathbb{R}$ and $\mathbb{R}^2$ under $L^2$ distance, monotone functions from $\mathbb{R}$ into a compact subset of $\mathbb{R}$, Lipschitz functions on $[0,1]$ under $L^2$ distance, and suitable bounded classes of convex functions on compact convex subsets of $\mathbb{R}^d$; see Theorem 2.6.9 and related examples of @vanDerVaartWellner1996.
## The Fréchet Variance CLT
::: {#thm-frechet-variance-clt .theorem name="CLT for the empirical Fréchet variance"}
Under (A1)–(A3) (Theorem 1 of @DubeyMuller2019),
$$
\sqrt{n}(\hat{V}_n - V) \Rightarrow N(0, \sigma^2).
$$
The paper also establishes consistency (indeed, a CLT) for
$$
\hat{\sigma}_n^2 = \frac{1}{n}\sum_{i=1}^n \{d^2(X_i, \hat{\mu}_n) - \hat{V}_n\}^2,
$$
so that $\hat{\sigma}_n^2\to\sigma^2$ in probability and, when $\sigma>0$,
$$
\frac{\sqrt n(\hat V_n-V)}{\hat\sigma_n}\Rightarrow N(0,1).
$$
:::
The key reduction is Proposition 1 of @DubeyMuller2019:
$$
\frac1n\sum_{i=1}^n
\{d^2(\hat\mu_n,X_i)-d^2(\mu,X_i)\}
=o_P(n^{-1/2}).
$$
Consequently,
$$
\sqrt n(\hat V_n-V)
=
\frac1{\sqrt n}\sum_{i=1}^n\{d^2(\mu,X_i)-V\}+o_P(1),
$$
and the leading term obeys the ordinary scalar central limit theorem.
::: {.callout-tip title="Why a scalar CLT?"}
In a general metric space, there is no vector $\hat{\mu}_n - \mu$ to form a multivariate CLT around — there is no subtraction and no linear structure. The scalar quantity $d^2(X_i, \mu)$ is always well-defined, and averages of scalars obey the classical CLT under weak dependence conditions. This is the "lowest common denominator" of metric-space inference: reduce everything to real-valued random variables.
:::
# Fréchet ANOVA for $k$ Independent Samples
@DubeyMuller2019 use the preceding scalar asymptotics to compare $k$ populations of random
objects in terms of their Fréchet means and variances.
## Setup
Continue to assume that $(\Omega,d)$ is bounded. For groups $j=1,\ldots,k$, observe mutually
independent samples, with the observations within group $j$ independently and identically
distributed according to $P_j$:
$$
Y_{j1},\ldots,Y_{jn_j}\overset{\mathrm{iid}}{\sim}P_j,\qquad
n = \sum_{j=1}^k n_j, \qquad
\lambda_{j,n} = n_j / n.
$$
For group $j$, define the population quantities:
$$
F_j(\omega) = \mathbb{E}_j\, d^2(Y_j, \omega), \qquad
\mu_j = \arg\min_{\omega \in \Omega} F_j(\omega), \qquad
V_j = F_j(\mu_j), \qquad
\sigma_j^2 = \operatorname{var}_j\{d^2(Y_j, \mu_j)\}.
$$
It is useful to separate the two component null hypotheses. The **Fréchet mean null** is
$$
H_0^\mu: \mu_1 = \cdots = \mu_k,
$$
and the **Fréchet variance null** is
$$
H_0^V: V_1 = \cdots = V_k.
$$
The combined Fréchet ANOVA procedure targets the **joint null**
$$
H_0 = H_0^\mu \cap H_0^V.
$$
::: {.callout-note}
This is **not** a test of $P_1=\cdots=P_k$. Distinct distributions can have the same Fréchet
mean and variance. Conversely, rejection of the omnibus null shows that at least one mean or
variance differs, but the omnibus rejection alone does not identify which component differs.
:::
## Sample Quantities
For each group, compute:
$$
\hat{F}_j(\omega) = \frac{1}{n_j}\sum_{i=1}^{n_j} d^2(Y_{ji}, \omega), \qquad
\hat{\mu}_j \in \arg\min_{\omega \in \Omega} \hat{F}_j(\omega), \qquad
\hat{V}_j = \hat{F}_j(\hat{\mu}_j),
$$
and
$$
\hat{\sigma}_j^2 = \frac{1}{n_j}\sum_{i=1}^{n_j} \{d^2(Y_{ji}, \hat{\mu}_j) - \hat{V}_j\}^2.
$$
The pooled sample Fréchet mean and variance are
$$
\hat{\mu}_p
=
\arg\min_{\omega\in\Omega}
\frac1n\sum_{j=1}^k\sum_{i=1}^{n_j}d^2(Y_{ji},\omega),
\qquad
\hat{V}_p = \frac{1}{n}\sum_{j=1}^k \sum_{i=1}^{n_j} d^2(Y_{ji}, \hat{\mu}_p).
$$
## Mean-Sensitive Contrast
The **between-group Fréchet contrast** is
$$
F_n = \hat{V}_p - \sum_{j=1}^k \lambda_{j,n} \hat{V}_j.
$$
This is always nonnegative and measures the loss from using one pooled center instead of separate group centers — the metric analogue of between-group variation. In classical Euclidean ANOVA, $F_n$ corresponds to the weighted variance of the group means (the numerator of the $F$-ratio).
Proposition 3 of @DubeyMuller2019 establishes the following null behavior.
::: {#prp-frechet-anova-prop3 .proposition name="Negligibility of the mean contrast under equal Fréchet means"}
Suppose that $\hat{\mu}_p$ and $\hat{\mu}_j$ exist and are unique almost surely for all $j = 1, \ldots, k$. Let $0 < \lambda_{j,n} < 1$ with $\lambda_{j,n} \to \lambda_j \in (0, 1)$, where $\sum_{j=1}^k \lambda_j = 1$. Then, under the null hypothesis $H_0^\mu$ of equal Fréchet means and under assumptions (A1)--(A2) for each group,
$$
n^{1/2} F_n = o_P(1) \qquad \text{as } n \to \infty.
$$
:::
The normalized mean-sensitive component used in the omnibus statistic is
$$
T_n^\mu = \frac{n F_n^2}{\sum_{j=1}^k \lambda_{j,n}^2 \hat{\sigma}_j^2}.
$$
Under $H_0^\mu$, consistency of the $\hat{\sigma}_j^2$ and @prp-frechet-anova-prop3 give
$T_n^\mu=o_P(1)$. Thus $T_n^\mu$ has a **degenerate** asymptotic null distribution: the
paper uses it as the mean-sensitive component of the omnibus statistic, not as a separately
calibrated asymptotic level-$\alpha$ test.
::: {.callout-tip title="What the mean-contrast proposition tells us"}
If the Fréchet means are equal, the loss from using a pooled center rather than separate group
centers is negligible relative to the $n^{-1/2}$ scale. Under the additional existence,
uniqueness, and separation assumptions of Proposition 5 of @DubeyMuller2019, $F_n$ converges to a population
contrast that is positive exactly when the group Fréchet means are not all equal.
:::
## Test of Equality of Fréchet Variances
The **variance contrast** is
$$
U_n = \sum_{1 \le j < \ell \le k} \frac{\lambda_{j,n}\lambda_{\ell,n}}{\hat{\sigma}_j^2 \hat{\sigma}_\ell^2} (\hat{V}_j - \hat{V}_\ell)^2.
$$
This is a studentized pairwise contrast of the group Fréchet variances — in the Euclidean case, it is a modified version of Levene's test, using squared distances from the group Fréchet means rather than absolute deviations. The next proposition gives its asymptotic distribution under $H_0^V$.
::: {#prp-frechet-anova-prop4 .proposition name="Asymptotic distribution of the variance contrast"}
Under the assumptions of @prp-frechet-anova-prop3 and under the null hypothesis $H_0^V$ of equal population Fréchet variances,
$$
\frac{n U_n}{\sum_{j=1}^k \lambda_{j,n} / \hat{\sigma}_j^2} \Rightarrow \chi^2_{k-1}
\qquad \text{as } n \to \infty.
$$
:::
Consequently, a level-$\alpha$ test of $H_0^V$ rejects when
$$
\frac{n U_n}{\sum_{j=1}^k \lambda_{j,n} / \hat{\sigma}_j^2} > \chi^2_{k-1, 1-\alpha},
$$
where $\chi^2_{k-1, 1-\alpha}$ is the $(1-\alpha)$ quantile of the $\chi^2_{k-1}$ distribution.
::: {.callout-tip title="Connection to Levene's test"}
In classical analysis of variance, Levene's test assesses equality of group variances by applying ANOVA to the absolute deviations $|Y_{ji} - \bar{Y}_j|$. Here, the role of the absolute deviation is played by the squared distance $d^2(Y_{ji}, \hat{\mu}_j)$ from the group Fréchet mean. The statistic $U_n$ studentizes pairwise differences $\hat{V}_j - \hat{V}_\ell$ by the estimated variances $\hat{\sigma}_j^2$ and $\hat{\sigma}_\ell^2$ of these squared distances, yielding a $\chi^2_{k-1}$ limit.
:::
## Combined Omnibus Test
The paper combines the variance statistic with the mean-sensitive component:
$$
T_n = \frac{n U_n}{\sum_{j=1}^k \lambda_{j,n}/\hat{\sigma}_j^2} + \frac{n F_n^2}{\sum_{j=1}^k \lambda_{j,n}^2 \hat{\sigma}_j^2}.
$$
The two nonnegative terms receive equal weight. The paper notes that, if prior information
favors one type of departure, a modified convex weighting may instead be considered; such a
modification requires its own calibration.
::: {#thm-frechet-anova .theorem name="Fréchet ANOVA limit"}
Suppose the assumptions of @prp-frechet-anova-prop3 hold: the group and pooled sample Fréchet means
exist and are unique almost surely, $\lambda_{j,n}\to\lambda_j\in(0,1)$, and (A1)–(A2)
hold for each group. Under the joint null $H_0=H_0^\mu\cap H_0^V$,
$$
T_n \Rightarrow \chi^2_{k-1}.
$$
Consequently, the test that rejects $H_0$ when $T_n > \chi^2_{k-1, 1-\alpha}$ has asymptotic level $\alpha$.
:::
| Role | Null Hypothesis | Statistic | Asymptotic Null Behavior | Paper Result |
|---|---|---|---|---|
| **Mean-sensitive component** | $H_0^\mu: \mu_1 = \cdots = \mu_k$ | $\frac{nF_n^2}{\sum_j\lambda_{j,n}^2\hat{\sigma}_j^2}$ | $o_P(1)$ (degenerate; no separate asymptotic level-$\alpha$ test is supplied) | @prp-frechet-anova-prop3 |
| **Variance test** | $H_0^V: V_1 = \cdots = V_k$ | $\frac{n U_n}{\sum_j \lambda_{j,n} / \hat{\sigma}_j^2}$ | $\chi^2_{k-1}$ | @prp-frechet-anova-prop4 |
| **Omnibus** | $H_0 = H_0^\mu \cap H_0^V$ | $T_n$ (sum of both) | $\chi^2_{k-1}$ | @thm-frechet-anova |
Under the joint null, the mean-sensitive term is $o_P(1)$, while the variance term supplies
the $\chi^2_{k-1}$ limit. The behavior under alternatives is developed next.
## Power Analysis
The power analysis in @DubeyMuller2019 covers both fixed alternatives and alternatives that
approach the joint null as the total sample size increases. It is based on population counterparts
of the sample contrasts $F_n$ and $U_n$.
### Population Pooled Quantities and Contrasts
Let $\lambda_{j,n}\to\lambda_j\in(0,1)$, with $\sum_{j=1}^k\lambda_j=1$, and let
$Y_j\sim P_j$. Define the pooled population Fréchet mean and its Fréchet variance by
$$
\mu_p
=
\arg\min_{\omega\in\Omega}
\sum_{j=1}^k\lambda_j\,\mathbb E_j\{d^2(\omega,Y_j)\},
\qquad
V_p
=
\sum_{j=1}^k\lambda_j\,\mathbb E_j\{d^2(\mu_p,Y_j)\}.
$$
The population mean and variance contrasts are
$$
F
=
V_p-\sum_{j=1}^k\lambda_jV_j,
\qquad
U
=
\sum_{1\le j<\ell\le k}
\frac{\lambda_j\lambda_\ell}{\sigma_j^2\sigma_\ell^2}
(V_j-V_\ell)^2.
$$
The mean contrast can also be written as
$$
F
=
\sum_{j=1}^k\lambda_j
\left[
\mathbb E_j\{d^2(\mu_p,Y_j)\}
-
\mathbb E_j\{d^2(\mu_j,Y_j)\}
\right].
$$
Each summand is nonnegative because $\mu_j$ minimizes the $j$th population Fréchet function.
Under uniqueness, equality holds in every summand exactly when $\mu_p=\mu_j$ for every $j$.
Thus
$$
F=0
\quad\Longleftrightarrow\quad
\mu_1=\cdots=\mu_k,
$$
whereas, provided the $\sigma_j^2$ are positive,
$$
U=0
\quad\Longleftrightarrow\quad
V_1=\cdots=V_k.
$$
The two population contrasts therefore measure the two types of departure targeted by the
omnibus statistic.
### Consistency Under Alternatives
::: {#prp-anova-population-consistency .proposition name="Consistency of the pooled mean and mean contrast"}
Suppose that $\hat\mu_p,\hat\mu_j,\mu_p$, and $\mu_j$ exist and are unique for
$j=1,\ldots,k$, almost surely for the sample-based quantities. Suppose also that, for every
$\varepsilon>0$,
$$
\inf_{d(\omega,\mu_p)>\varepsilon}
\sum_{j=1}^k\lambda_j\mathbb E_j\{d^2(\omega,Y_j)\}
>
\sum_{j=1}^k\lambda_j\mathbb E_j\{d^2(\mu_p,Y_j)\},
$$
and, for every $j$,
$$
\inf_{d(\omega,\mu_j)>\varepsilon}
\mathbb E_j\{d^2(\omega,Y_j)\}
>
\mathbb E_j\{d^2(\mu_j,Y_j)\}.
$$
If $0<\lambda_{j,n}<1$ and $\lambda_{j,n}\to\lambda_j\in(0,1)$, then
$$
d(\hat\mu_p,\mu_p)=o_P(1),
\qquad
|F_n-F|=o_P(1).
$$
Moreover, $F\ge0$, and $F=0$ if and only if the population Fréchet means
$\mu_1,\ldots,\mu_k$ are all equal.
:::
This is Proposition 5 of @DubeyMuller2019. Its proof establishes
$\hat V_p\to V_p$ in probability; the earlier groupwise results similarly give
$\hat V_j\to V_j$ and $\hat\sigma_j^2\to\sigma_j^2$. Consequently, under a fixed
mean alternative, $F_n$ converges to a positive constant. Under a fixed variance alternative,
$U_n$ converges to a positive constant, subject to the corresponding nondegeneracy conditions.
The appropriate component of $T_n$ therefore diverges.
### Shrinking Alternatives and the Power Function
Let
$$
c_\alpha=\chi^2_{k-1,1-\alpha},
\qquad
R_{n,\alpha}=\{T_n>c_\alpha\}
$$
be the rejection region of the asymptotic level-$\alpha$ omnibus test. For nonnegative
sequences $a_n$ and $b_n$, define the component-specific alternative classes
$$
\mathcal H_n^F=\{(U,F):F\ge a_n\},
\qquad
\mathcal H_n^U=\{(U,F):U\ge b_n\},
$$
and their union
$$
\mathcal H_n=\mathcal H_n^F\cup\mathcal H_n^U
=
\{(U,F):F\ge a_n\ \text{or}\ U\ge b_n\}.
$$
When $a_n\to0$ and $b_n\to0$, these classes describe alternatives that approach the joint
null. For $\mathcal A\in\{\mathcal H_n^F,\mathcal H_n^U,\mathcal H_n\}$, write
$$
\beta_{\mathcal A}
=
\inf_{(U,F)\in\mathcal A}
P(R_{n,\alpha}).
$$
This is the worst-case power over the indicated class. The paper uses
$\beta_{\mathcal H_n}$ for the relevant sequence of alternatives; separating the two
component classes makes the two rate conclusions explicit.
::: {#thm-anova-power .theorem name="Consistency under shrinking alternatives"}
Under the assumptions of @prp-anova-population-consistency and the global entropy condition
(A3), let $a_n\to0$ and $b_n\to0$. For every $\alpha>0$:
1. If $\sqrt n\,a_n\to\infty$, then
$\beta_{\mathcal H_n^F}\to1$.
2. If $n b_n\to\infty$, then
$\beta_{\mathcal H_n^U}\to1$.
:::
This is Theorem 3 of @DubeyMuller2019. It supplies sufficient separation rates for uniform
consistency of the omnibus test. Fixed alternatives are included as an immediate special case:
if $F$ or $U$ is bounded away from zero, the corresponding condition holds automatically.
If both rate conditions hold, then the worst-case power over the union also satisfies
$\beta_{\mathcal H_n}\to1$.
# Application: Fréchet ANOVA for Market Risk Regimes
In ordinary ANOVA, we compare group means of scalar or vector-valued observations. In **Fréchet ANOVA**, the observations may be general random objects in a metric space. This section applies the Fréchet ANOVA framework of @DubeyMuller2019 to a concrete finance problem: comparing covariance structures across market regimes.
**The random objects** are covariance matrices estimated from independent blocks of financial
asset returns:
$$
\Sigma_t \in \mathcal{S}_{++}^p,
$$
where $\mathcal{S}_{++}^p$ is the space of $p\times p$ symmetric positive-definite
(SPD) matrices. Each object summarizes the marginal variances and cross-asset covariances in
one return block.
::: {.callout-note title="Main question"}
Do calm, normal, and stress market periods have the same average covariance structure, or do their Fréchet mean covariance matrices differ on the SPD manifold?
:::
## Data Objects
Each simulated sample covariance matrix is treated as **one observation** in
$\mathcal{S}_{++}^p$. The six coordinates represent:
| Ticker | Interpretation |
|---|---|
| SPY | US equity market |
| TLT | Long-term US Treasury |
| GLD | Gold |
| XLF | Financial sector |
| XLE | Energy sector |
| QQQ | Technology-heavy equity |
The regime labels in this simulation are fixed by the data-generating mechanism; they are not
estimated from the same covariance matrices being tested.
## Data-Generating Mechanism
Let $g\in\{\mathrm{calm},\mathrm{normal},\mathrm{stress}\}$ and let $n_g=m$ be the
number of covariance objects per regime. For object $i=1,\ldots,m$, generate a fresh,
non-overlapping block of $L$ returns,
$$
R_{gi1},\ldots,R_{giL}
\overset{\mathrm{iid}}{\sim}
N_p\{0,\Sigma_g(s)\},
$$
independently across both $g$ and $i$. The observed random object is the unbiased block sample
covariance
$$
Y_{gi}
=
\frac1{L-1}
\sum_{\ell=1}^L
(R_{gi\ell}-\bar R_{gi})(R_{gi\ell}-\bar R_{gi})^\top.
$$
Because the return blocks do not overlap, $Y_{g1},\ldots,Y_{gm}$ are i.i.d. within each
regime and the three samples are mutually independent, matching the sampling structure used
in the Fréchet ANOVA theory.
::: {.callout-warning title="Real financial returns are generally dependent"}
The independence above is a simplifying feature of the simulation, not a realistic assumption
for a return series such as SPY. Returns may have weak serial dependence, and their magnitudes
and squares commonly exhibit volatility clustering. Consequently, even covariance matrices
computed from **non-overlapping** blocks can remain dependent across time; non-overlap removes
the mechanical sharing of observations but does not create independence.
For real time-series applications, the i.i.d. Fréchet ANOVA calibration used here is therefore
not automatically valid. Dependence-aware methods may require mixing or other weak-dependence
conditions, long-run variance estimation, block bootstrap or subsampling procedures, or an
explicit time-series model such as a multivariate GARCH or stochastic-volatility model.
Regime labels should also be defined exogenously or handled with methods that account for
data-dependent regime estimation.
:::
The slider value $s\in[1,4]$ is converted to
$\tau=(s-1)/3\in[0,1]$. Starting from a common baseline covariance $\Sigma_0$, each
regime covariance follows the log-Euclidean path
$$
\Sigma_g(s)
=
\exp\!\left[
(1-\tau)\log\Sigma_0+\tau\log\Sigma_g^{\mathrm{target}}
\right].
$$
The normal-regime target is $\Sigma_0$. The calm target has lower marginal volatility, while
the stress target has higher volatility and a stronger common equity factor. Therefore:
- at $s=1$, $\Sigma_{\mathrm{calm}}(s)=\Sigma_{\mathrm{normal}}(s)
=\Sigma_{\mathrm{stress}}(s)=\Sigma_0$, so the joint null is exactly true;
- increasing $s$ moves smoothly toward the regime-specific targets; and
- every population covariance remains SPD—no entrywise correlation clipping is used.
The marginal standard deviations used to construct the baseline and targets are:
| Asset | Baseline | Calm target | Stress target |
|---|---:|---:|---:|
| SPY | 0.16 | 0.08 | 0.30 |
| TLT | 0.18 | 0.12 | 0.22 |
| GLD | 0.15 | 0.10 | 0.18 |
| XLF | 0.20 | 0.11 | 0.35 |
| XLE | 0.22 | 0.14 | 0.38 |
| QQQ | 0.19 | 0.13 | 0.32 |
The calm target retains the baseline correlation matrix. For the stress target, each asset has
a two-factor loading $\ell_a$ and the off-diagonal correlations are
$\rho_{ab}=\ell_a^\top\ell_b$. The diagonal idiosyncratic variance is
$1-\|\ell_a\|^2>0$, so this construction is positive definite by design. The loadings are
$$
\ell_{\mathrm{SPY}}=(0.90,0),\quad
\ell_{\mathrm{TLT}}=(-0.15,0.35),\quad
\ell_{\mathrm{GLD}}=(0.05,0.45),
$$
$$
\ell_{\mathrm{XLF}}=(0.88,0),\quad
\ell_{\mathrm{XLE}}=(0.80,0.15),\quad
\ell_{\mathrm{QQQ}}=(0.92,0).
$$
The condition $L>p$ in the controls makes each sample covariance positive definite with
probability one under the Gaussian model.
## Geometry on SPD Matrices
We use the **log-Euclidean metric**, which is simple, stable, and admits a closed-form Fréchet mean:
$$
d_{\mathrm{LE}}(A, B) = \|\log A - \log B\|_F,
$$
where $\|\cdot\|_F$ is the Frobenius norm and $\log$ is the matrix logarithm. Under this metric, the sample Fréchet mean of a collection of SPD matrices has the closed form
$$
\hat{\mu}_g = \exp\!\left(\frac{1}{n_g}\sum_{i=1}^{n_g} \log \Sigma_{gi}\right),
$$
and the sample Fréchet variance is
$$
\hat{V}_g = \frac{1}{n_g}\sum_{i=1}^{n_g} d_{\mathrm{LE}}^2(\Sigma_{gi}, \hat{\mu}_g).
$$
## Fréchet ANOVA Hypothesis
The Fréchet ANOVA tests whether the group Fréchet means and variances are equal:
$$
H_0^\mu: \mu_{\text{calm}} = \mu_{\text{normal}} = \mu_{\text{stress}}, \qquad
H_0^V: V_{\text{calm}} = V_{\text{normal}} = V_{\text{stress}}.
$$
A rejection of $H_0^\mu$ indicates that the typical **covariance matrices** differ; the
difference may come from marginal variances, correlations, or both. A rejection of $H_0^V$
indicates that the dispersion of covariance matrices around their regime centers differs.
The implementation computes $F_n$, $U_n$, the variance statistic from
@prp-frechet-anova-prop4, and the omnibus statistic $T_n$ from
@thm-frechet-anova. Since there are $k=3$ groups, the variance and omnibus statistics are
compared with $\chi^2_2$. The normalized mean-sensitive term is displayed as a component of
$T_n$, not as a separately calibrated test.
::: {.callout-note title="Scope of the asymptotic illustration"}
The complete log-Euclidean SPD space is unbounded, whereas @DubeyMuller2019 formulate their
theorems for a bounded object space. The simulation illustrates the sample statistics and
their $\chi^2_2$ calibration. A literal application of the paper's theorem can be obtained by
restricting covariance eigenvalues to a fixed compact interval, or by establishing suitable
moment and entropy extensions.
:::
## Interactive Exploration: Simulated Market Regime Data
The demo generates independent return blocks, converts each block into one covariance object,
and computes the Fréchet ANOVA statistics.
```{ojs}
//| label: fig-frechet-anova-spd
//| fig-cap: "Interactive: Fréchet ANOVA for market risk regimes on SPD matrices"
viewof simulation_parameters = {
const form = Inputs.form([
Inputs.range([20, 160], {
step: 10, value: 100,
label: "Independent covariance objects per regime"
}),
Inputs.range([10, 40], {
step: 2, value: 20,
label: "Returns per covariance object L"
}),
Inputs.range([1, 4], {
step: 0.1, value: 2.0,
label: "Regime separation factor s (1 = exact null)"
}),
Inputs.range([1, 100], {
step: 1, value: 42,
label: "Random seed"
})
]);
form.style.display = "grid";
form.style.gridTemplateColumns =
"repeat(auto-fit, minmax(min(100%, 360px), 1fr))";
form.style.gap = "0.8rem 1.5rem";
form.style.width = "100%";
form.style.maxWidth = "900px";
form.style.alignItems = "end";
return form;
}
n_objects = simulation_parameters[0]
block_len = simulation_parameters[1]
regime_sep = simulation_parameters[2]
seed = simulation_parameters[3]
// SPD matrix utilities
function matSqrt(A) {
// Symmetric matrix square root via eigendecomposition
const n = A.length;
// Simple power iteration for eigendecomposition of small matrices
// Use a numerical approach: compute eigenvalues and eigenvectors
// For small 6x6 matrices we use a simple Jacobi-like approach
// Build identity for initialization
function eye(n) { return Array.from({length: n}, (_, i) => Array.from({length: n}, (_, j) => i === j ? 1 : 0)); }
// QR decomposition for eigendecomposition iteration
function qrDecomp(M) {
const m = M.length;
const Q = eye(m);
let R = M.map(r => r.slice());
for (let j = 0; j < m; j++) {
for (let i = m - 1; i > j; i--) {
const a = R[i-1][j], b = R[i][j];
const r = Math.sqrt(a*a + b*b);
if (r < 1e-15) continue;
const c = a / r, s = -b / r;
// Apply rotation to R
for (let k = j; k < m; k++) {
const r1 = R[i-1][k], r2 = R[i][k];
R[i-1][k] = c * r1 - s * r2;
R[i][k] = s * r1 + c * r2;
}
// Apply rotation to Q
for (let k = 0; k < m; k++) {
const q1 = Q[k][i-1], q2 = Q[k][i];
Q[k][i-1] = c * q1 - s * q2;
Q[k][i] = s * q1 + c * q2;
}
}
}
return {Q, R};
}
// Compute eigenvalues/eigenvectors via QR iteration
let Ak = A.map(r => r.slice());
let V = eye(n);
for (let iter = 0; iter < 100; iter++) {
const {Q, R} = qrDecomp(Ak);
Ak = R.map((row, i) => row.map((_, j) => {
let sum = 0;
for (let k = 0; k < n; k++) sum += row[k] * Q[k][j];
return sum;
}));
const Vnew = V.map((row, i) => row.map((_, j) => {
let sum = 0;
for (let k = 0; k < n; k++) sum += row[k] * Q[k][j];
return sum;
}));
V = Vnew;
// Check convergence (off-diagonal elements small)
let offDiag = 0;
for (let i = 0; i < n; i++)
for (let j = 0; j < n; j++)
if (i !== j) offDiag += Ak[i][j] * Ak[i][j];
if (offDiag < 1e-12) break;
}
// Extract eigenvalues and eigenvectors
const evals = Array.from({length: n}, (_, i) => Ak[i][i]);
const sqrtEvals = evals.map(v => Math.sqrt(Math.max(1e-15, v)));
// Compute sqrt: V * diag(sqrt(evals)) * V^T
const sqrtM = Array.from({length: n}, () => Array(n).fill(0));
for (let i = 0; i < n; i++)
for (let j = 0; j < n; j++)
for (let k = 0; k < n; k++)
sqrtM[i][j] += V[i][k] * sqrtEvals[k] * V[j][k];
return sqrtM;
}
function matLog(A) {
const n = A.length;
// Eigendecomposition via same QR iteration approach
function eye(n) { return Array.from({length: n}, (_, i) => Array.from({length: n}, (_, j) => i === j ? 1 : 0)); }
function qrDecomp(M) {
const m = M.length;
const Q = eye(m);
let R = M.map(r => r.slice());
for (let j = 0; j < m; j++) {
for (let i = m - 1; i > j; i--) {
const a = R[i-1][j], b = R[i][j];
const r = Math.sqrt(a*a + b*b);
if (r < 1e-15) continue;
const c = a / r, s = -b / r;
for (let k = j; k < m; k++) {
const r1 = R[i-1][k], r2 = R[i][k];
R[i-1][k] = c * r1 - s * r2;
R[i][k] = s * r1 + c * r2;
}
for (let k = 0; k < m; k++) {
const q1 = Q[k][i-1], q2 = Q[k][i];
Q[k][i-1] = c * q1 - s * q2;
Q[k][i] = s * q1 + c * q2;
}
}
}
return {Q, R};
}
let Ak = sym(A);
let V = eye(n);
let converged = false;
for (let iter = 0; iter < 100; iter++) {
const {Q, R} = qrDecomp(Ak);
Ak = R.map((row, i) => row.map((_, j) => {
let sum = 0;
for (let k = 0; k < n; k++) sum += row[k] * Q[k][j];
return sum;
}));
const Vnew = V.map((row, i) => row.map((_, j) => {
let sum = 0;
for (let k = 0; k < n; k++) sum += row[k] * Q[k][j];
return sum;
}));
V = Vnew;
let offDiag = 0;
for (let i = 0; i < n; i++)
for (let j = 0; j < n; j++)
if (i !== j) offDiag += Ak[i][j] * Ak[i][j];
if (offDiag < 1e-12) {
converged = true;
break;
}
}
if (!converged) throw new Error("Symmetric eigendecomposition failed to converge in matLog");
const evals = Array.from({length: n}, (_, i) => Ak[i][i]);
const logEvals = evals.map(v => {
if (!(v > 1e-12)) throw new Error(`matLog requires an SPD matrix; eigenvalue ${v}`);
return Math.log(v);
});
const logM = Array.from({length: n}, () => Array(n).fill(0));
for (let i = 0; i < n; i++)
for (let j = 0; j < n; j++)
for (let k = 0; k < n; k++)
logM[i][j] += V[i][k] * logEvals[k] * V[j][k];
return logM;
}
function matExp(A) {
const n = A.length;
function eye(n) { return Array.from({length: n}, (_, i) => Array.from({length: n}, (_, j) => i === j ? 1 : 0)); }
function qrDecomp(M) {
const m = M.length;
const Q = eye(m);
let R = M.map(r => r.slice());
for (let j = 0; j < m; j++) {
for (let i = m - 1; i > j; i--) {
const a = R[i-1][j], b = R[i][j];
const r = Math.sqrt(a*a + b*b);
if (r < 1e-15) continue;
const c = a / r, s = -b / r;
for (let k = j; k < m; k++) {
const r1 = R[i-1][k], r2 = R[i][k];
R[i-1][k] = c * r1 - s * r2;
R[i][k] = s * r1 + c * r2;
}
for (let k = 0; k < m; k++) {
const q1 = Q[k][i-1], q2 = Q[k][i];
Q[k][i-1] = c * q1 - s * q2;
Q[k][i] = s * q1 + c * q2;
}
}
}
return {Q, R};
}
let Ak = sym(A);
let V = eye(n);
let converged = false;
for (let iter = 0; iter < 100; iter++) {
const {Q, R} = qrDecomp(Ak);
Ak = R.map((row, i) => row.map((_, j) => {
let sum = 0;
for (let k = 0; k < n; k++) sum += row[k] * Q[k][j];
return sum;
}));
const Vnew = V.map((row, i) => row.map((_, j) => {
let sum = 0;
for (let k = 0; k < n; k++) sum += row[k] * Q[k][j];
return sum;
}));
V = Vnew;
let offDiag = 0;
for (let i = 0; i < n; i++)
for (let j = 0; j < n; j++)
if (i !== j) offDiag += Ak[i][j] * Ak[i][j];
if (offDiag < 1e-12) {
converged = true;
break;
}
}
if (!converged) throw new Error("Symmetric eigendecomposition failed to converge in matExp");
const evals = Array.from({length: n}, (_, i) => Ak[i][i]);
const expEvals = evals.map(v => Math.exp(v));
const expM = Array.from({length: n}, () => Array(n).fill(0));
for (let i = 0; i < n; i++)
for (let j = 0; j < n; j++)
for (let k = 0; k < n; k++)
expM[i][j] += V[i][k] * expEvals[k] * V[j][k];
return expM;
}
function logEuclideanDist(A, B) {
const logA = matLog(A);
const logB = matLog(B);
return Math.sqrt(frobeniusSq(logA, logB));
}
function logEuclideanMean(mats) {
const n = mats.length;
const p = mats[0].length;
// Average of matrix logs
const avgLog = Array.from({length: p}, () => Array(p).fill(0));
for (let k = 0; k < n; k++) {
const L = matLog(mats[k]);
for (let i = 0; i < p; i++)
for (let j = 0; j < p; j++)
avgLog[i][j] += L[i][j] / n;
}
return matExp(avgLog);
}
function frechetVar(mats, mean) {
const meanLog = matLog(mean);
const logs = mats.map(matLog);
return frechetVarLog(logs, meanLog);
}
function frobeniusSq(A, B) {
let sum = 0;
for (let i = 0; i < A.length; i++)
for (let j = 0; j < A.length; j++)
sum += (A[i][j] - B[i][j]) ** 2;
return sum;
}
function meanLogMatrix(logMats) {
const n = logMats.length;
const p = logMats[0].length;
const avgLog = Array.from({length: p}, () => Array(p).fill(0));
for (let k = 0; k < n; k++)
for (let i = 0; i < p; i++)
for (let j = 0; j < p; j++)
avgLog[i][j] += logMats[k][i][j] / n;
return avgLog;
}
function frechetVarLog(logMats, meanLog) {
let sum = 0;
for (let k = 0; k < logMats.length; k++)
sum += frobeniusSq(logMats[k], meanLog);
return sum / logMats.length;
}
// Matrix symmetrization
function sym(A) {
const n = A.length;
const S = Array.from({length: n}, () => Array(n).fill(0));
for (let i = 0; i < n; i++)
for (let j = 0; j < n; j++)
S[i][j] = 0.5 * (A[i][j] + A[j][i]);
return S;
}
// Symmetric Jacobi eigendecomposition with an explicit convergence check.
function symmetricEigenJacobi(A) {
const n = A.length;
const D = sym(A);
const V = Array.from({length: n}, (_, i) =>
Array.from({length: n}, (_, j) => i === j ? 1 : 0)
);
const tolerance = 1e-12;
const maxIterations = 100 * n * n;
let maxOffDiagonal = Infinity;
for (let iter = 0; iter < maxIterations; iter++) {
let p = 0, q = 1;
maxOffDiagonal = 0;
for (let i = 0; i < n; i++) {
for (let j = i + 1; j < n; j++) {
if (Math.abs(D[i][j]) > maxOffDiagonal) {
maxOffDiagonal = Math.abs(D[i][j]);
p = i;
q = j;
}
}
}
if (maxOffDiagonal < tolerance) break;
const app = D[p][p], aqq = D[q][q], apq = D[p][q];
const theta = 0.5 * Math.atan2(2 * apq, aqq - app);
const c = Math.cos(theta), s = Math.sin(theta);
for (let k = 0; k < n; k++) {
if (k === p || k === q) continue;
const dkp = D[k][p], dkq = D[k][q];
D[k][p] = D[p][k] = c * dkp - s * dkq;
D[k][q] = D[q][k] = s * dkp + c * dkq;
}
D[p][p] = c*c*app - 2*s*c*apq + s*s*aqq;
D[q][q] = s*s*app + 2*s*c*apq + c*c*aqq;
D[p][q] = D[q][p] = 0;
for (let k = 0; k < n; k++) {
const vkp = V[k][p], vkq = V[k][q];
V[k][p] = c * vkp - s * vkq;
V[k][q] = s * vkp + c * vkq;
}
}
if (maxOffDiagonal >= tolerance) {
throw new Error(`Jacobi eigendecomposition failed; residual ${maxOffDiagonal}`);
}
return {values: Array.from({length: n}, (_, i) => D[i][i]), vectors: V};
}
function symmetricSpectralMap(A, transform, requirePositive = false) {
const {values, vectors} = symmetricEigenJacobi(A);
const mapped = values.map(value => {
if (requirePositive && !(value > 1e-12)) {
throw new Error(`SPD matrix required; eigenvalue ${value}`);
}
return transform(value);
});
const n = A.length;
const result = Array.from({length: n}, () => Array(n).fill(0));
for (let i = 0; i < n; i++)
for (let j = 0; j < n; j++)
for (let k = 0; k < n; k++)
result[i][j] += vectors[i][k] * mapped[k] * vectors[j][k];
return sym(result);
}
function spdLog(A) {
return symmetricSpectralMap(A, Math.log, true);
}
function symmetricExp(A) {
return symmetricSpectralMap(A, Math.exp, false);
}
// Data generation and Fréchet ANOVA simulation
function runSPDAnova(nObjects, blockLen, sepFactor, rngSeed) {
function mulberry32(a) {
return function() {
a |= 0; a = a + 0x6D2B79F5 | 0;
var t = Math.imul(a ^ a >>> 15, 1 | a);
t = t + Math.imul(t ^ t >>> 7, 61 | t) ^ t;
return ((t ^ t >>> 14) >>> 0) / 4294967296;
};
}
const rng = mulberry32(rngSeed);
// Box-Muller for normal samples
function randn() {
const u1 = rng(), u2 = rng();
return Math.sqrt(-2 * Math.log(Math.max(u1, 1e-15))) * Math.cos(2 * Math.PI * u2);
}
const NAMES = ["SPY", "TLT", "GLD", "XLF", "XLE", "QQQ"];
const p = 6;
// Common null covariance and regime-specific targets
const vols_baseline = [0.16, 0.18, 0.15, 0.20, 0.22, 0.19];
const vols_calm_target = [0.08, 0.12, 0.10, 0.11, 0.14, 0.13];
const vols_stress_target = [0.30, 0.22, 0.18, 0.35, 0.38, 0.32];
// Correlation matrix base
const base_corr = [
[1.0, -0.2, 0.0, 0.7, 0.5, 0.8],
[-0.2, 1.0, 0.3, -0.1, -0.1, -0.15],
[0.0, 0.3, 1.0, 0.0, 0.1, 0.05],
[0.7, -0.1, 0.0, 1.0, 0.6, 0.65],
[0.5, -0.1, 0.1, 0.6, 1.0, 0.55],
[0.8, -0.15, 0.05, 0.65, 0.55, 1.0]
];
// A factor-model correlation matrix is automatically positive definite.
function correlationFromLoadings(loadings) {
const C = Array.from({length: p}, () => Array(p).fill(0));
for (let i = 0; i < p; i++) {
const normSq = loadings[i].reduce((s, x) => s + x*x, 0);
if (!(normSq < 1)) throw new Error("Every factor loading must have squared norm below one");
for (let j = 0; j < p; j++) {
C[i][j] = i === j
? 1
: loadings[i].reduce((s, x, h) => s + x * loadings[j][h], 0);
}
}
return C;
}
const stress_corr_target = correlationFromLoadings([
[ 0.90, 0.00], // SPY
[-0.15, 0.35], // TLT
[ 0.05, 0.45], // GLD
[ 0.88, 0.00], // XLF
[ 0.80, 0.15], // XLE
[ 0.92, 0.00] // QQQ
]);
function makeCov(vols, corr) {
const D = Array.from({length: p}, () => Array(p).fill(0));
for (let i = 0; i < p; i++) D[i][i] = vols[i];
const result = Array.from({length: p}, () => Array(p).fill(0));
for (let i = 0; i < p; i++)
for (let j = 0; j < p; j++)
for (let k = 0; k < p; k++)
result[i][j] += D[i][k] * corr[k][j];
// D * corr * D
const result2 = Array.from({length: p}, () => Array(p).fill(0));
for (let i = 0; i < p; i++)
for (let j = 0; j < p; j++)
for (let k = 0; k < p; k++)
result2[i][j] += result[i][k] * D[k][j];
return sym(result2);
}
function logEuclideanPath(A, B, tau) {
const logA = spdLog(A);
const logB = spdLog(B);
const pathLog = logA.map((row, i) =>
row.map((x, j) => (1 - tau) * x + tau * logB[i][j])
);
return symmetricExp(pathLog);
}
const baseline_cov = makeCov(vols_baseline, base_corr);
const calm_target_cov = makeCov(vols_calm_target, base_corr);
const stress_target_cov = makeCov(vols_stress_target, stress_corr_target);
const tau = (sepFactor - 1) / 3;
const true_cov_calm = logEuclideanPath(baseline_cov, calm_target_cov, tau);
const true_cov_normal = baseline_cov.map(row => row.slice());
const true_cov_stress = logEuclideanPath(baseline_cov, stress_target_cov, tau);
// Cholesky decomposition
function cholesky(A) {
const n = A.length;
const L = Array.from({length: n}, () => Array(n).fill(0));
for (let i = 0; i < n; i++) {
for (let j = 0; j <= i; j++) {
let sum = 0;
for (let k = 0; k < j; k++) sum += L[i][k] * L[j][k];
if (i === j) {
const pivot = A[i][i] - sum;
if (!(pivot > 1e-12)) {
throw new Error(`Cholesky requires an SPD matrix; pivot ${pivot}`);
}
L[i][j] = Math.sqrt(pivot);
} else {
L[i][j] = (A[i][j] - sum) / L[j][j];
}
}
}
return L;
}
function generateReturns(cov, nObs) {
const L = cholesky(cov);
const returns = [];
for (let t = 0; t < nObs; t++) {
const z = Array.from({length: p}, () => randn());
const r = Array(p).fill(0);
for (let i = 0; i < p; i++)
for (let j = 0; j <= i; j++)
r[i] += L[i][j] * z[j];
returns.push(r);
}
return returns;
}
function sampleCovariance(returns) {
const L = returns.length;
const mean = Array(p).fill(0);
for (let i = 0; i < L; i++)
for (let j = 0; j < p; j++)
mean[j] += returns[i][j] / L;
const cov = Array.from({length: p}, () => Array(p).fill(0));
for (let i = 0; i < L; i++)
for (let a = 0; a < p; a++)
for (let b = 0; b < p; b++)
cov[a][b] += (returns[i][a] - mean[a]) *
(returns[i][b] - mean[b]) / (L - 1);
return sym(cov);
}
// Every covariance object uses a fresh, non-overlapping return block.
function generateCovarianceObjects(cov, count, L) {
const objects = [];
for (let i = 0; i < count; i++) {
objects.push(sampleCovariance(generateReturns(cov, L)));
}
return objects;
}
const covs_calm = generateCovarianceObjects(true_cov_calm, nObjects, blockLen);
const covs_normal = generateCovarianceObjects(true_cov_normal, nObjects, blockLen);
const covs_stress = generateCovarianceObjects(true_cov_stress, nObjects, blockLen);
const logs_calm = covs_calm.map(spdLog);
const logs_normal = covs_normal.map(spdLog);
const logs_stress = covs_stress.map(spdLog);
// Group Fréchet means
const meanLog_calm = meanLogMatrix(logs_calm);
const meanLog_normal = meanLogMatrix(logs_normal);
const meanLog_stress = meanLogMatrix(logs_stress);
const mean_calm = symmetricExp(meanLog_calm);
const mean_normal = symmetricExp(meanLog_normal);
const mean_stress = symmetricExp(meanLog_stress);
// Group Fréchet variances
const var_calm = frechetVarLog(logs_calm, meanLog_calm);
const var_normal = frechetVarLog(logs_normal, meanLog_normal);
const var_stress = frechetVarLog(logs_stress, meanLog_stress);
function squaredDistanceVariance(logs, meanLog, variance) {
return logs.reduce((sum, L) => {
const centered = frobeniusSq(L, meanLog) - variance;
return sum + centered * centered;
}, 0) / logs.length;
}
const sigma2_calm = squaredDistanceVariance(logs_calm, meanLog_calm, var_calm);
const sigma2_normal = squaredDistanceVariance(logs_normal, meanLog_normal, var_normal);
const sigma2_stress = squaredDistanceVariance(logs_stress, meanLog_stress, var_stress);
// Pooled Fréchet mean
const all_logs = [...logs_calm, ...logs_normal, ...logs_stress];
const meanLog_pool = meanLogMatrix(all_logs);
const var_pool = frechetVarLog(all_logs, meanLog_pool);
// Fréchet ANOVA statistics from Dubey and Müller (2019)
const n_total = all_logs.length;
const lambdas = [nObjects/n_total, nObjects/n_total, nObjects/n_total];
const variances = [var_calm, var_normal, var_stress];
const sigma2 = [sigma2_calm, sigma2_normal, sigma2_stress];
if (sigma2.some(x => !(x > 1e-14))) {
throw new Error("The squared-distance variance is numerically degenerate");
}
const Fn_obs = var_pool -
lambdas.reduce((sum, w, j) => sum + w * variances[j], 0);
let Un = 0;
for (let j = 0; j < 3; j++) {
for (let ell = j + 1; ell < 3; ell++) {
Un += lambdas[j] * lambdas[ell] /
(sigma2[j] * sigma2[ell]) *
(variances[j] - variances[ell]) ** 2;
}
}
const varianceDenom = lambdas.reduce((sum, w, j) => sum + w / sigma2[j], 0);
const meanDenom = lambdas.reduce((sum, w, j) => sum + w*w*sigma2[j], 0);
const variance_stat = n_total * Un / varianceDenom;
const mean_component = n_total * Fn_obs * Fn_obs / meanDenom;
const omnibus_stat = variance_stat + mean_component;
// With k=3 groups, the asymptotic reference law is chi-square with df=2.
const p_variance = Math.exp(-variance_stat / 2);
const p_omnibus = Math.exp(-omnibus_stat / 2);
// Pairwise distances between group means
const d_calm_normal = Math.sqrt(frobeniusSq(meanLog_calm, meanLog_normal));
const d_calm_stress = Math.sqrt(frobeniusSq(meanLog_calm, meanLog_stress));
const d_normal_stress = Math.sqrt(frobeniusSq(meanLog_normal, meanLog_stress));
// Convert means to 2D arrays for display
function matTo2D(m) { return m.map(r => r.slice()); }
return {
NAMES,
nObjects, blockLen, sepFactor, tau,
mean_calm: matTo2D(mean_calm),
mean_normal: matTo2D(mean_normal),
mean_stress: matTo2D(mean_stress),
var_calm, var_normal, var_stress,
sigma2_calm, sigma2_normal, sigma2_stress,
var_pool,
Fn_obs, Un, variance_stat, mean_component, omnibus_stat,
p_variance, p_omnibus,
d_calm_normal, d_calm_stress, d_normal_stress
};
}
spd_result = runSPDAnova(n_objects, block_len, regime_sep, seed);
// Display results
function formatPvalue(p) {
if (p < 0.001) return '< 0.001';
if (p < 0.01) return p.toFixed(4);
return p.toFixed(3);
}
function renderHeatmap(mat, names, title, maxVal) {
const n = mat.length;
const size = 300;
const margin = {top: 30, right: 10, bottom: 50, left: 60};
const cellW = (size - margin.left - margin.right) / (n + 1);
const cellH = (size - margin.top - margin.bottom) / (n + 1);
// Value range for color scale (blue-white-red)
const absMax = maxVal || Math.max(...mat.flat().map(Math.abs));
function color(val) {
if (Math.abs(val) < 1e-8) return '#f8f9fa';
const t = Math.max(-1, Math.min(1, val / absMax));
if (t > 0) return `rgb(${Math.round(255*t)}, ${Math.round(255*(1-t))}, ${Math.round(255*(1-t))})`;
else return `rgb(${Math.round(255*(1+t))}, ${Math.round(255*(1+t))}, ${Math.round(255*(-t))})`;
}
const rects = [];
for (let i = 0; i < n; i++) {
for (let j = 0; j < n; j++) {
const x = margin.left + (j + 0.5) * cellW;
const y = margin.top + (i + 0.5) * cellH;
rects.push(`<rect x="${x}" y="${y}" width="${cellW - 2}" height="${cellH - 2}"
fill="${color(mat[i][j])}" rx="2"/>`);
rects.push(`<text x="${x + (cellW - 2)/2}" y="${y + (cellH - 2)/2 + 3}"
text-anchor="middle" font-size="${Math.max(7, 30/n)}px" fill="#212529">${mat[i][j].toFixed(3)}</text>`);
}
}
// Row/col labels
for (let i = 0; i < n; i++) {
rects.push(`<text x="${margin.left - 5}" y="${margin.top + (i + 1) * cellH + 2}"
text-anchor="end" font-size="10px" fill="#495057">${names[i]}</text>`);
rects.push(`<text x="${margin.left + (i + 1) * cellW}" y="${margin.top + n * cellH + 15}"
text-anchor="end" font-size="10px" fill="#495057"
transform="rotate(-45, ${margin.left + (i + 1) * cellW}, ${margin.top + n * cellH + 15})">${names[i]}</text>`);
}
return `<svg viewBox="0 0 ${size} ${size}" role="img"
style="display:block;width:100%;max-width:${size}px;height:auto;"
xmlns="http://www.w3.org/2000/svg">
<text x="${size/2}" y="${margin.top - 8}" text-anchor="middle" font-size="12px" font-weight="bold">${title}</text>
${rects.join('')}
</svg>`;
}
function renderAnovaComponents(meanComponent, varianceStat) {
const width = 560, height = 155;
const margin = {top: 18, right: 25, bottom: 35, left: 135};
const plotW = width - margin.left - margin.right;
const critical = 5.991; // chi-square_2,0.95
const xmax = Math.max(critical * 1.2, meanComponent, varianceStat, meanComponent + varianceStat);
const x = value => margin.left + Math.min(value / xmax, 1) * plotW;
const rows = [
{label: "Mean component", value: meanComponent, color: "#f59f00", y: 35},
{label: "Variance statistic", value: varianceStat, color: "#339af0", y: 78},
{label: "Omnibus Tₙ", value: meanComponent + varianceStat, color: "#845ef7", y: 121}
];
return `<svg width="${width}" height="${height}" xmlns="http://www.w3.org/2000/svg">
${rows.map(row => `
<text x="${margin.left - 8}" y="${row.y + 13}" text-anchor="end" font-size="11px" fill="#495057">${row.label}</text>
<rect x="${margin.left}" y="${row.y}" width="${Math.max(1, x(row.value)-margin.left)}" height="18"
fill="${row.color}" opacity="0.78" rx="3"/>
<text x="${Math.min(x(row.value)+5, width-30)}" y="${row.y + 13}" font-size="10px" fill="#343a40">${row.value.toFixed(3)}</text>
`).join('')}
<line x1="${x(critical)}" y1="${margin.top}" x2="${x(critical)}" y2="${height-margin.bottom+15}"
stroke="#c92a2a" stroke-width="2" stroke-dasharray="4,3"/>
<text x="${x(critical)}" y="12" text-anchor="middle" font-size="10px" fill="#c92a2a">χ²₂ 5% critical value</text>
</svg>`;
}
// Assemble HTML display
r = spd_result;
html`
<div style="font-family: system-ui, sans-serif; width:100%; max-width: 1000px;">
<div style="display:flex;gap:16px;flex-wrap:wrap;justify-content:center;margin-bottom:24px;">
<div style="flex:1 1 280px;max-width:310px;min-width:0;">
${renderHeatmap(r.mean_calm, r.NAMES, "Calm: Fréchet mean cov", 0.15)}
</div>
<div style="flex:1 1 280px;max-width:310px;min-width:0;">
${renderHeatmap(r.mean_normal, r.NAMES, "Normal: Fréchet mean cov", 0.15)}
</div>
<div style="flex:1 1 280px;max-width:310px;min-width:0;">
${renderHeatmap(r.mean_stress, r.NAMES, "Stress: Fréchet mean cov", 0.15)}
</div>
</div>
<div style="margin-bottom: 24px;">
<h4 style="margin: 0 0 8px 0;">Fréchet ANOVA Components</h4>
${renderAnovaComponents(r.mean_component, r.variance_stat)}
</div>
<div style="display: grid; grid-template-columns: 1fr 1fr; gap: 20px;">
<div>
<h4 style="margin: 0 0 8px 0;">Results</h4>
<table style="width: 100%; border-collapse: collapse;">
<tr><td style="padding: 4px 8px; border-bottom: 1px solid #dee2e6;">Mean contrast F<sub>n</sub></td>
<td style="padding: 4px 8px; border-bottom: 1px solid #dee2e6; text-align: right;">${r.Fn_obs.toFixed(5)}</td></tr>
<tr><td style="padding: 4px 8px; border-bottom: 1px solid #dee2e6;">Normalized mean component</td>
<td style="padding: 4px 8px; border-bottom: 1px solid #dee2e6; text-align: right;">${r.mean_component.toFixed(4)}</td></tr>
<tr><td style="padding: 4px 8px; border-bottom: 1px solid #dee2e6;">Variance statistic</td>
<td style="padding: 4px 8px; border-bottom: 1px solid #dee2e6; text-align: right;">${r.variance_stat.toFixed(4)}</td></tr>
<tr><td style="padding: 4px 8px; border-bottom: 1px solid #dee2e6;"><b>Omnibus T<sub>n</sub></b></td>
<td style="padding: 4px 8px; border-bottom: 1px solid #dee2e6; text-align: right;"><b>${r.omnibus_stat.toFixed(4)}</b></td></tr>
<tr><td style="padding: 4px 8px; border-bottom: 1px solid #dee2e6;">Variance-test p-value</td>
<td style="padding: 4px 8px; border-bottom: 1px solid #dee2e6; text-align: right;">${formatPvalue(r.p_variance)}</td></tr>
<tr><td style="padding: 4px 8px; border-bottom: 1px solid #dee2e6;"><b>Omnibus p-value</b></td>
<td style="padding: 4px 8px; border-bottom: 1px solid #dee2e6; text-align: right; font-weight: bold;
color: ${r.p_omnibus < 0.05 ? '#c92a2a' : '#2b8a3e'};">${formatPvalue(r.p_omnibus)}</td></tr>
<tr><td style="padding: 4px 8px;"><b>Omnibus conclusion</b></td>
<td style="padding: 4px 8px; text-align: right; font-weight: bold;">
${r.p_omnibus < 0.05 ? 'Reject joint H₀' : 'Fail to reject joint H₀'}</td></tr>
</table>
<p style="font-size: 0.82em; color: #868e96;">
The mean component has no separate χ² calibration; it is included in the omnibus statistic.
</p>
</div>
<div>
<h4 style="margin: 0 0 8px 0;">Group Fréchet Variances</h4>
<table style="width: 100%; border-collapse: collapse;">
<tr><td style="padding: 4px 8px; border-bottom: 1px solid #dee2e6;">Calm</td>
<td style="padding: 4px 8px; border-bottom: 1px solid #dee2e6; text-align: right;">${r.var_calm.toFixed(5)}</td></tr>
<tr><td style="padding: 4px 8px; border-bottom: 1px solid #dee2e6;">Normal</td>
<td style="padding: 4px 8px; border-bottom: 1px solid #dee2e6; text-align: right;">${r.var_normal.toFixed(5)}</td></tr>
<tr><td style="padding: 4px 8px; border-bottom: 1px solid #dee2e6;">Stress</td>
<td style="padding: 4px 8px; border-bottom: 1px solid #dee2e6; text-align: right;">${r.var_stress.toFixed(5)}</td></tr>
<tr><td style="padding: 4px 8px; border-bottom: 1px solid #dee2e6;"><b>Pooled</b></td>
<td style="padding: 4px 8px; border-bottom: 1px solid #dee2e6; text-align: right;"><b>${r.var_pool.toFixed(5)}</b></td></tr>
</table>
<h4 style="margin: 16px 0 8px 0;">Pairwise Distances Between Fréchet Means</h4>
<table style="width: 100%; border-collapse: collapse;">
<tr><td style="padding: 4px 8px; border-bottom: 1px solid #dee2e6;">Calm ↔ Normal</td>
<td style="padding: 4px 8px; border-bottom: 1px solid #dee2e6; text-align: right;">${r.d_calm_normal.toFixed(4)}</td></tr>
<tr><td style="padding: 4px 8px; border-bottom: 1px solid #dee2e6;">Calm ↔ Stress</td>
<td style="padding: 4px 8px; border-bottom: 1px solid #dee2e6; text-align: right;">${r.d_calm_stress.toFixed(4)}</td></tr>
<tr><td style="padding: 4px 8px;">Normal ↔ Stress</td>
<td style="padding: 4px 8px; text-align: right;">${r.d_normal_stress.toFixed(4)}</td></tr>
</table>
</div>
</div>
<p style="margin-top: 16px; font-size: 0.85em; color: #868e96;">
<b>Parameters:</b> ${r.nObjects} independent covariance objects per regime ×
${r.blockLen} returns per object, separation factor ${r.sepFactor.toFixed(1)}
(τ = ${r.tau.toFixed(2)}).
</p>
</div>
`
```
::: {.callout-tip title="Try these experiments"}
- **Set $s=1$.** All three population covariance matrices are identical, so the joint null is exactly true. Across repeated seeds, the $p$-values—not necessarily any single $p$-value—should follow the null calibration approximately.
- **Increase $s$.** The calm and stress covariances move smoothly away from the common baseline, so the mean-sensitive component and omnibus power generally increase.
- **Increase the number of covariance objects.** This increases the number of independent random objects in each group.
- **Increase $L$.** Each covariance object is estimated more precisely from its own return block. Blocks remain non-overlapping for every setting.
:::
# Key Takeaways
- **Inference is harder than estimation** in metric spaces because there is no subtraction, no tangent space, and no covariance matrix.
- The **scalar Fréchet-variance CLT** (@thm-frechet-variance-clt) reduces the leading term to an ordinary average of $d^2(X_i,\mu)$; assumptions (A1)–(A3) control existence, consistency, and the substitution of $\hat\mu_n$ for $\mu$.
- The quantities $V=\mathbb E\,d^2(X,\mu)$ and $\sigma^2=\operatorname{var}\{d^2(X,\mu)\}$ play different roles: $V$ is the Fréchet variance, whereas $\sigma^2$ is the CLT variance.
- In $k$ samples, $F_n$ is the **mean-sensitive pooled-loss contrast** and is negligible under $H_0^\mu$ as stated in @prp-frechet-anova-prop3. Its normalized square is an omnibus component, not a separately calibrated test.
- The studentized $U_n$ component tests $H_0^V$ and converges to $\chi^2_{k-1}$ under @prp-frechet-anova-prop4. Adding the normalized mean component gives the omnibus $T_n$ for the joint null, with the limit in @thm-frechet-anova.
- The population contrasts $F$ and $U$ identify mean and variance departures. Theorem 3 of @DubeyMuller2019 gives uniform consistency when $\sqrt n\,a_n\to\infty$ for mean separation or $n b_n\to\infty$ for variance separation.
- In the SPD application, independent non-overlapping return blocks are essential for the sampling assumptions. The log-Euclidean path makes $s=1$ an exact null and preserves positive definiteness for every separation value.
# Exercises
1. **Scalar CLT derivation sketch**: Explain why $\sqrt{n}(\hat{V}_n - V)$ is asymptotically normal even though $\hat{\mu}_n$ is an $M$-estimator in a metric space. What roles do the local entropy condition (A2) and global entropy condition (A3) play in the proof? <a href="javascript:void(0)" onclick="showSolution('l4-sol-1')" class="solution-link">📝 Show Solution</a>
2. **Two-sample Fréchet ANOVA**: For $k = 2$ groups, write down the mean-sensitive component, the variance statistic, and the omnibus statistic $T_n$. Simplify the variance statistic and show that it takes the form of a squared two-sample $z$-statistic comparing $\hat{V}_1$ and $\hat{V}_2$. What does @prp-frechet-anova-prop3 imply about the mean-sensitive component under $H_0^\mu$? <a href="javascript:void(0)" onclick="showSolution('l4-sol-2')" class="solution-link">📝 Show Solution</a>
3. **Population contrasts and fixed alternatives**: Define the population contrasts $F$ and $U$. Explain why $F=0$ identifies equality of Fréchet means and $U=0$ identifies equality of Fréchet variances. What happens to the relevant component of $T_n$ under a fixed alternative? <a href="javascript:void(0)" onclick="showSolution('l4-sol-3')" class="solution-link">📝 Show Solution</a>
4. **Shrinking alternatives**: Let $F\ge a_n$ or $U\ge b_n$, with $a_n,b_n\to0$. State the two sufficient rates from @thm-anova-power and explain why the mean and variance contrasts have different detection rates. <a href="javascript:void(0)" onclick="showSolution('l4-sol-4')" class="solution-link">📝 Show Solution</a>
<style>
.solution-link {
font-size: 0.9em;
text-decoration: none;
white-space: nowrap;
margin-left: 0.3em;
}
.solution-link:hover {
text-decoration: underline;
}
.solution-dialog {
padding: 0;
max-width: 720px;
}
.solution-dialog-header {
display: flex;
justify-content: space-between;
align-items: flex-start;
border-bottom: 1px solid #dee2e6;
padding: 1.25rem 1.5rem 1rem;
background: #f8f9fa;
border-radius: 8px 8px 0 0;
}
.solution-dialog-header h4 {
margin: 0;
font-size: 1.15rem;
}
.solution-dialog-close {
background: none;
border: 1px solid #adb5bd;
border-radius: 4px;
padding: 0.2rem 0.75rem;
cursor: pointer;
font-size: 0.9rem;
color: #495057;
white-space: nowrap;
flex-shrink: 0;
}
.solution-dialog-close:hover {
background: #e9ecef;
}
.solution-original {
padding: 1rem 1.5rem;
background: #f1f3f5;
border-left: 4px solid #868e96;
margin: 1rem 1.5rem;
border-radius: 4px;
font-size: 0.95rem;
}
.solution-original strong {
color: #495057;
}
.solution-answer {
padding: 0.5rem 1.5rem 1.5rem;
}
.solution-answer strong {
color: #2b8a3e;
}
dialog {
border: none;
border-radius: 8px;
box-shadow: 0 8px 32px rgba(0,0,0,0.22);
padding: 0;
max-width: 750px;
width: 90vw;
}
dialog::backdrop {
background: rgba(0,0,0,0.45);
}
</style>
<dialog id="l4-sol-1">
<div class="solution-dialog">
<div class="solution-dialog-header">
<h4>Exercise 1: Scalar CLT Derivation Sketch</h4>
<button onclick="closeSolution('l4-sol-1')" class="solution-dialog-close">✕ Close</button>
</div>
<div class="solution-original">
<strong>Exercise:</strong> Explain why $\sqrt{n}(\hat{V}_n - V)$ is asymptotically normal even though $\hat{\mu}_n$ is an M-estimator in a metric space. What roles do the local entropy condition (A2) and global entropy condition (A3) play?
</div>
<div class="solution-answer">
<strong>Solution:</strong>
Write the decomposition:
$$
\sqrt{n}(\hat{V}_n - V) = \underbrace{\sqrt{n}\left(\frac{1}{n}\sum_{i=1}^n d^2(X_i, \mu) - V\right)}_{\text{Term 1}}
+ \underbrace{\sqrt{n}\left(\frac{1}{n}\sum_{i=1}^n [d^2(X_i, \hat{\mu}_n) - d^2(X_i, \mu)]\right)}_{\text{Term 2}}.
$$
**Term 1** is an ordinary sample average of the i.i.d. scalar random variables $d^2(X_i, \mu)$. By the classical CLT, it converges to $N(0, \sigma^2)$ where $\sigma^2 = \operatorname{var}(d^2(X, \mu))$.
**Term 2** is the "correction" for using $\hat{\mu}_n$ instead of $\mu$. The key is to show that Term 2 is $o_P(1)$ — asymptotically negligible. This requires controlling the empirical process
$$
\mathbb{G}_n(\omega) = \sqrt{n}\left(\frac{1}{n}\sum_{i=1}^n d^2(X_i, \omega) - \mathbb{E}\, d^2(X, \omega)\right)
$$
uniformly over $\omega$ in a neighborhood of $\mu$.
**Roles of the entropy conditions:** The local condition (A2),
$\delta J(\delta,\mu)\to0$, controls the oscillation of the squared-distance empirical process
over shrinking neighborhoods of $\mu$. The global condition (A3) supplies the whole-space
entropy control assumed in Proposition 1 of @DubeyMuller2019. Together with (A1) and
$d(\hat{\mu}_n,\mu)=o_P(1)$, these conditions yield
$$
\frac1n\sum_{i=1}^n
\{d^2(X_i,\hat\mu_n)-d^2(X_i,\mu)\}
=o_P(n^{-1/2}),
$$
which is precisely the statement needed to make Term 2 $o_P(1)$. The asymptotic
distribution is therefore driven entirely by Term 1.
**Key insight:** The scalar CLT works because $d^2(X_i, \mu)$ is a real-valued random variable — the metric-space complexity only enters through the need to control the substitution error.
</div>
</div>
</dialog>
<dialog id="l4-sol-2">
<div class="solution-dialog">
<div class="solution-dialog-header">
<h4>Exercise 2: Two-Sample Fréchet ANOVA</h4>
<button onclick="closeSolution('l4-sol-2')" class="solution-dialog-close">✕ Close</button>
</div>
<div class="solution-original">
<strong>Exercise:</strong> For $k = 2$ groups, write down the mean-sensitive component, the variance statistic, and the omnibus statistic $T_n$. Simplify the variance statistic and show that it takes the form of a squared two-sample $z$-statistic comparing $\hat{V}_1$ and $\hat{V}_2$. What does @prp-frechet-anova-prop3 imply about the mean-sensitive component under $H_0^\mu$?
</div>
<div class="solution-answer">
<strong>Solution:</strong>
**Components and statistics for $k=2$:** Let $\lambda_{1,n}=n_1/n$ and $\lambda_{2,n}=n_2/n$.
- **Mean-sensitive component:** $T_n^\mu = \dfrac{n F_n^2}{\lambda_{1,n}^2 \hat{\sigma}_1^2 + \lambda_{2,n}^2 \hat{\sigma}_2^2}$, where $F_n = \hat{V}_p - \lambda_{1,n}\hat{V}_1 - \lambda_{2,n}\hat{V}_2$. This is a component of the omnibus statistic, not a separately calibrated asymptotic test in the paper.
- **Variance statistic:** the variance contrast simplifies to a single pairwise term:
$$
U_n = \frac{\lambda_{1,n}\lambda_{2,n}}{\hat{\sigma}_1^2 \hat{\sigma}_2^2} (\hat{V}_1 - \hat{V}_2)^2.
$$
The variance-focused statistic is therefore
$$
\frac{n U_n}{\frac{\lambda_{1,n}}{\hat{\sigma}_1^2} + \frac{\lambda_{2,n}}{\hat{\sigma}_2^2}}
= \frac{n \lambda_{1,n}\lambda_{2,n} (\hat{V}_1 - \hat{V}_2)^2}{\lambda_{1,n}\hat{\sigma}_2^2 + \lambda_{2,n}\hat{\sigma}_1^2}
= \frac{n_1 n_2 (\hat{V}_1 - \hat{V}_2)^2}{n_1 \hat{\sigma}_2^2 + n_2 \hat{\sigma}_1^2}.
$$
When $\hat{\sigma}_1^2 \approx \hat{\sigma}_2^2 \approx \hat{\sigma}^2$, this becomes $\frac{n_1 n_2}{n_1 + n_2} \frac{(\hat{V}_1 - \hat{V}_2)^2}{\hat{\sigma}^2}$, the square of a two-sample $z$-statistic applied to the per-group Fréchet variances.
- **Omnibus:** $T_n$ is the sum of the variance statistic and the mean-sensitive component.
**@prp-frechet-anova-prop3** states that under $H_0^\mu: \mu_1 = \mu_2$, $n^{1/2}F_n = o_P(1)$. Consequently, $T_n^\mu = o_P(1)$ — the mean-focused term is asymptotically negligible when the Fréchet means are truly equal. Under the joint null $H_0 = H_0^\mu \cap H_0^V$, the omnibus statistic $T_n$ is therefore driven by the variance-focused term and converges to $\chi^2_1$.
**Power interpretation:** Under the additional assumptions of Proposition 5 of @DubeyMuller2019, if $H_0^\mu$ is false but $H_0^V$ holds, then the population mean contrast is positive, $T_n^\mu$ diverges, and the variance term stays bounded. Thus the omnibus test can detect mean differences even when the Fréchet variances are equal.
**Connection to classical ANOVA:** In Euclidean one-way ANOVA, the $F$-statistic compares between-group to within-group variability. Here, $F_n$ measures between-group variability through the **Fréchet variance loss** from pooling, while $U_n$ measures heterogeneity of the **Fréchet variances** themselves. Their combination gives power against either type of departure.
</div>
</div>
</dialog>
<dialog id="l4-sol-3">
<div class="solution-dialog">
<div class="solution-dialog-header">
<h4>Exercise 3: Population Contrasts and Fixed Alternatives</h4>
<button onclick="closeSolution('l4-sol-3')" class="solution-dialog-close">✕ Close</button>
</div>
<div class="solution-original">
<strong>Exercise:</strong> Define $F$ and $U$, state what their zeros identify, and
describe the behavior of $T_n$ under fixed alternatives.
</div>
<div class="solution-answer">
<strong>Solution:</strong>
Let $\mu_p$ and $V_p$ be the pooled population Fréchet mean and variance. The population
mean contrast is
$$
F=V_p-\sum_{j=1}^k\lambda_jV_j
=
\sum_{j=1}^k\lambda_j
\left[
\mathbb E_j d^2(\mu_p,Y_j)-\mathbb E_j d^2(\mu_j,Y_j)
\right].
$$
Every summand is nonnegative because $\mu_j$ minimizes the $j$th Fréchet function. Under
unique means, $F=0$ if and only if $\mu_1=\cdots=\mu_k$.
The variance contrast is
$$
U
=
\sum_{j<\ell}
\frac{\lambda_j\lambda_\ell}{\sigma_j^2\sigma_\ell^2}
(V_j-V_\ell)^2.
$$
Provided $\sigma_j^2>0$, $U=0$ if and only if $V_1=\cdots=V_k$.
Under the consistency conditions in @prp-anova-population-consistency,
$F_n\to F$ and the sample variance quantities converge to their population counterparts.
Thus a fixed mean alternative gives $F>0$, making the normalized mean component diverge.
A fixed variance alternative gives $U>0$, making the variance component diverge. In either
case the omnibus rejection probability tends to one.
</div>
</div>
</dialog>
<dialog id="l4-sol-4">
<div class="solution-dialog">
<div class="solution-dialog-header">
<h4>Exercise 4: Shrinking Alternatives</h4>
<button onclick="closeSolution('l4-sol-4')" class="solution-dialog-close">✕ Close</button>
</div>
<div class="solution-original">
<strong>Exercise:</strong> State the sufficient mean- and variance-separation rates under
shrinking alternatives and explain why they differ.
</div>
<div class="solution-answer">
<strong>Solution:</strong>
For alternative classes satisfying $F\ge a_n$ or $U\ge b_n$, with
$a_n,b_n\to0$, @thm-anova-power gives
$$
\sqrt n\,a_n\to\infty
\quad\Longrightarrow\quad
\beta_{\mathcal H_n^F}\to1,
$$
and
$$
n b_n\to\infty
\quad\Longrightarrow\quad
\beta_{\mathcal H_n^U}\to1.
$$
The difference comes from how the two contrasts enter $T_n$. The mean-sensitive term contains
$nF_n^2$, so a mean contrast must dominate $n^{-1/2}$. The variance component contains
$nU_n$, so a variance contrast must dominate $n^{-1}$. Fixed alternatives are immediate
special cases because a positive constant dominates either shrinking threshold.
</div>
</div>
</dialog>
<dialog id="l4-sol-5">
<div class="solution-dialog">
<div class="solution-dialog-header">
<h4>Exercise 5: SPD Simulation Design</h4>
<button onclick="closeSolution('l4-sol-5')" class="solution-dialog-close">✕ Close</button>
</div>
<div class="solution-original">
<strong>Exercise:</strong> Explain the non-overlapping block design, the exact null at
$s=1$, positive-definite interpolation, and the interpretation of rejection.
</div>
<div class="solution-answer">
<strong>Solution:</strong>
For every regime $g$ and object $i$, the simulation generates a fresh block of $L$ Gaussian
returns and computes one block sample covariance $Y_{gi}$. Distinct blocks share no returns,
so the covariance objects are independent within and across groups, as assumed by the
$k$-sample theory. Overlapping rolling windows would violate that structure and make the
nominal sample size misleading.
With $\tau=(s-1)/3$, the population covariance path is
$$
\Sigma_g(s)
=
\exp\{(1-\tau)\log\Sigma_0+\tau\log\Sigma_g^{\mathrm{target}}\}.
$$
At $s=1$, $\tau=0$ and every group covariance equals $\Sigma_0$, so both the mean and
variance nulls are true. Matrix exponentials of symmetric matrices are SPD, so the full path
remains positive definite without clipping individual correlation entries.
An omnibus rejection says that at least one regime Fréchet mean or Fréchet variance differs.
Because the objects are covariance matrices, the difference can arise from marginal
variances, correlations, or both. It does not establish a correlation-only difference; that
claim would require applying the analysis to correlation-matrix objects.
</div>
</div>
</dialog>
<script>
function showSolution(id) {
const dialog = document.getElementById(id);
if (dialog) {
dialog.showModal();
if (window.MathJax && MathJax.typesetPromise) {
MathJax.typesetPromise([dialog]).catch(function(err) {
console.log('MathJax typeset error:', err);
});
}
}
}
function closeSolution(id) {
const dialog = document.getElementById(id);
if (dialog) {
dialog.close();
}
}
document.addEventListener('click', function(e) {
if (e.target.tagName === 'DIALOG') {
e.target.close();
}
});
</script>
# Further Reading
- @DubeyMuller2019 — The Fréchet ANOVA framework: $k$-sample testing for equality of Fréchet means and variances.
- @DubeyMuller2020 — Fréchet change-point detection, building on the Fréchet variance CLT.
# Self-Assessment Quiz
Test your understanding of this lecture with the interactive MCQ quiz:
👉 **[Lecture 4 Quiz — 10 Multiple-Choice Questions](../quizzes/lecture-04-quiz.qmd)**