<figure><img src="image_1.png" /></figure>

HAL
open science

# Three-Cornered Hat and Groslambert Covariance: A First Attempt to Assess the Uncertainty Domains

François Vernotte, Eric Lantz

► To cite this version:

François Vernotte, Eric Lantz. Three-Cornered Hat and Groslambert Covariance: A First Attempt to Assess the Uncertainty Domains. IEEE Transactions on Ultrasonics, Ferroelectrics and Frequency Control, 2019, 66 (3), pp.643 - 653. 10.1109/TUFFC.2018.2889703. hal-02382643

HAL Id: hal-02382643

https://hal.science/hal-02382643v1

Submitted on 29 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Three-Cornered Hat and Groslambert Covariance: A first attempt to assess the uncertainty domains

François Vernotte and Éric Lantz

**Abstract**—The three-cornered hat method and the Groslambert Covariance are very often used to estimate the frequency stability of each individual oscillator in a set of three oscillators by comparing them in pairs. However, no rigorous method to assess the uncertainties over their estimates has yet been formulated. In order to overcome this lack, this paper will first study the direct problem, i.e. the calculation of the statistics of the clock stability estimates by assuming known values of the true clock stabilities and then will propose a first attempt to solve the inverse problem, i.e. the assessment of a confidence interval over the true clock stabilities by assuming known values of the clock stability estimates. We show that this method is reliable from 5 Equivalent Degrees of Freedom (EDF) and beyond.

## I. INTRODUCTION

The “Three-Cornered Hat method” was introduced by Gray and Allan in 1974 [1] to estimate the frequency stability of each individual oscillator in a set of three oscillators by comparing them in pairs. This method has been proved to be very efficient if the underlying assumption of the oscillator independence is fulfilled. Another approach, based on covariances, was proposed by Fest, Groslambert and Gagnepain in 1983 [2]. This method gives very similar results to such an extent that both methods were considered to be perfectly equivalent until very recently. However, we showed in 2016 that the latter, renamed the Groslambert Covariance (GCov), has the advantage to reject the noise measurement, i.e. the counter noise [3]. The use of GCov is becoming widespread and dedicated measuring instruments are beginning to appear [4], [5].

However, except approximated error bars valid only for small integration times [6] or limited guidelines [3], no rigorous method to assess confidence intervals over the three-cornered hat or GCov estimates has yet been formulated. Nevertheless, it is an important issue, especially since negative variance estimates may be obtained by these methods, because these estimates are computed by calculating differences. Moreover, the uncertainty over each clock stability estimate strongly depends on the stabilities of all the clocks.

The aim of this paper is to study the statistics of the three-cornered hat/GCov method in order to find a way for computing confidence intervals over the true clock stabilities. After a quick reminder of both three-cornered hat and GCov methods, we will first address the so-called “direct problem”, i.e. the

calculation of the probability distribution of the clock stability estimates by knowing the true clock stabilities. Then, we will propose a first attempt of solving the “inverse problem”, i.e. assessing a confidence interval over the true clock stabilities by knowing the clock stability estimates obtained by the three-cornered hat/GCov method.

## II. STATEMENT OF THE PROBLEM

### A. Clock comparison

Three-cornered hat as well as Groslambert Covariance rely both on simultaneous comparisons of 3 clocks A, B and C in pairs [3]. The clocks are assumed to be uncorrelated.

1) **Time and frequency quantities:** In this paper we refer to time and frequency quantities such as phase time x, fractional frequency y, z (see below)...of oscillators. We indicate the signals from oscillators with capital letters (A, B, C...) that are used as subscripts of the related quantities, i.e. yA is the fractional frequency of the oscillator A. The measurement number k can also be used as subscript.

2) **Variances and covariances:** AVAR is generally defined as

$$ \sigma_y^2 = \frac{1}{2} E \left[ (\bar{y}_{k+1} - \bar{y}_k)^2 \right] $$

where $E[\cdot]$ is the mathematical expectation, and $\bar{y}_k = \frac{1}{\tau} \int_{t_k}^{t_{k+1}} y(t) dt$ with $t_{k+1} = t_k + \tau$.

In order to simplify the notations, let us define the quantity $z_k = (\bar{y}_{k+1} - \bar{y}_k)/\sqrt{2}$. AVAR may then be written as $\sigma_y^2(\tau) = E[z_k^2]$. The associated estimator is then<fnref n="1" />

$$ \hat{\sigma}_y^2(\tau) = \frac{1}{M} \sum_{k=1}^{M} z_k^2 $$

where M is the number of different $z_k$ in a data run of length T, whether AVAR is calculated with or without overlapping.

3) **GCov vs 3-cornered hat:** The classical 3-cornered hat relies on the assumption of independence of the channel noises and of the oscillators A, B, C. The intercomparison of A and B is measured by $z_{AB} = z_B - z_A$ and its variance is $\sigma_{AB}^2 = \sigma_A^2 + \sigma_B^2$ where all dependences on $\tau$ are omitted. Similarly, $\sigma_{BC}^2 = \sigma_B^2 + \sigma_C^2$ and $\sigma_{CA}^2 = \sigma_C^2 + \sigma_A^2$. The three-cornered hat method uses the following property: $\sigma_A^2 = \frac{1}{2}(\sigma_{AB}^2 - \sigma_{BC}^2 + \sigma_{CA}^2)$.

On the other hand, the variance of oscillator A may be estimated by using the covariance of the inter-comparison of A and B as well as of the inter-comparison of A and

<footnote>F. Vernotte is with UTINAM, UMR 6213 CNRS, Observatory THETA, Université Bourgogne Franche-Comté, 41 bis avenue de l'Observatoire, 25010 Besançon - France</footnote>

<footnote>E. Lantz is with Institut FEMTO-ST, Département d'Optique P. M. Duffieux, UMR 6174 CNRS, Université Bourgogne Franche-Comté, 15b Avenue des Montboucons, 25030 Besançon - France</footnote>

<footnote n="1">¹In this paper, the symbol $\hat{\cdot}$ stands for the estimate of the quantity which is below.</footnote>
C. This leads to the Groslambert covariance $GCov_A = E[(z_A - z_B)(z_A - z_C)]$.

As mentioned above, these two approaches are almost equivalent, but $GCov$ is not polluted by the measurement noises since all cross-covariances are zero-mean [3]. However, in this paper we will not distinguish these two approaches and will use one or the other method for mathematical derivations.

## B. Definition of the problem

1) **Measurements, model parameters and estimates:** In the following, in order to simplify the notation, we will refer to $z_{PQ}$ and $\hat{\sigma}_{PQ}^2$ with $P, Q \in \{A, B, C\}$ to ensure the consistency of the notation regardless of the clock pair. The **"elementary estimates"** of the variance resulting from the comparison of the clocks $P$ and $Q$ is then

$$ \hat{\sigma}_{PQ}^2 = \frac{1}{M} \sum_{k=1}^{M} z_{PQ,k}^2 \quad (1) $$

Similarly, the Groslambert covariance estimates will be denoted by

$$ \widehat{GCov}_{PO,PQ} = \frac{1}{M} \sum_{k=1}^{M} z_{PO,k} \cdot z_{PQ,k} \quad (2) $$

where $O, P, Q \in \{A, B, C\}$ and $M$ is the number of different $z_{PO,k}$ (or $z_{PQ,k}$) in a data run of length $T$.

Finally, from these covariances or from the three-cornered hat, we can compute the **"final estimates"** $\hat{\sigma}_P^2$ with $P \in \{A, B, C\}$.

2) **Direct problem and inverse problem:** In order to assess the uncertainties over the estimation of the individual clock stabilities, we will have to distinguish two main issues:

* The **direct problem** consists in calculating the statistics of the elementary estimates $\hat{\sigma}_{PQ}^2$, or final estimates $\hat{\sigma}_P^2$, knowing the model parameters $\sigma_A^2, \sigma_B^2$ and $\sigma_C^2$.

* The **inverse problem**, conversely, consists in calculating a confidence interval over each model parameter $\sigma_P^2$, from the final estimates $\hat{\sigma}_A^2, \hat{\sigma}_B^2$ and $\hat{\sigma}_C^2$. Obviously, this last step is the true purpose of this paper.

This distinction corresponds to the main sections of this paper.

## III. DIRECT PROBLEM

### A. Calculation of two independent $\chi^2$ distributions

1) **Statistics of the $\hat{\sigma}_{PQ}^2$ estimates:** The two approaches described above are strictly equivalent when the counter noises are negligible. This condition occurs for large integration time $\tau$. We will then alternatively use the 3-cornered hat or the $GCov$ formalism for our demonstrations. Here we begin with the 3-cornered hat but the mathematical derivations below remain valid for the $GCov$ method regardless of $\tau$.

In the definition of the elementary estimates given by (1), the number $M$ is a key element since it determines the number of Equivalent Degrees of Freedom (EDF) of these random variables (r.v.). If we use AVAR "without overlapping" on a White FM noise, the EDF is simply $\nu = T/\tau - 1$ where $T$ is the total duration of the data-run. The use of Allan

Variance with overlapping or its application to other noises than White FM will describe an identical statistics but with a different number of EDF. Since the number of EDFs is an independent parameter in this study, it will suffice to enter the number of EDF corresponding to the case treated in the obtained relationships, regardless of the choice of the variance or the type of noise.

The $z_{PQ,k}$ measurements are Gaussian centered r.v. and Equation (1) shows that $\hat{\sigma}_{PQ}^2$ is a r.v. that follows a $\chi^2$ distribution with $\nu$ degrees of freedom, which we will write $\chi^2_\nu$. The number of EDFs is here such that $1 \le \nu \le M$ according to the correlation between the $z_{PQ,k}$ measurements. The precise determination of the number of EDFs does not fit into this study and is described in other publications (see for example [7]).

Therefore, the 3 estimates $\hat{\sigma}_{AB}^2, \hat{\sigma}_{BC}^2$ and $\hat{\sigma}_{CA}^2$ follow different $\chi^2$ distributions: $\hat{\sigma}_{AB}^2 = k_{AB}\hat{x}_\nu^2, \hat{\sigma}_{BC}^2 = k_{BC}\hat{x}_\nu^2$ and $\hat{\sigma}_{CA}^2 = k_{CA}\hat{x}_\nu^2$ where the coefficients $k_{AB}, k_{BC}, k_{CA} \in \mathbb{R}^+$. Since an oscillator is involved in two distributions, these distributions are correlated.

2) **Vector formalization of the problem:** To determine the statistics followed by the final estimates, we must write them in the form of linear combinations of independent r.v.. Therefore, we have to determine an orthonormal basis of $N$ linearly independent vectors from $N'$ linearly dependent vectors, with $N < N'$. In our case, $N' = 3$ (3 correlated $\chi^2$ r.v.) and we will show in the following that $N = 2$ (2 independent $\chi^2$ r.v.).

To simplify the problem, let us express (1) in the case $M = 1$ and thus $\nu = 1$ EDF: $\hat{\sigma}_{PQ}^2 = z_{PQ}^2$. We can also consider that $z_P = (\bar{y}_{P,2} - \bar{y}_{P,1})/\sqrt{2}$ and therefore $z_{PQ} = z_Q - z_P$. Each $z_P$ follows a centered normal (Gaussian or Laplace-Gauss) law $LG(0, \mathcal{Z}_P)$. Hence $z_{PQ} = z_Q - z_P$ follows a centered normal law of variance $\hat{\sigma}_{PQ}^2 = \hat{z}_P^2 + \hat{z}_Q^2$.

Since the 3 clocks are independent, the quantities $\mathcal{Z}_P$ can be considered as the 3 coordinates of a vector in a 3-dimensional space of basis $(\vec{e}_A, \vec{e}_B, \vec{e}_C)$. It is a vector space of normal laws since each of these normed basis vectors is characterized by a centered and reduced normal $LG(0, 1)$ independent in such a way that $E[\vec{e}_P \cdot \vec{e}_Q] = \delta_{P,Q}$, where $\delta_{P,Q}$ represents the Kronecker symbol.

This vector space is thus endowed with a scalar product, denoted ‘.’, and a norm denoted ‘$||\cdot||^2$’, defined by

$$ \left\{ \begin{array}{l} \vec{e}_P \cdot \vec{e}_Q = \mathrm{L}\mathrm{G}(0, 1) \cdot \mathrm{L}\mathrm{G}(0, 1) \quad \text{if } P \neq Q \\ \vec{e}_P \cdot \vec{e}_P = \|\vec{e}_P\|^2 = \chi_1^2 \end{array} \right. $$

where $\chi_1^2$ is a r.v. following a $\chi^2$ law with 1 degree of freedom and $\mathrm{L}\mathrm{G}(0, 1)$ and $\mathrm{L}\mathrm{G}(0, 1)$ represent 2 independent normal r.v. Their product therefore follows a Bessel distribution of mathematical expectation 0 and variance 1 [8], which is a special case of the variance-gamma distribution defined in Appendix A-1, with $A = B$ and $\nu = 1$, and then $\theta = 0, \eta = \lambda = \kappa = 1/2$ in (17).

Thus, the quantities $z_A, z_B$ and $z_C$ become vectors that we will write $\vec{A}, \vec{B}$ and $\vec{C}$ and which are defined by: $\vec{A} = \mathcal{Z}_A \vec{e}_A = (\mathcal{Z}_A, 0, 0)^T$, $\vec{B} = \mathcal{Z}_B \vec{e}_B = (0, \mathcal{Z}_B, 0)^T$, and $\vec{C} = \mathcal{Z}_C \vec{e}_C = (0, 0, \mathcal{Z}_C)^T$.

We can now express the variances and covariances of the primary oscillators through norms and scalar products of
vectors: $||\vec{B}||^2 = Z_B^2||\vec{e}_B||^2 = Z_B^2\chi_A^2$. We see in particular that $\sigma_B^2 = E[||\vec{B}||^2] = Z_B^2$. As for the scalar product between two different basis vectors, its expectation is null and equal to the covariance between two independent oscillators. However, it is described by a variance-gamma distribution: $\vec{A} \cdot \vec{B} = Z_A Z_B \vec{e}_A \cdot \vec{e}_B = Z_A Z_B V\Gamma(1, 0, 1/2)$.

The clock comparisons $z_{AB}$, $z_{BC}$ and $z_{CA}$ also become vectors formed according to the following model: $z_{AB} = z_B - z_A = \vec{A}\vec{B} = \vec{B} - \vec{A} = -Z_A\vec{e}_A + Z_B\vec{e}_B$.

Similarly, the variances $\hat{\sigma}_{AB}^2$, $\hat{\sigma}_{BC}^2$ and $\hat{\sigma}_{CA}^2$ can be rewritten with these notations:

$$ 
\begin{aligned}
\hat{\sigma}_{AB}^2 &= ||\overrightarrow{AB}||^2 = Z_A^2||\vec{e}_A||^2 + Z_B^2||\vec{e}_B||^2 \\
&\quad -2Z_A Z_B \vec{e}_A \cdot \vec{e}_B \\
&= Z_A^2 \dot{\chi}_1^2 + Z_B^2 \ddot{\chi}_1^2 - 2Z_A Z_B V\Gamma(1, 0, 1/2)
\end{aligned}
$$

It is now easy to see that the 3 vectors $\overrightarrow{AB}$, $\overrightarrow{BC}$ and $\overrightarrow{CA}$ are linearly dependent since $\overrightarrow{AB} + \overrightarrow{BC} + \overrightarrow{CA} = \vec{0}$.

Therefore, they all belong to a 2-dimensional subspace of the above defined 3-dimensional space. It is thus necessary to look for an orthonormal basis ($\vec{u}_1, \vec{u}_2$) in which we can rewrite these 3 vectors.

3) Search for an independent basis: The Gram-Schmidt algorithm can be used to solve this problem. Let us choose $\vec{u}_2 = \gamma \vec{C} \vec{A}$ in such a way that $E[||\vec{u}_2||^2] = 1$:

$$ ||\vec{u}_2||^2 = \gamma^2 Z_A^2 \dot{\chi}_1^2 + \gamma^2 Z_C^2 \ddot{\chi}_1^2 - 2\gamma^2 Z_A Z_C V\Gamma(1, 0, 1/2) $$

and then $E[||\vec{u}_2||^2] = \gamma^2 Z_A^2 + \gamma^2 Z_C^2 = 1$ where

$$ \gamma = \frac{1}{\sqrt{Z_A^2 + Z_C^2}}. \qquad (3) $$

We now need to find a second vector, $\vec{u}_1$, which is a linear combination of $\overrightarrow{AB}$ and $\overrightarrow{BC}$, perpendicular to $\vec{u}_2$, with a norm of unity mathematical expectation:

$$ 
\left\{ 
\begin{aligned}
\vec{u}_1 &= \alpha \overrightarrow{AB} - \beta \overrightarrow{BC} \\
E[\vec{u}_1 \cdot \vec{u}_2] &= -\alpha \gamma Z_A^2 + \beta \gamma Z_C^2 = 0 \\
E[||\vec{u}_1||^2] &= \alpha^2 Z_A^2 + (\alpha - \beta)^2 Z_B^2 + \beta^2 Z_C^2 = 1.
\end{aligned}
\right. 
$$

The solution is:

$$ 
\left\{ 
\begin{aligned}
\alpha &= \frac{Z_C^2}{\sqrt{Z_A^4(Z_B^2+Z_C^2)+Z_C^4(Z_A^2+Z_B^2)+2Z_A^2Z_B^2Z_C^2}} \\
\beta &= \frac{Z_A^2}{\sqrt{Z_A^4(Z_B^2+Z_C^2)+Z_C^4(Z_A^2+Z_B^2)+2Z_A^2Z_B^2Z_C^2}}
\end{aligned}
\right. 
\qquad (4) 
$$

We can now write the three vectors $\overrightarrow{AB}$, $\overrightarrow{BC}$ and $\overrightarrow{CA}$ as linear combinations of the basis vectors ($\vec{u}_1, \vec{u}_2$):

$$ \vec{C} \vec{A} = \frac{1}{\gamma} \vec{u}_2. \qquad (5) $$

Using the property $\overrightarrow{AB} + \overrightarrow{BC} + \overrightarrow{CA} = 0$, we obtain:

$$ 
\begin{aligned}
\vec{u}_1 &= \alpha \overrightarrow{AB} - \beta \overrightarrow{BC} = \alpha \overrightarrow{AB} - \beta (-\overrightarrow{AB} - \overrightarrow{CA}) \\
&= (\alpha + \beta) \overrightarrow{AB} + \frac{\beta}{\gamma} \vec{u}_2
\end{aligned}
$$

leading finally to:

$$ \overrightarrow{AB} = \frac{1}{\alpha + \beta} \vec{u}_1 - \frac{\beta}{\gamma(\alpha + \beta)} \vec{u}_2. \qquad (6) $$

Similarly, we find:

$$ \overrightarrow{BC} = -\frac{1}{\alpha + \beta} \vec{u}_1 - \frac{\alpha}{\gamma(\alpha + \beta)} \vec{u}_2. \qquad (7) $$

Thanks to the equations (5), (6) and (7), we can estimate $\sigma_B^2$, for example by using $\overline{\text{GRev}}_{BA,BC}$:

$$ 
\begin{aligned}
\overline{\text{GRev}}_{BA,BC} &= \overrightarrow{BA} \cdot \overrightarrow{BC} \\
&= \frac{1}{(\alpha + \beta)^2} \dot{\chi}_1^2 - \frac{\alpha\beta}{\gamma^2(\alpha + \beta)^2} \ddot{\chi}_1^2 \\
&\quad + \frac{\alpha - \beta}{\gamma(\alpha + \beta)^2} V\Gamma(1, 0, 1/2)
\end{aligned}
$$

where now the $2 \chi^2$ r.v. are totally independent. We note, however, that there is a third term of null mathematical expectation, but of which variance is non-zero except in the particular case where $\alpha = \beta$, i.e. $Z_A = Z_C$. In this case, we will have access to the distribution of $\hat{\sigma}_B^2$ since we know that the difference of two $\chi^2$ r.v. corresponds to a r.v. with a probability density given by the equation (17) of Appendix A-1.

In the general case $\alpha \neq \beta$, a Bessel distribution is added to the difference between the $\chi^2$ r.v., preventing us from calculating the density of $\hat{\sigma}_B^2$. This problem can be solved by rotating the basis ($\vec{u}_1, \vec{u}_2$) to get a new basis ($\vec{V}_1, \vec{V}_2$) so that the multiplicative factor of the scalar product $\vec{V}_1 \cdot \vec{V}_2$ is identically zero.

## B. Rotation of the basis vectors

It is thus necessary to find the eigenvector basis ($\vec{V}_1, \vec{V}_2$) such that

$$ 
\begin{aligned}
\hat{\sigma}_B^2 &= \frac{1}{(\alpha + \beta)^2} ||\vec{u}_1||^2 - \frac{\alpha\beta}{\gamma^2(\alpha + \beta)^2} ||\vec{u}_2||^2 \\
&\quad + \frac{\alpha - \beta}{\gamma(\alpha + \beta)^2} \vec{u}_1 \cdot \vec{u}_2 \\
&= \mathcal{A} ||\vec{V}_1||^2 - \mathcal{B} ||\vec{V}_2||^2 + \mathcal{C} \vec{V}_1 \cdot \vec{V}_2
\end{aligned}
\qquad (8) 
$$

with $C = 0$. Let us simplify the notation of Equation (8): $\hat{\sigma}_B^2 = a||\vec{u}_1||^2 - b||\vec{u}_2||^2 + c\vec{u}_1 \cdot \vec{u}_2$ with

$$ 
\left\{ 
\begin{aligned}
a &= \frac{1}{(\alpha + \beta)^2} \\
b &= \frac{\alpha\beta}{\gamma^2(\alpha + \beta)^2} \\
c &= \frac{\alpha - \beta}{\gamma(\alpha + \beta)^2}.
\end{aligned}
\right. 
\qquad (9) 
$$

It is a quadratic form which associates a scalar $w_0$ to any vector $\vec{w} = w_1\vec{u}_1 + w_2\vec{u}_2$ according to

$$ w_0 = \vec{w}^T Q \vec{w} \qquad (10) $$

with

$$ Q = \begin{pmatrix} a & c/2 \\ c/2 & -b \end{pmatrix}. $$

Diagonalizing the matrix $Q$ gives the eigenvalues:

$$ 
\left\{ 
\begin{aligned}
\mathcal{A} &= \frac{a-b+\sqrt{(a+b)^2+c^2}}{2} \\
\mathcal{B} &= \frac{a-b-\sqrt{(a+b)^2+c^2}}{2}.
\end{aligned}
\right. 
\qquad (11) 
$$
and eigenvectors $\vec{E}_A = (c/2, A - a)^T$ and $\vec{E}_B = (B + b, c/2)^T$. The eigenvector matrix is $E = (\vec{E}_A \vec{E}_B)$ and the quadratic form of (10) may be rewritten as $w_0 = \vec{w}^T E \Lambda E^T \vec{w}$. Thus, this new formulation can be interpreted as a basis change in which a vector $\vec{w}$ is transformed into a vector $\vec{W}$ according to:

$$ \vec{W} = E^T \vec{w}. \quad (12) $$

The basis change of relation (12) is a rotation of angle $\varphi$ since

$$ E^T = \begin{pmatrix} \cos \varphi & \sin \varphi \\ -\sin \varphi & \cos \varphi \end{pmatrix} $$

with

$$ \varphi = \arctan \left[ \frac{(-a - b + \sqrt{(a + b)^2 + c^2})}{c} \right]. \quad (13) $$

This angle, although not essential for the calculation of basis vectors, proves to be a very useful indicator to analyze the data (see Section III-C3).

### C. Validation of the theoretical probability laws by Monte Carlo simulations

1) **Validation principle:** According to the preceding section, the probability density of $\hat{\sigma}_B^2$, equal to the difference of two independent $\chi^2$ r.v., can now be calculated using the function $p(x)$ of the equation (17) of Appendix A-1 by assigning the following values to the parameters of this function:

$$ \left\{ \begin{array}{l} \eta = \frac{A+B}{4AB} \\ \theta = \frac{A-B}{4AB} \\ \kappa = \sqrt{\eta^2 - \theta^2} = \frac{1}{2\sqrt{AB}} \\ \lambda = \nu/2 \\ \mu = 0 \end{array} \right. \quad (14) $$

where $A$ and $B$ are the values given in (11) and $\nu$ the number of EDF of the considered $\chi^2$ laws. Indeed, if the above study was done for $\nu = 1$ for reasons of simplification of the formalism, it remains perfectly valid regardless of the number of EDFs.

To verify the results of this model, we compared them with those given by a simulation that seems realistic: we chose to return to the frequency deviations of the individual clocks by simulating the quantities $\bar{y}_{A,k}$, $\bar{y}_{B,k}$ and $\bar{y}_{C,k}$ (see § II-A2). They were simulated by a centered normal r.v. (randn function of Octave, the Matlab clone). It might be objected that we were simulating white noise while frequency deviations are much more likely to be red noise samples. But what matters in this study is the Gaussian character of the probability law more than the power law of its spectral density. The only effect of the latter is reducing the number of EDFs of the $\chi^2$ laws.

### – Simulation algorithm

The simulation algorithm follows these 6 steps

S1: Assignment of the 3 noise levels $Z_P$

S2: Drawing of $3 \times \nu$ samples $\bar{y}_{P,k}$

S3: Computation of the $3 \times \nu$ differences $\bar{y}_{PQ,k} = \bar{y}_{Q,k} - \bar{y}_{P,k}$

S4: Computation of the 3 estimates $\hat{\sigma}_P^2 = \text{GCov}_{PO,PQ}$

<figure><img src="image_3.png" /><figcaption>Fig. 1. Probability densities of the final estimate $\hat{\sigma}_B^2$ obtained by the model without rotation (solid lines) and by simulations (boxes) for comparable stability clocks: $\sigma_A^2 = \sigma_B^2 = \sigma_C^2 = 1$ (red); B. $\sigma_A^2 = 0.5$, $\sigma_B^2 = 1$, $\sigma_C^2 = 2$ (green). The chosen number of EDFs was $\nu = 5$.</figcaption></figure>

S5: Repetition $N = 10^7$ times of the steps S2 to S4 of this sequence.

S6: Drawing of the 3 histograms of the $\hat{\sigma}_P^2$.

In all simulations, we chose a number of EDF $\nu = 5$.

### – Modeling algorithm

The modeling algorithm follows these 6 steps for each estimate $\hat{\sigma}_P^2$

S1: Assignment of the noise levels of the 3 clocks $Z_A$, $Z_B$ and $Z_C$

S2: **Independent basis**

* Computation of the coefficients $\alpha, \beta$ and $\gamma$ according to the relationships (3) and (4)

* Computation of the coefficients $a, b$ and $c$ according to (9)

S3: **Vector rotation** (option)

* Computation of the rotation angle $\varphi$ with (13)

* Computation of the coefficients $A, B$ and $C$ with (11)

S4: Computation of the coefficients $\eta, \theta, \kappa$ and $\lambda$ with (14)

S5: Plotting the probability density with (17).

S6: Repeat the steps S2 to S5 for the other 2 clocks.

2) **Modeling without the basis rotation:** All the histograms shown in the following figures have been normalized by the total number of draws and the width of the bins to get the average probability density over each bin.

As expected, when the noise levels of the 3 clocks are the same, the match between the model and the simulations is excellent (figure 1 in red).

On the other hand, when the differences in noise level are significant ($\sigma_C^2 = 2\sigma_B^2 = 4\sigma_A^2$), the gap becomes obvious, especially with regard to the measurement $\hat{\sigma}_B^2$ (see Figure 1 in green). It is indeed for this one that the disparity between the 2 other levels is the most important ($\sigma_C^2 = 8\sigma_A^2$) and the approximation $A - B \ll A + B$, necessary for the model without rotation to be valid, is no longer verified.

3) **Modeling with the rotation:** Now, let's put the model with the rotation of the basis vectors on the test bench by
<table><caption>TABLE I
COMPARISON OF 3 CLOCKS OF COMPARABLE STABILITIES USING THE ROTATION OF THE EIGENVECTORS. THE CHOSEN NUMBER OF EDFs WAS $\nu = 5$.</caption><thead><tr><th rowspan="2">Angle</th><th colspan="2">$\sigma_A^2 = 0.1$</th><th colspan="2">$\sigma_B^2 = 1$</th><th colspan="2">$\sigma_C^2 = 10$</th></tr><tr><th>Model</th><th>Simul.</th><th>Model</th><th>Simul.</th><th>Model</th><th>Simul.</th></tr></thead><tbody><tr><td>Fractile</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>2.5 %</td><td>-2.894</td><td>-2.893</td><td>-1.773</td><td>-1.775</td><td>1.428</td><td>1.428</td></tr><tr><td>97.5 %</td><td>3.190</td><td>3.193</td><td>4.715</td><td>4.715</td><td>26.09</td><td>26.08</td></tr><tr><td>$P(\hat{\sigma} < 0)$</td><td>47.5 %</td><td>47.4 %</td><td>26.6 %</td><td>26.6 %</td><td>0.06 %</td><td>0.06 %</td></tr></tbody></table>

comparing it to the simulations.

### Clocks with comparable stabilities

At first, we took again the levels of stability tested with the model without rotation (see Section III-C2) and we found in each case an excellent agreement between model and simulation.

We can treat with the same success the case where there is a factor 100 between $\sigma_A^2$ and $\sigma_C^2$ (see Table I and Figure 2). To be able to compare the results given by the model and the simulations, we give them with 4 significant digits in Table I. Whether it be for the fractiles (lines 4 and 5) or for the probability of getting a negative measurement (line 6), we can notice that the results agree at least up to the significant 3rd digit.

We obtain measurements that are almost symmetrically distributed around 0 for $\hat{\sigma}_A^2$, meaning that these measurements are masked by the measurement fluctuations of the other clocks. Therefore, the stability of this clock is not measurable using the two clocks B and C which are much less stable. On the other hand, the stability of the least stable clock, $\sigma_C^2$, is relatively well determined (let us remind that with only 5 EDFs, the confidence intervals remain nevertheless very wide) since the rate of negative measurements is only 6 out of 10 000. Finally, the determination of the stability $\sigma_B^2$ is intermediate between those of $\sigma_A^2$ and of $\sigma_C^2$.

<figure><img src="image_4.png" /><figcaption>Fig. 2. Probability densities of the final estimate $\hat{\sigma}_A^2$ (red), $\hat{\sigma}_B^2$ (green), $\hat{\sigma}_C^2$ (blue), obtained by the model with rotation (solid lines) and by the simulations (boxes). The stability of the clock ensemble was given by the following parameters: $\sigma_A^2 = 0.1$, $\sigma_B^2 = 1$, $\sigma_C^2 = 10$. The chosen number of EDFs was $\nu = 5$.</figcaption></figure>

### One clock is much less stable or much more stable than the others

We then examined the case where two of the clocks are of identical stabilities and the other is either 100 times more stable or 100 times less stable. Here also we found a full agreement between model and simulation.

To summarize, the model with the rotation of the basis vectors fits perfectly whatever the noise levels of the clocks. Therefore, we can consider that we have solved the direct problem.

## IV. INVERSE PROBLEM

### A. Principle of the method

1) Inverse problem and parameter uncertainties: The metrologist has to solve the inverse problem, i.e. the determination of a confidence interval for the true variances $\sigma_A^2$, $\sigma_B^2$ and $\sigma_C^2$, given a set of measurements and a priori information, i.e. any information known before the measurements. In this experimental world, the true variances appear as random variables, of which the a posteriori probability density, i.e. probability density that takes into account the measurements, is determined by using the Bayes theorem:

$$ \left\{ \begin{array}{l} p(\sigma_A^2, \sigma_B^2, \sigma_C^2 | \hat{\sigma}_A^2, \hat{\sigma}_B^2, \hat{\sigma}_C^2) \\ \propto \pi(\sigma_A^2, \sigma_B^2, \sigma_C^2) p(\hat{\sigma}_A^2, \hat{\sigma}_B^2, \hat{\sigma}_C^2 | \sigma_A^2, \sigma_B^2, \sigma_C^2) \\ \iiint_0^\infty p(\sigma_A^2, \sigma_B^2, \sigma_C^2 | \hat{\sigma}_A^2, \hat{\sigma}_B^2, \hat{\sigma}_C^2) d\sigma_A^2 d\sigma_B^2 d\sigma_C^2 = 1 \end{array} \right. \quad (15) $$

where $\pi(\theta)$ is the a priori probability density, named prior, of a value $\theta$ [9], here the variances. Even if we have no a priori information on the variances, $\pi(\sigma_A^2, \sigma_B^2, \sigma_C^2)$ can be defined: a variance is a positive scale parameter[10], [11] and a prior reflecting no a priori knowledge (total ignorance) is proportional to $\frac{1}{\theta}$, meaning that all orders of magnitudes have the same a priori probability. Because the oscillators are assumed to be independent, the prior of a triplet of variances is simply given by the product of the individual priors:

$$ \pi(\sigma_A^2, \sigma_B^2, \sigma_C^2) \propto \frac{1}{\sigma_A^2 \sigma_B^2 \sigma_C^2} $$

To compute the a posteriori probability of a variance triplet given by Eq. (15), it remains to calculate the probability $p(\hat{\sigma}_A^2, \hat{\sigma}_B^2, \hat{\sigma}_C^2 | \sigma_A^2, \sigma_B^2, \sigma_C^2)$. This is a direct problem, that has been solved in the first part of this paper for a single estimate, for example $p(\hat{\sigma}_A^2 | \sigma_A^2, \sigma_B^2, \sigma_C^2)$. Unfortunately, the three final estimates $\hat{\sigma}_A^2$, $\hat{\sigma}_B^2$, $\hat{\sigma}_C^2$ are not independent and the probability of a triplet of estimate is not given by the product of the probabilities determined in the preceding section. This issue cannot be solved by using the elementary estimates defined in section (II-B1): $\hat{\sigma}_{AB}^2$, $\hat{\sigma}_{AC}^2$, $\hat{\sigma}_{BC}^2$, since they are neither independent (each oscillator participates in two of them). Before making probability products, we must find linear combinations of our three estimates that are independent, or, at least, uncorrelated. It is well known that these combinations are obtained by applying the Karhunen-Loève (K.L.) transform to our original estimates, either elementary or final [12]. In the
following lines, the principle of the K.L transform is recalled and applied to our specific case. Since we are supposed to know at this step the true variances, it is possible to calculate the true covariance matrix of our estimates (either $\hat{\sigma}_A^2, \hat{\sigma}_B^2, \hat{\sigma}_C^2$ or $\hat{\sigma}_{AB}^2, \hat{\sigma}_{AC}^2, \hat{\sigma}_{BC}^2$), and to find the rotation that renders diagonal this covariance matrix. The calculation of this true covariance matrix is detailed in Section A-2.1. The coefficients of the eigenvectors of this diagonalisation process are used as weighting numbers to compute a new triplet of uncorrelated estimates. Actually, these K.L estimates do not obey Gaussian statistics, since the original estimates follow the not trivial probability density law exposed in the preceding sections. However, we will approximate in the following this distribution by a Gaussian one. If we assume that this approximation is correct, the three K.L. estimates obey each Gaussian statistics, with their three variances given by the diagonalized covariance matrix. Moreover, in the frame of this approximation, the K.L. estimates are independent and the probability density of a K.L. triplet is simply given by the product of the three probability densities. Because of the one to one correspondence between the K.L. triplet and the original estimates, this probability density is also proportional to the probability density of the triplet of original estimates. We expect that this Gaussian approximation becomes more accurate for a large number of measurements (more EDF). We will see in the following that this is correct, but the approximation can be used even for a small number of measurements, at least 5 (see Section IV-B3), at the price of a reasonable inaccuracy in the limits of the confidence intervals. This inaccuracy will be assessed in Section IV-B.

2) Algorithm: We have to calculate the triple integral of Eq. (15) on several order of magnitudes. A direct calculation would lead to prohibitive computation times. We have preferred to use a Monte-Carlo scheme with random sampling. This sampling ensures the observance of the total ignorance a priori law: the samples are chosen at random on a logarithmic scale in such a way that all orders of magnitude have the same probability (see the concept of total ignorance in §IV-A1), independently for each variance. With a computation on 8 decades (between $10^{-5}$ and $10^3$) with $10^4$ sampling steps, $10^7$ samples proved to be sufficient to render negligible the sampling error, in comparison with other inaccuracies. With the same sampling step, a direct calculation would involve $10^{12}$ cells.

We work in the experimenter point of view: we assume that a triplet of estimates (either final or elementary) has been calculated from the 3 $\cdot m$ elementary measurements. These three numbers have three definite values that will be used in the calculations detailed below. The different steps of the calculation can be summarized as

S1: Choose at random a triplet of true variances, with a uniform probability on a logarithmic scale for each variance and independence between the three variances.

S2: Calculate for this triplet the covariance matrix of the estimate triplet, either using the final or the elementary estimates.

S3: Calculate the eigenvectors and eigenvalues of this covariance matrix

S4: Multiply the vector of the estimates by the matrix of these eigenvectors (K.L. transform).

S5: Perform the same operation for the vector of true variances.

S6: Calculate the probability density of each K.L estimate given the triplet of K.L. true variances: for each of the three K.L. variables, we assume a Gaussian normal law of mean the K.L true variance and of variance the corresponding eigenvalue of the covariance matrix.

S7: Perform the product of these three probability densities. This only number will be associated in the following to the triplet of variances chosen at the first step of the algorithm.

S8: Repeat $10^7$ times the entire process.

S9: For each of the three variables, normalize the probability densities by dividing by their sum (sum of $10^7$ values).

S10: Also for each of the three variables, sort the true variance values and calculate the cumulative density function by a partial sum on the associated normalized probability densities.

S11: Determine a confidence interval at 95% on each true variance from the corresponding cumulative density functions.

S12: Verify that the low limit of the confidence interval is meaningful. For a Gaussian distribution, 99.7% of data are included in a confidence interval at $\pm 3\sigma$. If the low limit of this $\pm 3\sigma$ confidence interval (in logarithmic scale) is smaller than the low limit of the a priori range (here $10^{-5}$), we suspect (and have verified) that the low limit of the smaller confidence interval calculated in the preceding step will depend on the low limit of the a priori range. If it occurs, we replace the low limit of the confidence interval by 0.

We have verified that employing the elementary or the final estimates gives exactly the same results. Moreover, the uncertainties due to the random character of the Monte-Carlo integration are less than 1% in relative value. Hence, the only non negligible cause of error in the algorithm is the Gaussian approximation. This error is assessed in the following paragraph (see Section IV-B). But before, it is useful to find the expected properties of the final estimates.

3) Properties of the final estimates: – General properties

P1: Only one final estimate may be negative.

Demonstration: All elementary estimates are positive $\hat{\sigma}_{AB}^2 = \hat{\sigma}_A^2 + \hat{\sigma}_B^2 > 0$ and $\hat{\sigma}_{BC}^2 = \hat{\sigma}_B^2 + \hat{\sigma}_C^2 > 0$ and $\hat{\sigma}_{CA}^2 = \hat{\sigma}_C^2 + \hat{\sigma}_A^2 > 0$. Therefore, if two final estimates would be negative, at least one of the elementary estimate would be negative and this is impossible.

P2: If a final estimate is negative, its absolute value is smaller than the absolute values of the other two final estimates of the triplet.

Demonstration: if $\hat{\sigma}_A^2 < 0$, $\hat{\sigma}_{AB}^2 = \hat{\sigma}_A^2 + \hat{\sigma}_B^2 = -|\hat{\sigma}_A^2| + |\hat{\sigma}_B^2| > 0$ ⇒ $|\hat{\sigma}_B^2| > |\hat{\sigma}_A^2|$. Similarly, $\hat{\sigma}_{CA}^2 = \hat{\sigma}_C^2 + \hat{\sigma}_A^2 = |\hat{\sigma}_C^2| - |\hat{\sigma}_A^2| > 0$ ⇒ $|\hat{\sigma}_C^2| > |\hat{\sigma}_A^2|$.

– Case of $\nu = 1$ EDF
P3: $\hat{\sigma}_P^2 = -\frac{\hat{\sigma}_O^2 \hat{\sigma}_Q^2}{\hat{\sigma}_O^2 + \hat{\sigma}_Q^2}$ with 1 EDF and $\{O, P, Q\}$ any circular permutation of $\{A, B, C\}$.

Demonstration: With 1 EDF, Equation (2) becomes: $\hat{\sigma}_P^2 = -z_{OP} \cdot z_{PQ} = -(z_P - z_O)(z_Q - z_P)$. Similarly, with 1 EDF Equation (1) becomes: $\hat{\sigma}_{PQ}^2 = z_{PQ}^2 = (z_Q - z_P)^2$ Moreover $\hat{\sigma}_{PQ}^2 = \hat{\sigma}_P^2 + \hat{\sigma}_Q^2$. Therefore

$$ 
\begin{aligned}
-\frac{\hat{\sigma}_O^2 \hat{\sigma}_Q^2}{\hat{\sigma}_O^2 + \hat{\sigma}_Q^2} &= -\frac{\hat{\sigma}_O^2 \hat{\sigma}_Q^2}{\hat{\sigma}_{OQ}^2} \\
&= -\frac{(z_O - z_Q)(z_O - z_P)(z_Q - z_P)(z_Q - z_O)}{(z_Q - z_O)^2} \\
&= -(z_P - z_O)(z_Q - z_P) = \hat{\sigma}_P^2.
\end{aligned}
$$

This latter property is very important since the case of $\nu = 1$ EDF occurs for the largest integration time, i.e. $\tau = T/2$ for AVAR and a dataset of duration $T$. For this largest $\tau$, Property 3 implies two major consequences:

* one of the final estimates of the triplet and only one is necessarily negative

* the knowledge of two final estimates of the triplet is enough since the third one can be deduced from them by using Property 3 (e.g. if $\hat{\sigma}_B^2 = \hat{\sigma}_C^2 = 1$ then $\hat{\sigma}_A^2 = -1/2$).

## B. Accuracy of the inverse algorithm

1) **Principle of the simulations:** In order to assess the accuracy of the method we propose, we compared it to Monte-Carlo simulations. The principle consists in randomly drawing parameter triplets $(\sigma_A^2, \sigma_B^2, \sigma_C^2)$, computing the corresponding final estimate triplet $(\hat{\sigma}_A^2, \hat{\sigma}_B^2, \hat{\sigma}_C^2)$ and keeping only the parameter triplets which yield a previously given final estimate triplet $(\hat{\sigma}_A^2, \hat{\sigma}_B^2, \hat{\sigma}_C^2) = (A_0, B_0, C_0)$. Obviously, the random generation is of importance: as previously, we choose each element of each triplet independently according to an uniform probability law on a logarithmic scale. The pseudorandom numbers may vary within a large interval depending on the $(A_0, B_0, C_0)$ triplet and of the EDF, typically between $B_l = 10^{-3}$ and $B_h = 10^{+2}$ for $(1, 1, 1)$. Each simulation run stops when 10,000 achievements have been obtained.

2) **Preliminary observations on the simulations: - Case of several EDF**

Before using the simulation results to assess the method accuracy, let us observe them in a few cases: keeping $\hat{\sigma}_B^2 = \hat{\sigma}_C^2 = 1$, we will successively vary $\hat{\sigma}_A^2 \in \{0.1, 1, 10\}$, i.e. corresponding to one final estimate lower, equal to or higher than the other two ones. Figure 3 shows the histograms of the parameter $\sigma_A^2$ in each of these cases with 20 EDF.

$\hat{\sigma}_A^2 = 0.1$: The histogram of the corresponding $\sigma_A^2$ parameter (see Figure 3 in red) is constant for the lower values, exhibits a very slight bulge between 0.1 and 1, and tends to 0 after 1. Obviously, the histogram is limited to $10^{-4}$ at its left-hand side because we limited the random generation to $B_l = 10^{-4}$ but this trend should continue down to $\sigma_A^2 = 0$. The lower bound of the confidence interval should then be 0.

<figure><img src="image_5.png" /><figcaption>Fig. 3. Normalized histograms of 10000 $\sigma_A^2$ parameters giving the final estimate triplets (0.1, 1, 1) in red, (1,1,1) in green, (10,1,1) in blue. The number of EDF is 20.</figcaption></figure>

$\hat{\sigma}_A^2 = 1$: The histogram plotted in green in Figure 3 shows an important bump between 0.1 and 1 and tends to 0 after 1. Nevertheless, there is also a constant tail, although much lower than in the previous case, for the $\sigma_A^2$ values below 0.1 down to 0. Here also this tail is limited by the lower bound of the random generation ($B_l = 10^{-3}$). The lower bound of the confidence interval is then still 0.

$\hat{\sigma}_A^2 = 10$: In this case (see Figure 3 in blue), $\sigma_A^2$ is well constrained around 10 and the histogram seems to be almost Gaussian. No doubt that the 95 % confidence interval will be defined for the lower and the upper bound.

Another representation is given in Figure 4. This 3D plot was built by associating a dimension to each parameter of the triplets from the dataset obtained with a final estimate triplet equal to (1, 1, 1) (the same dataset as for Figure 3). It has the advantage of showing the relationships between these parameters. This type of plot exhibits a structure with 3 perpendicular branches, more or less dense depending on the final estimate triplet, converging to (1, 1, 1).

<figure><img src="image_6.png" /><figcaption>Fig. 4. 3D-plot of 10000 parameter triplets giving the final estimate triplets (1,1,1). The number of EDF is 20. The solid red lines represent respectively the lines $[\sigma_A^2 = \hat{\sigma}_A^2 \& \sigma_B^2 = \hat{\sigma}_B^2]$, $[\sigma_B^2 = \hat{\sigma}_B^2 \& \sigma_C^2 = \hat{\sigma}_C^2]$ and $[\sigma_C^2 = \hat{\sigma}_C^2 \& \sigma_A^2 = \hat{\sigma}_A^2]$.</figcaption></figure>
### 3) Influence of the number of EDF
In order to check the validity of our method to assess the confidence intervals over the parameters, let us first assign the final estimate triplet to (1, 1, 1) and vary the number of EDF.

#### - number of EDF varying from 2 to 1000

Figure 5 shows a pretty good agreement between the 95 % confidence intervals estimated by our method (colored area) and the one obtained from a set of 10 000 Monte-Carlo simulations (error bars), except for $v = 2$ EDF where our method seriously underestimates the 97.5 % bound. On the other hand, the method is fully reliable from $v = 5$ EDF and above.

For $v$ varying from 2 to 20, there is also a noticeable discrepancy between the 2.5 % bound obtained by our method and by simulations. However, the 2.5 % bound of the simulation error bars is almost exactly the lower bound $B_l$ which limits the pseudo-random excursion of our simulations. Therefore, this 2.5 % bound is only due to a computational artifact and the 0 result of our method is more reliable.

For $v = 50$, our method seems to slightly overestimate the confidence interval, but for higher EDF the Monte-Carlo simulations give the same confidence intervals as our method. This result is not surprising since the Gaussian approximation we adopted in our algorithm is perfectly justified for large EDFs.

#### - Case of 1 EDF

Table II shows that the 97.5 % bound obtained by our method is underestimated by a factor of approximately 100! It should then not be used for $v \le 2$.

### 4) Influence of relative values of the final estimates: - Low EDF: $v = 5$

Figure 6 shows a very good agreement between our method and the simulations for the parameter $\hat{\sigma}_A^2$ (in red). There is only one slight underestimation of the 2.5 % bound for $\hat{\sigma}_A^2 = 4$. Similarly, for the parameter $\hat{\sigma}_B^2$ or $\hat{\sigma}_C^2$ (see Figure 6 in

<figure><img src="image_7.png" /><figcaption>Fig. 5. Influence of the number of EDF $v$ on the 95 % confidence interval over the $\hat{\sigma}_A^2$ parameter. The final estimate triplet was set at $\hat{\sigma}_A^2 = \hat{\sigma}_B^2 = \hat{\sigma}_C^2 = 1$. The colored area corresponds to the assessment of the confidence interval by our method whereas the error bars correspond to the empirical estimation from 10 000 Monte-Carlo simulations.</figcaption></figure>

<caption>TABLE II</caption>

CASE OF 1 EDF - COMPARISON OF THE 95 % CONFIDENCE INTERVAL BOUNDS OBTAINED BY OUR METHOD (HEREAFTER LV METHOD STANDING FOR LANTZ-VERNOTTE METHOD) AND BY 10 000 MONTE-CARLO SIMULATIONS. THE FINAL ESTIMATE TRIPLET WAS SET AT $(-1/2, 1, 1)$.

<table><thead><tr><th rowspan="2">Parameter</th><th colspan="2">2.5 % bound</th><th colspan="2">97.5 % bound</th></tr><tr><th>LV method</th><th>Simulation</th><th>LV method</th><th>Simulation</th></tr></thead><tbody><tr><td>$\hat{\sigma}_A^2$</td><td>0</td><td>2.01 · 10<sup>-6</sup></td><td>1.39</td><td>159</td></tr><tr><td>$\hat{\sigma}_B^2$</td><td>0</td><td>3.96 · 10<sup>-6</sup></td><td>5.28</td><td>474</td></tr><tr><td>$\hat{\sigma}_C^2$</td><td>0</td><td>3.50 · 10<sup>-6</sup></td><td>5.31</td><td>464</td></tr></tbody></table>

<figure><img src="image_8.png" /><figcaption>Fig. 6. Influence of the final estimate $\hat{\sigma}_A^2$ on the 95 % confidence interval over the $\hat{\sigma}_A^2$ parameter (red) and over the $\hat{\sigma}_B^2$ or $\hat{\sigma}_C^2$ parameters (green). The $\hat{\sigma}_B^2$ and $\hat{\sigma}_C^2$ final estimates were set at 1 and the number of EDF is 5. The colored area corresponds to the assessment of the confidence interval by our method whereas the error bars correspond to the empirical estimation from 10 000 Monte-Carlo simulations. The graphs are plotted using a log-log plot.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 7. Influence of the final estimate $\hat{\sigma}_A^2$ on the 95 % confidence interval over the $\hat{\sigma}_A^2$ parameter (red) and over the $\hat{\sigma}_B^2$ or $\hat{\sigma}_C^2$ parameters (green). The $\hat{\sigma}_B^2$ and $\hat{\sigma}_C^2$ final estimates were set at 1 and the EDF is 5. The colored area corresponds to the assessment of the confidence interval by our method whereas the error bars correspond to the empirical estimation from 10 000 Monte-Carlo simulations. The graphs are plotted using a linear X-scale.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Fig. 8. Influence of the final estimate $\hat{\sigma}_A^2$ on the 95 % confidence interval over the $\sigma_A^2$ parameter (red) and over the $\sigma_B^2$ or $\sigma_C^2$ parameters (green). The $\hat{\sigma}_B^2$ and $\hat{\sigma}_C^2$ final estimates was set at 1 and the EDF is 100. The colored area corresponds to the assessment of the confidence interval by our method whereas the error bars correspond to the empirical estimation from 10 000 Monte-Carlo simulations.</figcaption></figure>

green), there is a good agreement between our method and the simulations. However, our method shows a 2.5 % bound which increases for a final estimate $\hat{\sigma}_A^2 < 0.1$ which is, unexpectedly, not confirmed by the simulations.

Figure 7 shows the same parameters versus the same final estimate $\hat{\sigma}_A^2$ but by using a linear X-scale allowing us to plot the negative values of the final estimate. For $\sigma_A^2$ (red plot on Figure 7), our method gives good results down to $\hat{\sigma}_A^2 = -0.3$ but, surprisingly, the 97.5 % bound increases for $\hat{\sigma}_A^2 < -0.3$. Here also, it is not confirmed by the simulations. The same odd behavior may be observed for $\sigma_B^2$ and $\sigma_C^2$ (green plot on Figure 7) for the 2.5 % bound around 0 as well as for the 97.5 % bound below $-0.3$.

These discrepancies between our method and the simulations are due to the difference between the true $\chi^2$ distribution with a small number of EDF and the normal law approximation in step S6 of the algorithm (see §IV-A2). However, it does not really matter since, on one hand, these discrepancies are limited, and on the other hand, they occur only for final estimates close to the low limit (-0.5, see Section IV-A3) and are therefore rare.

$$- \text{High EDF: } \nu = 100$$

Figures 8 shows that our method works perfectly for 100 EDF. The only noticeable discrepancies concern the 2.5 % bound and is, once again, clearly due to the computational artifact of the simulation error bars mentioned above.

## V. CONCLUSION

We have performed a thorough theoretical study of the statistics of the 3 cornered hat or Groslambert covariance estimates (direct problem) showing that they are formed by random variables which are the differences of $2 \chi^2$ r.v. Massive simulations have been performed and the agreement between these simulations and the results of the theoretical study is very convincing.

We have also proposed a first attempt to assess a confidence interval over the true clock stabilities (inverse problem) relying upon the determination of the posterior probability density function by a Monte-Carlo computation. Here also, we have compared the results of this method to massive simulations. This study shows that this method may be used from EDF larger than 2 but is fully reliable from 5 EDF and beyond. This means that for a data run of duration $T$, the uncertainty domain is valid for an integration time $\tau \le T/10$ in the case of a white FM noise and below for flicker or random walk FM.

On the other hand, this method is relatively slow since its algorithm relies on a Monte-Carlo scheme involving $10^7$ random draws. This causes a computation time of the order of 1 minute per confidence interval (per error bar).

However, since the analysis of the direct problem presented in this paper seems to be convincingly achieved, it could constitute the basis of a new method of solving the inverse problem, i.e. of estimating confidence intervals around the true stability of the clocks knowing the results obtained by the three cornered hat/GCov method.

## ACKNOWLEDGMENT

This work was partially funded by the ANR Programme d'Investissement d'Avenir (PIA) under the Oscillator IMP project and the FIRST-TF network. The authors are very grateful to Professor Enrico Rubiola for suggesting us the subject of this work.

## REFERENCES

[1] J. Gray and D. Allan, "A method for estimating the frequency stability of an individual oscillator," in Proceedings of the 28th Annual Symposium on Frequency Control, May 1974, pp. 243-246.

[2] D. Fest, J. Groslambert, and J. Gagnepain, "Individual characterization of an oscillator by means of cross-correlation or cross-variance method," IEEE Transactions on Instrumentation and Measurement, vol. 32, no. 3, pp. 447-450, September 1983.

[3] F. Vernotte, C. E. Calosso, and E. Rubiola, "Three-Cornered Hat versus Allan Covariance," in IEEE International Frequency Control Symposium (IFCS), May 2016, pp. 1-6.

[4] M. Caligaris, G. Costanzo, and C. Calosso, "6/12-channel synchronous digital phasemeter for ultrastable signal characterization and use," in Proc. of Joint European Frequency and Time Forum and International Frequency Control Symposium, Dever (Co, USA), April 2015.

[5] C. E. Calosso, F. Vernotte, V. Giordano, C. Fluhr, B. Dubois, and E. Rubiola, "Frequency stability measurement of cryogenic sapphire oscillators with a multichannel tracking DDS and the two-sample covariance," IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, 2018 (in press).

[6] C. R. Ekstrom and P. A. Koppang, "Error bars for three-cornered hats," IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, vol. 53, no. 5, pp. 876-879, 2006.

[7] C. Greenhall and W. Riley, "Uncertainty of stability variances based on finite differences," in 35th annual Precise Time and Time Interval meeting, San Diego (California, USA), December 2003, pp. 267-280.

[8] B. Sorin and P. Thionet, "Lois de probabilités de bessel," Revue de statistique appliquée, vol. 16, no. 4, pp. 65-72, 1968.

[9] J. Bernardo and A. Smith, Bayesian Theory, ser. Wiley Series in Probability and Statistics. Wiley & Sons, 1994.

[10] D. V. Lindley, "Fiducial distributions and Bayes' theorem," Journal of the Royal Statistical Society, vol. 20, no. 1, pp. 102-107, 1958.

[11] E. Lantz and F. Vernotte, "Can we define a best estimator in simple 1-D cases?" IEEE Signal Processing Magazine, vol. 30, no. 6, pp. 151-156, 2013.

[12] F. Vernotte, D. Charraut, E. Lantz, and D. A. Courjon, "Karhunen-loève transform in optical microscopy: application to linewidth measurement," Optical Engineering, vol. 29, 1990.
# APPENDIX

## A. Difference of two random variables following two independent $\chi^2$ distributions

1) **Definition:** We consider a random variable $X$ which is the difference of two independent random variables following each a $\chi^2$ law with the same number of degrees of freedom $\nu$:

$$X = \mathcal{A}\chi_{\nu}^{2} - \mathcal{B}\ddot{\chi}_{\nu}^{2} \quad (16)$$

with $\mathcal{A}, \mathcal{B}, \nu \in \mathbb{R}^+$ and $\nu \ge 1$. The upper dots of $\dot{\chi}_{\nu}^2$ and $\ddot{\chi}_{\nu}^2$ make it possible to distinguish 2 independent random variables following the same $\chi^2$ distribution.

The probability density function (PDF) $p(x)$ of the random variable $X$ is the variance-gamma distribution:

$$p(x) = \frac{\kappa^{2\lambda}|x - \mu|^{\lambda-1/2} K_{\lambda-1/2}(\eta|x - \mu|)}{\sqrt{\pi\Gamma(\lambda)(2\eta)^{\lambda-1/2}}} e^{\theta(x-\mu)} \quad (17)$$

with $\eta = (\mathcal{A} + \mathcal{B})/4\mathcal{A}\mathcal{B}$, $\theta = (\mathcal{A} - \mathcal{B})/4\mathcal{A}\mathcal{B}$, $\lambda = \nu/2$, $\mu = 0$, $\kappa = \sqrt{\eta^2 - \theta^2}$ and $K_\omega(z)$ is a hyperbolic Bessel function of second kind ($\omega \in \mathbb{R}$ and $z \in \mathbb{C})^2$.

2) **Simulation:** The random variable $X$ following the distribution defined in (16) was simulated with $\mathcal{A} = 1$, $\mathcal{B} = 1/3$ and $\nu = 5$:

$$X = \dot{\chi}_{5}^{2} - \frac{1}{3}\ddot{\chi}_{5}^{2}.$$

A number of $N = 10\,000\,000$ random draws was realized.

The probability density expressed in (17) was compared to the histogram obtained from $N$ draws. The agreement is almost perfect.

Finally, the CDF fractiles were calculated for 2.5% and 97.5% to achieve a 95% confidence interval. Here too, the agreement is excellent.

A similar study was conducted for the $(\mathcal{A}, \mathcal{B}, \nu)$ triplets in $\{(1, 1/3, 5), (1, 1/2, 1), (1, 1/4, 1), (2, 2, 1), (1, 1, 1)\}$. The results were equally concordant.

## B. Calculation of the true covariance matrix of the estimates

We suppose in this appendix that the true variances are known. With this hypothesis of the model world, the estimates, either final or elementary, become random variables, with means equal either to the true variances (final estimates) or to the sum of two variances (elementary estimates). We calculate in this appendix their covariance matrix, that is used to solve the inverse problem.

1) **Covariance matrix of the elementary estimates:** From now on, we will denote $E(\cdot)$, $\text{Var}(\cdot)$ and $\text{Cov}(\cdot)$ respectively the mathematical expectation, the variance and the covariance of the quantity $\cdot$ which stands here for $z$ measurements as well as elementary or final estimates (variance of variances!).

To be specific, we calculate in the following a diagonal element, $\text{Var}(\hat{\sigma}_{AB}^2)$, and a non diagonal element, $\text{Cov}(\hat{\sigma}_{AB}^2, \hat{\sigma}_{AC}^2)$. We assume that successive measurements of the same quantity are independent, meaning that the variances and covariances after $m$ measurements are equal to the

(co)variances after one measurement divided by $m$. Hence, the following calculations are presented with $m = 1$. For a diagonal element, we obtain:

$$\begin{aligned} \text{Var}(\hat{\sigma}_{AB}^2) &= E\left[(z_A - z_B)^4\right] - \{E\left[(z_A - z_B)^2\right]\}^2 \\ &= \text{Var}(z_A^2) + \text{Var}(z_B^2) + 4\text{Var}(z_A) \cdot \text{Var}(z_B) \end{aligned}$$

We assume that $z_A$ is centered Gaussian, meaning that $z_A^2$ follows a $\chi^2$ law of mean $\text{Var}(z_A)$. For such a law, we have $\text{Var}(z_A^2) = 2[\text{Var}(z_A)]^2$, which allows a more compact formulation of the above result:

$$\begin{aligned} \text{Var}(\hat{\sigma}_{AB}^2) &= 2[\text{Var}(z_A)]^2 + 2[\text{Var}(z_B)]^2 \\ & \quad +4\text{Var}(z_A) \cdot \text{Var}(z_B) \\ &= 2[\text{Var}(z_A) + \text{Var}(z_B)]^2. \end{aligned}$$

We pass now to a non diagonal element:

$$\begin{aligned} \text{Cov}(\hat{\sigma}_{AB}^2, \hat{\sigma}_{AC}^2) &= E\left[(z_A - z_B)^2(z_A - z_C)^2\right] \\ & \quad -E\left[(z_A - z_B)^2\right] \cdot E\left[(z_A - z_C)^2\right] \\ &= 2[\text{Var}(z_A)]^2. \end{aligned}$$

2) **Covariance matrix of the final estimates:** We use the Grandjean covariance definition of the final estimates, $\hat{\sigma}_P^2 = \text{G}\text{Cov}_{PO,PQ}$ and calculate, as above, the covariance of these estimates: $\text{Cov}(\hat{\sigma}_P^2, \hat{\sigma}_Q^2)$. We assume no measurement noise. As in the previous subsection, we consider the case $m = 1$ measurement and we give explicit clock symbols to $P$ and $Q$: first, we calculate the diagonal elements of the estimates covariance matrix by choosing $P = Q = A$, then the non diagonal elements by choosing $P = A, Q = B$. In the first case, we obtain:

$$\begin{aligned} \text{Var}(\hat{\sigma}_A^2) &= E\left\{[(z_A - z_B)(z_A - z_C)]^2\right\} \\ & \quad - \{E\left[(z_A - z_B)(z_A - z_C)\right]\}^2 \\ &= \text{Var}(z_A^2) + \text{Var}(z_B) \cdot \text{Var}(z_C) \\ & \quad + \text{Var}(z_A) \cdot \text{Var}(z_B) + \text{Var}(z_A) \cdot \text{Var}(z_C) \end{aligned}$$

The non diagonal case gives:

$$\begin{aligned} \text{Cov}(\hat{\sigma}_A^2, \hat{\sigma}_B^2) &= E\left[(z_A - z_B)(z_A - z_C) \cdot (z_B - z_A)(z_B - z_C)\right] \\ & \quad -E\left[(z_A - z_B)(z_A - z_C)\right] \cdot E\left[(z_B - z_A)(z_B - z_C)\right] \\ &= \text{Var}(z_A) \cdot \text{Var}(z_B) \\ & \quad -\text{Var}(z_C)[\text{Var}(z_A) + \text{Var}(z_B)] \end{aligned}$$

In both cases, we have used the independence of the oscillators, giving: $E(z_P^2 z_Q^2) = E(z_P^2) \cdot E(z_Q^2)$ and their zero mean, giving $\text{Var}(z_P) = E(z_P^2)$.

Using the elementary or the final estimates gives exactly the same results, for the same measurements and the same Monte-Carlo set.

<footnote n="2">²References : https://math.stackexchange.com/questions/85249/distribution-of-difference-of-chi-squared-variables and https://en.wikipedia.org/wiki/Variance-gamma_distribution. The relationships between the $\eta, \theta$ coefficients and the $\mathcal{A}, \mathcal{B}$ parameters have been empirically determined.</footnote>