# Statistical Biases and Very-Long-Term Time Stability Analysis

F. Vernotte, E. Lantz

► To cite this version:

F. Vernotte, E. Lantz. Statistical Biases and Very-Long-Term Time Stability Analysis. IEEE Transactions on Ultrasonics, Ferroelectrics and Frequency Control, 2012, 59 (3), pp.523-530. 10.1109/TUFFC.2012.2223. hal-00708070

HAL Id: hal-00708070

https://hal.science/hal-00708070

Submitted on 29 Apr 2021

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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Statistical Biases and Very-Long-Term Time Stability Analysis

François Vernotte and Éric Lantz

**Abstract**—The prediction of very-long-term time stability is a key issue in various fields, such as time keeping, obviously, but also navigation and spatial applications. This is usually performed by extrapolating the measurement data obtained by estimators such as the Allan variance, modified Allan variance, Hadamard variance, etc.

This extrapolation may be assessed from a fit over the variance estimates. However, this fit should be performed on the log-log graph of the estimates, which corresponds to a least-squares minimization of the relative difference between the variance estimates and the fitting curve. However, a bias exists between the average of the log of the estimates and the log of the true value of the estimated variance.

This paper presents the theoretical calculation of this log-log bias based on the number of equivalent degrees of freedom of the estimates, shows simulations over a large number of realizations, and provides a reliable method of unbiased logarithmic fit. Extrapolating this fit yields a more confident assessment of the very-long-term time stability.

## I. INTRODUCTION

TIME stability of clocks and oscillators is generally assessed by estimators such as the Allan deviation, modified Allan deviation, Hadamard deviation, etc. However, unlike the variance estimates, the deviation estimates are severely biased and should not be used for further calculations, such as fit.

On the other hand, the prediction of very-long-term time stability is usually performed by fitting and extrapolating an Allan variance (AVAR) curve. Such a curve may be modeled as

$$ \sigma_y^2(\tau) = \sum_{i=0}^{4} C_i \Phi_i(\tau), \quad \text{with } \Phi_i(\tau) = \tau^{i-2}. \quad (1) $$

An extrapolation may be achieved by extending the asymptotes of such a curve beyond the largest $\tau$ value (see Fig. 1).

The $C_i$ coefficient are related to the noise levels $h_\alpha$ of the power law noise model by Table I, [1]:

$$ S_y(f) = \sum_{\alpha=-2}^{+2} h_\alpha f^\alpha, \quad (2) $$

F. Vernotte is with UTINAM, Observatory Terre-Homme-Environnement-Temps-Astronomie (THETA) of Franche-Comté, University of Franche-Comté/Centre National de la Recherche Scientifique (CNRS), Besançon, France (e-mail: francois.vernotte@obs-besancon.fr).

E. Lantz is with Department of Optics, FEMTO-ST, University of Franche-Comté/CNRS, Unité de Formation et de Recherche (UFR) Sciences, Besançon, France.

and to the linear frequency drift $D_1$:

$$ y(t) = D_0 + D_1 t + \eta(t), \quad (3) $$

where $\eta(t)$ represents the pure random component of the frequency deviation $y(t)$.

However such an extrapolation is very sensitive to estimation errors. It is then of importance to perform a reliable estimation of the $C_i$ coefficients. The weighted relative least squares are an effective method [2] which is weighted, to take into account the uncertainty $\Delta_j$ over each estimate $\hat{\sigma}_y^2(\tau_j)$, and relative, to be equivalent to a classical least square fit on the log-log plot by minimizing the relative distance between the curve and the estimate $[(\text{estimate} - \text{curve})/\text{estimate}]$,

$$ \sum_{j=1}^{N} \left[ \frac{1}{\Delta_j} \times \frac{\hat{\sigma}_y^2(\tau_j) - \sum_{i=0}^{4} C_i \Phi_i(\tau_j)}{\hat{\sigma}_y^2(\tau_j)} \right]^2 \quad \text{minimum.} \quad (4) $$

It is important to emphasize that for the relative fit, i.e., the fit of $\Delta\hat{\sigma}^2/\hat{\sigma}^2$ and not of $\Delta\hat{\sigma}^2$, that 1) it ensures that the fit will minimize the relative differences between the curve and the estimates (e.g., 10% for an estimate at the $10^{-14}$ level as well as 10% for an estimate at the $10^9$ level), whereas a classical fit would minimize the absolute differences, overwhelming the small estimates; and 2) because $dx/x = d\ln x$, it is exactly equivalent to a classical fit on the log of the estimates.

As a consequence, this extrapolation method is based on the logarithm of the estimates and not the estimates themselves. We know that the estimate $\hat{\sigma}_y^2(\tau)$ is an unbiased estimator of $\sigma_y^2(\tau)$, but it is not known if $\ln(\hat{\sigma}_y^2(\tau))$ is an unbiased estimator of $\ln(\sigma_y^2(\tau))$.

## II. STATISTICAL APPROACH

### A. Statistics of the Allan Variance and the Allan Deviation

The definition of the Allan variance is [3], [4]

$$ \sigma_y^2(\tau) = \frac{1}{2} \langle (\bar{y}_2 - \bar{y}_1)^2 \rangle. \quad (5) $$

The estimate of the Allan variance may be calculated as
<table><caption>TABLE I. THEORETICAL RESPONSES OF THE ALLAN VARIANCE FOR THE DIFFERENT NOISE TYPES OF THE POWER LAW MODEL AND FOR A LINEAR FREQUENCY DRIFT.</caption><thead><tr><th>S<sub>y</sub>(f)</th><th>h<sub>+2</sub>f<sup>+2</sup></th><th>h<sub>+1</sub>f<sup>+1</sup></th><th>h<sub>0</sub>f<sup>0</sup></th><th>h<sub>-1</sub>f<sup>-1</sup></th><th>h<sub>-2</sub>f<sup>-2</sup></th><th>y(t) = D<sub>1</sub>t</th></tr></thead><tbody><tr><th>σ<sub>y</sub><sup>2</sup>(τ)</th><td><sup>3h<sub>+2</sub>f<sub>h</sub></sup>&frasl;<sub>4π<sup>2</sup>τ<sup>2</sup></sub></td><td><sup>[1.038 + 3ln(2πf<sub>h</sub>τ)]h<sub>+1</sub></sup>&frasl;<sub>4π<sup>2</sup>τ<sup>2</sup></sub></td><td><sup>h<sub>0</sub></sup>&frasl;<sub>2τ</sub></td><td>2ln(2)h<sub>-1</sub></td><td><sup>2π<sup>2</sup>h<sub>-2</sub>τ</sup>&frasl;<sub>3</sub></td><td><sup>D<sub>1</sub><sup>2</sup>τ<sup>2</sup></sup>&frasl;<sub>2</sub></td></tr></tbody></table>

For the high frequency noises ($f^{+2}$ FM and $f^{+1}$ FM), it is necessary to introduce a high cut-off frequency $f_h$ [10].

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

Fig. 1. Extrapolation of a least square fit. The graph shows the Allan deviation (i.e., the square root of the Allan variance) but the fit was performed over the Allan variance estimates according to the criterion (4) using SigmaTheta-1.2 [15]. The equation of the curve is $[(4.7 \cdot 10^{-13})^2\tau^{-2} + (3.1 \cdot 10^{-13})^2\tau^{-1} + (4.8 \cdot 10^{-14})^2 + (1.2 \cdot 10^{-15})^2\tau^{+1} + (5.2 \cdot 10^{-18})^2\tau^{+2}]^{1/2}$.

$$ \hat{\sigma}_{y}^{2}(\tau) = \frac{1}{2N} \sum_{i=1}^{N} (\bar{y}_{2} - \bar{y}_{1})^{2}. \qquad (6) $$

Because $\bar{y}_2$ and $\bar{y}_1$ are Gaussian values<fnref n="1" />, $\bar{y}_2 - \bar{y}_1$ is a Gaussian centered value, so $(\bar{y}_2 - \bar{y}_1)^2$ obeys a $\chi_1^2$ distribution; therefore, $1/(2N) \sum_{i=1}^N (\bar{y}_2 - \bar{y}_1)^2$ is a $\chi_\nu^2$ distribution with $\nu = N$ if the $\bar{y}_i$ samples are independent, with $\nu < N$ otherwise. $\nu$ is the number of equivalent degrees of freedom (EDF).

The $\chi_\nu^2$ distributions have the following properties: $E(\chi_\nu^2) = \nu$, and $\text{var}(\chi_\nu^2) = 2\nu$ ($\Rightarrow \text{std}(\chi_\nu^2) = \sqrt{2\nu}$).

As shown by Fig. 2, the number of EDF reflects the AVAR estimates dispersion. A method for assessing the EDF may be found in [6].

## B. Model World and Measurement World

Let us define [see Fig. 3(a)] $\theta$ as the model parameter and $\xi$ as the measured data of the parameter. For example, we could consider the true value $\sigma_y^2(\tau = 10 \text{ s}) = h_0/20$.

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

Fig. 2. Dispersion of Allan variance (AVAR) estimates for a pure white frequency noise and their equivalent degrees of freedom, $\nu$ (65 536 samples at 1 s).

(where $h_0$ is the white FM level) as a parameter; $\hat{\sigma}_y^2(\tau)$ is a measurement data of $\sigma_y^2(\tau = 10 \text{ s})$.

Two problems are generally addressed: the direct problem, which attempts to forecast the measurement data knowing the parameter and the inverse problem, which aims to estimate the parameter knowing the measurement data. In the same vein, Tarantola distinguishes the model space, i.e., the space in which the parameter is given, from the data space, i.e., the space in which the measurement data are given [7]. In the following, we will use the terms model world and measurement world.

1) Model World (Direct Problem): In the direct problem, the task is to determine how the measurement data $\xi$ are distributed when the parameter $\theta_0$ is known [see Fig. 3(b)]. However, the parameter $\theta_0$ is precisely the unknown quantity that we want to estimate. This parameter is only known in theory and simulations.

2) Measurement World (Inverse Problem): In the inverse problem, the task is to estimate a confidence interval over $\theta$ from the known measurement data $\xi_0$ [see Fig. 3(c)]. This is the usual task of the metrologist.

## C. Study of a $\chi^2$ Distribution With $\nu = 1$ Degree of Freedom

1) Properties: The properties of a $\chi_1^2$ distribution are discussed in this section [8]. The probability density function (PDF) is

<footnote n="1">¹It is generally accepted that the frequency deviation samples follow statistics that can be approximated by a Gaussian model [5], although the non-stationary effects of low-frequency noises could limit the relevance of this model.</footnote>
<figure><img src="image_4.png" /><figcaption>Fig. 3. (a) Model parameter and measurement data for a $\chi_1^2$ distribution ($10^6$ realizations of $\theta \cdot \chi_1^2$ for $\theta = 10^k$ and $-3 < k < +9$). (b) Measurement dispersion for a parameter $\theta_0 \approx 1$. (c) Model parameter values for a measurement data $\xi_0 \approx 1$.</figcaption></figure>

$$p(X) = \frac{e^{-X/2}}{\sqrt{2\pi X}}. \qquad (7)$$

The PDF is normalized:

$$\int_{0}^{\infty} p(X)dX = 1. \qquad (8)$$

The mathematical expectation is

$$E(X) = \int_{0}^{\infty} X \cdot p(X)dX = \nu = 1. \qquad (9)$$

We can define the logarithmic mean as the exponential of the mathematical expectation of $\ln(X)$

$$\begin{aligned} e^{E(\ln(X))} &= \exp\left[\int_{0}^{\infty} \ln(X) \cdot p(X)dX\right] \\ &= \frac{e^{-C}}{2} \approx 0.2807, \end{aligned} \qquad (10)$$

where $C$ is the Euler constant; $C \approx 0.57722$. The cumulative distribution function (CDF) is

$$P(X) = \int_{0}^{X} p(y)dy = 1 - \frac{\Gamma_{\text{inc}}(1/2, X/2)}{\sqrt{\pi}}, \qquad (11)$$

and the inverse CDF is

$$P^{-1}(\alpha) = 2\Gamma_{\text{inc}}^{-1}(1/2, \sqrt{\pi}(1-\alpha)). \qquad (12)$$

This relationship is useful for getting the bounds of a confidence interval.

How is the reduced variable $X$ of the relationships (7) to (11) related to the measurement data $\xi$ and the model parameter $\theta$?

2) Parameter, Measurement, and Reduced Variable: Let us consider the standard $\chi_1^2$ variable $X = x^2$ where $x$ is a Gaussian centered standard random variable:

$$E(X) = \nu = 1 \Rightarrow E\left(\frac{X}{\nu}\right) = 1. \qquad (13)$$

It is well known that $\hat{\sigma}_y^2(\tau) = \xi$ is $\chi_1^2$ distributed; $\xi$ is an unbiased estimator of the parameter $\sigma_y^2(\tau) = \theta$. Therefore,

$$E\left(\frac{\xi}{\theta}\right) = 1. \qquad (14)$$

We can then define the reduced variable $X$ as

$$X = \nu \frac{\xi}{\theta} = \frac{\xi}{\theta}. \qquad (15)$$

The differential $dX$ is then

$$dX = \frac{\partial X}{\partial \xi} d\xi + \frac{\partial X}{\partial \theta} d\theta. \qquad (16)$$

3) Conditional Probabilities in the Model World (Direct Problem): The parameter $\theta$ is a constant. Let $\theta = \theta_0$; the differential $dX$ is then:
$$dX = \frac{\nu}{\theta_0} d\xi = \frac{1}{\theta_0} d\xi. \qquad (17)$$

From (17) and the relationships (7) to (9), we find the conditional probability density function:

$$p(\xi | \theta_0) = \frac{e^{-\xi/(2\theta_0)}}{\sqrt{2\pi\theta_0\xi}}. \qquad (18)$$

$\xi$ is unbiased:

$$E(\xi) = \int_{0}^{\infty} \xi \cdot p(\xi | \theta_0) d\xi = \theta_0. \qquad (19)$$

The logarithmic mean is

$$\begin{aligned} e^{E[\ln(\xi)]} &= \exp\left[\int_{0}^{\infty} \ln(\xi) \cdot p(\xi | \theta_0) d\xi\right] \qquad (20) \\ &= \frac{e^{-C}}{2} \theta_0 \approx 0.2807 \cdot \theta_0. \end{aligned}$$

Thus, the mathematical expectation of the log of the measurement data does not converge to the log of the parameter. Therefore, there is a significant bias for log-log representation (and fit); let us call it the log-log bias.

### 4) Conditional Probabilities in the Measurement World (Inverse Problem)

The measurement data $\xi$ is a constant. Let $\xi = \xi_0$; the differential $dX$ is then

$$dX = - \frac{\xi_0}{\theta^2} d\theta. \qquad (21)$$

From (21) and the relationships (7)–(9) comes the conditional probability density function:

$$p(\theta | \xi_0) = \sqrt{\frac{\xi_0}{2\pi\theta^3}} e^{-\xi_0/(2\theta)}. \qquad (22)$$

There is no mean value,

$$E(\theta) = \int_{0}^{\infty} \theta \cdot p(\theta | \xi_0) d\theta \qquad (23)$$

does not converge. The logarithmic mean is

$$\begin{aligned} e^{E[\ln(\theta)]} &= \exp\left[\int_{0}^{\infty} \ln(\theta) \cdot p(\theta | \xi_0) d\theta\right] \qquad (24) \\ &= 2e^C \xi_0 \approx 3.562 \cdot \xi_0. \end{aligned}$$

In the measurement world also, the mathematical expectation of the log of the parameter does not converge to the log of the measurement data. However, this bias is exactly the inverse of the log-log bias defined previously. (This property will be demonstrated later in a more general context, see Section II-D.) This implies that we just have to renormalize the measurement data by the log-log bias to get both an unbiased estimator in the model world

and a mathematical expectation of the parameter equal to the measurement data in the measurement world.

## D. Generalization to Any Positive-Valued Random Process

Let $\xi$ be a positive-valued estimator constructed with $M$ measurement data resulting from a random measurement process. $\xi$ is itself a random variable in the model world, which we assume to have a mean. The most obvious example is Allan variance. Another example is the estimation of the standard deviation. It is quite generally possible to define a new reduced variable,

$$Y = \frac{k\xi}{\theta_0}, \qquad (25)$$

which obeys (13), unlike $X$:

$$E(Y) = \int_{0}^{\infty} Y \cdot p(Y) dY = 1. \qquad (26)$$

Let $Z = f(Y)$ be a monotonic function of $Y$. A fundamental property of density probabilities is

$$p_Y(Y)dY = p_Z(Z)dZ, \qquad (27)$$

where the subscripts emphasize the fact that the different density probabilities are different functions (these subscripts are, as usual, omitted in the following). An immediate consequence is, for $Z = \xi$,

$$\int_{0}^{\infty} \xi \cdot p(\xi) d\xi = \theta_{0}/k. \qquad (28)$$

If $X$ is a standard $\chi_M^2$ variable with $\nu = M$ degrees of freedom, then $Y = X/\nu$; with $\xi$ a variance measurement data, we have, in the model world, $k = 1$ and (28) means that $\xi$ is an unbiased estimator of the true variance $\theta_0$. In the case of standard deviation estimation, the square root of the Allan variance is a biased estimator, but the value of $k$ can be easily calculated as a function of $M$ (the estimator is convergent: $k \to 1$ if $M \to \infty$) and (28) allows also the definition of an unbiased estimator.

The situation is worse in the measurement world. The variable is, in this case, $\theta$:

$$Y = \frac{k\xi_0}{\theta}, \qquad (29)$$

with, once more,

$$p(Y)dY = p(\theta | \xi_0)d\theta. \qquad (30)$$

Eqs. (29) and (30) give, for $k = 1$,

$$\int_{0}^{\infty} \theta \cdot p(\theta | \xi_0) d\theta = \xi_0 \cdot \int_{0}^{\infty} \frac{1}{Y} \cdot p(Y) dY \neq \xi_0. \qquad (31)$$

We have seen previously that this latter expression may not be convergent.
Let $Z = \ln(Y)$. We can easily calculate

$$ \int_{-\infty}^{\infty} Z \cdot p(Z) dZ = \int_{0}^{\infty} \ln(Y)p(Y)dY = B. \quad (32) $$

Hence, going back into the model world for $k=1$ and using (8):

$$ \begin{aligned} \int_{-\infty}^{\infty} \ln(\xi) \cdot p(\ln(\xi)) d(\ln(\xi)) &= \int_{0}^{\infty} \ln(\xi) p(Y) dY \quad (33) \\ &= B + \ln(\theta_0). \end{aligned} $$

Eq. (33) means that $\ln(\xi) - B$ is an unbiased estimator of $\ln(\theta_0)$. Hence, the log-unbiased estimator of $\theta_0$ appears to be, in the model world, $\xi'$, defined as

$$ \xi' = \xi \cdot e^{-B}. \quad (34) $$

In the measurement world, we have, from (29) with $k=1$,

$$ Y = \frac{\xi'_{0} \cdot e^{B}}{\theta}, \quad (35) $$

and then

$$ \ln(\theta) = \ln(\xi'_{0}) - \ln(Y) + B. \quad (36) $$

Unlike linear estimators, this estimator ensures that its mean is equal to the actual measurement data:

$$ \begin{aligned} E(\ln(\theta)) &= \int_{-\infty}^{\infty} \ln(\theta) \cdot p(\ln(\theta) | \xi'_{0}) d(\ln(\theta)) \quad (37) \\ &= \ln(\xi'_{0}) - B + B = \ln(\xi'_{0}). \end{aligned} $$

To conclude this section, a log-estimator can easily be rendered unbiased by multiplying the linear estimator by the right constant, with the important further advantage that, in the measurement world, the unknown true value has, for its mean, the actual measurement data. This equivalence between the model and measurement worlds for a scale parameter was already noticed 50 years ago [9].

## E. Generalization to a $\chi^2$ Distribution

The results of Section II-C can be extended to any $\chi^2$ distribution, regardless of the number of degrees of freedom $\nu$.

1) Model World: The reduced variable is

$$ X = \nu \frac{\xi}{\theta_0}, \quad (38) $$

and $E(X) = \nu$. The PDF is

$$ p(\xi | \theta_0) = \frac{\nu^{\nu/2} \xi^{\nu/2-1} e^{-(\nu\xi)/2\theta_0}}{2^{\nu/2} \Gamma(\nu/2) \theta_0^{\nu/2}}. \quad (39) $$

The mean is

$$ E(\xi) = \theta_0. \quad (40) $$

The logarithmic mean is

$$ e^{E[\ln(\xi)]} = \frac{2}{\nu} \exp\left[\psi^{(0)}\left(\frac{\nu}{2}\right)\right] \theta_0, \quad (41) $$

where $\psi^{(0)}(x)$ is the digamma function [10].

2) Measurement World: The reduced variable is

$$ X = \nu \frac{\xi_0}{\theta}, \quad (42) $$

and $E(X) = \nu$. The PDF is

$$ p(\theta | \xi_0) = \frac{\nu^{\nu/2} \xi_0^{\nu/2} e^{-(\nu\xi_0)/(2\theta)}}{2^{\nu/2} \Gamma(\nu/2) \theta^{\nu/2+1}}. \quad (43) $$

The mean ($\nu > 2$) is

$$ E(\theta) = \frac{\nu}{\nu - 2} \xi_0. \quad (44) $$

The logarithmic mean is

$$ e^{E[\ln(\theta)]} = \frac{\nu}{2} \exp\left[-\psi^{(0)}\left(\frac{\nu}{2}\right)\right] \xi_0. \quad (45) $$

Table II and Fig. 4 show the biases and the 95% confidence interval versus EDF.

## III. DISCUSSION

### A. A log-Unbiased Estimator

We know that $\hat{\sigma}_y^2(\tau)$ is an unbiased estimator of $\sigma_y^2(\tau)$. However, $E[\ln(\xi)] \neq \ln(\theta_0)$ and $e^{E[\ln(\xi)]} \neq \theta_0$. Similarly, $E[\ln(\theta)] \neq \ln(\xi_0)$ and $e^{E[\ln(\theta)]} \neq \xi_0$.

<table><caption>TABLE II. BIASES AND 95% CONFIDENCE INTERVAL BOUNDS VERSUS EDF $\nu$.</caption><thead><tr><th>ν</th><th>E(θ)/ξ<sub>0</sub></th><th>e<sup>E[ln(θ)]</sup>/ξ<sub>0</sub></th><th>B<sub>2.5%</sub></th><th>B<sub>97.5%</sub></th></tr></thead><tbody><tr><td>1</td><td>∞</td><td>3.56213</td><td>0.19905</td><td>1018.26086</td></tr><tr><td>2</td><td>∞</td><td>1.78108</td><td>0.27108</td><td>39.49783</td></tr><tr><td>3</td><td>3.00000</td><td>1.44625</td><td>0.32091</td><td>13.90209</td></tr><tr><td>4</td><td>2.00000</td><td>1.31043</td><td>0.35896</td><td>8.25735</td></tr><tr><td>5</td><td>1.66667</td><td>1.23754</td><td>0.38964</td><td>6.01530</td></tr><tr><td>6</td><td>1.50000</td><td>1.19222</td><td>0.41524</td><td>4.84911</td></tr><tr><td>7</td><td>1.40000</td><td>1.16137</td><td>0.43715</td><td>4.14232</td></tr><tr><td>8</td><td>1.33333</td><td>1.13902</td><td>0.45624</td><td>3.67017</td></tr><tr><td>9</td><td>1.28571</td><td>1.12210</td><td>0.47312</td><td>3.33285</td></tr><tr><td>10</td><td>1.25000</td><td>1.10885</td><td>0.48821</td><td>3.07978</td></tr></tbody></table>
<figure><img src="image_5.png" /><figcaption>Fig. 4. Biases and 95% confidence intervals over the model parameter $\theta$ versus equivalent degrees of freedom, $\nu$, in the measurement world where the measurement data are $\xi_0 = 1$. The circles represent the exponential of the mathematical expectation of the log of the model parameter and the crosses represent the mathematical expectation of the model parameter (which does not converge for $\nu \le 2$).</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 5. Example of log-log fit (obtained with SigmaTheta-1.2 [15]). As in Fig. 1, the graph shows the Allan deviation (ADEV) but the fit was performed over the Allan variance (AVAR) estimates. The crosses represent the ADEV estimates, whereas the circles represent the square root of the log-unbiased estimates.</figcaption></figure>

Let us define a new estimator $\hat{\mathcal{L}}_{\theta}$ as

$$ \hat{\mathcal{L}}_{\theta} = \ln\left\{\frac{\nu}{2} \exp\left[-\psi^{(0)}\left(\frac{\nu}{2}\right)\right]\xi\right\}, \quad (46) $$

and, with $\nu = 1$,

$$ \hat{\mathcal{L}}_{\theta} = \ln(3.5625 \cdot \xi) = \ln(\xi) + 1.2705. \quad (47) $$

The mathematical expectation of $\hat{\mathcal{L}}_{\theta}$ is then

$$ \langle \hat{\mathcal{L}}_{\theta} \rangle = \ln(\theta_0). \quad (48) $$

$\hat{\mathcal{L}}_{\theta}$ is then an unbiased estimator of $\ln(\theta_0)$.

On the other hand, from (45), for measurement data $\xi_0$:

$$ E[\ln(\theta)] = \ln\left\{\frac{\nu}{2} \exp\left[-\psi^{(0)}\left(\frac{\nu}{2}\right)\right]\xi_0\right\} = \hat{\mathcal{L}}_{\theta}, \quad (49) $$

and from (24), with $\nu = 1$,

$$ E[\ln(\theta)] = \ln(3.5625 \cdot \xi_0) = \hat{\mathcal{L}}_{\theta}. \quad (50) $$

The mathematical expectation of $\ln(\theta)$ is equal to $\hat{\mathcal{L}}_{\theta}$.

Therefore, $e^{\hat{\mathcal{L}}_{\theta}}$ is a log-unbiased estimator of $\theta$. For fitting a log-log plot, this log-unbiased estimator is then optimal.

## B. Fitting an Allan Variance Curve

Fig. 5 displays the bias estimates (crosses) as well as the log-unbiased estimates (circles) of a simulated oscillator. Fitting the AVAR estimates rather than the log-unbiased estimates leads to a significant underestimation of the long-term asymptotes ($\tau$ and $\tau^2$ asymptotes) and therefore of the very-long-term stability of the oscillator.

## C. Simulation Results

1) Pure Random-Walk FM: We first simulated pure random-walk FM, because it is often the dominating noise for the very-long-term.

a) Principle of the simulation: We performed one million simulations of a sequence of 1024 frequency deviation samples of pure random-walk FM with a noise level (model parameter) within the interval: $10^{-3} < h_{-2} < 10^{+9}$ with a uniform distribution of the exponent<fnref n="2" /> between -3 and +9. Thus, this variation range can cover all of the model parameters which may give a measurement data $\hat{h}_{-2}$ equal to 1. It is then possible to select either the realizations for which the model parameter $h_{-2}$ is close to 1 at 5% (i.e., within the interval [0.952, 1.05]), or the realizations for which the measurement data $\hat{h}_{-2}$ is close to 1 at 5%.

b) Study of AVAR for the higher $\tau$ value: It is well known that an AVAR measurement data for the higher $\tau$ value is $\chi_1^2$ distributed. Because each sequence contains 1024 samples, the higher $\tau$ argument of AVAR is $\tau_{max} = 1024/2 = 512$ arbitrary units (a.u.). From Table I, the theoretical AVAR result (model parameter) for $h_{-2} = 1$ a.u. and $\tau = 512$ a.u. is

$$ \sigma_y^2(\tau = 512) = \frac{2\pi^2 h_{-2}\tau}{3} = \frac{1024\pi^2}{3} = K_n \approx 3368.8. \quad (51) $$

Let us now consider the theoretical $\sigma_y^2(512)$ as the model parameter and the $\hat{\sigma}_y^2(512)$ computed value as the mean.

<footnote n="2">²This distribution represents, at best, the total ignorance of a value, i.e., a uniform law for the order of magnitude of this value. It was formalized as the Newcomb-Benford law [11], or first-digit law, and in particular explains the reason for the uneven wear of logarithm tables [12]. It is the basis of the concept of prior ignorance in Bayesian statistics [13].</footnote>
<figure><img src="image_7.png" /><figcaption>Fig. 6. Comparison of the histogram of the simulations and the theoretical conditional probability distribution function (PDF) (a) in the model world, for which the model parameter is $\theta_0 = 1$, and (b) in the measurement world, for which the measurement data are $\xi_0 = 1$.</figcaption></figure>

measurement data. To more easily compare these values, we normalize them by the coefficient $K_n$ given by (51):

$$ \left\{ \begin{aligned} \theta &= \frac{\sigma_y^2(512)}{K_n} = \frac{1024\pi^2 h_{-2}}{3K_n} \\ \xi &= \frac{\hat{\sigma}_y^2(512)}{K_n} = \frac{(\bar{Y}_2 - \bar{Y}_1)^2}{2K_n}, \end{aligned} \right. \qquad (52) $$

where $\bar{Y}_1$ and $\bar{Y}_2$ are, respectively, the average of the first 512 frequency deviation samples and the average of the last 512 frequency deviation samples of a considered simulated sequence.

Thus, we can select either the sequences for which $\theta = \theta_0 = 1 \pm 5\%$ (model world) or the sequences for which $\xi = \xi_0 = 1 \pm 5\%$ (measurement world). We can then easily compare the histogram of the $\xi$ measurement data, knowing that $\theta_0 = 1 \pm 5\%$, to the theoretical conditional probability given by (18) as well as the histogram of the $\theta$ parameter, knowing that $\xi_0 = 1 \pm 5\%$, to the theoretical conditional probability given by (22). Fig. 6 shows the

<figure><img src="image_8.png" /><figcaption>Fig. 7. Example of one realization among one million simulations of an Allan deviation (ADEV) curve for a pure random-walk FM with a unit noise level (obtained using SigmaTheta-1.2 [15]). The crosses represent the ADEV estimates, whereas the circles represent the square root of the log-unbiased estimates.</figcaption></figure>

very good agreement between histograms and theoretical conditional probability for both the model world and measurement world.

For the model world, we obtained 3568 sequences for which $\theta_0 = 1 \pm 5\%$ a.u., $\langle\xi\rangle = 0.977$ a.u., $\exp(\langle\ln(\xi)\rangle) = 0.276$ a.u., $\langle\exp(\hat{\mathcal{L}}_\theta)\rangle = 3.48$ a.u., $\exp(\langle\hat{\mathcal{L}}_\theta\rangle) = 0.983$ a.u. For the measurement world, we obtained 3523 sequences for which $\xi_0 = 1 \pm 5\%$ a.u., $\langle\theta\rangle = 19800$ a.u., $\exp(\langle\ln(\theta)\rangle) = 3.84$ a.u., and 3479 sequences for which $\exp(\hat{\mathcal{L}}_\theta) = 1 \pm 5\%$ a.u., $\langle\theta\rangle = 5630$ a.u., $\exp(\langle\ln(\theta)\rangle) = 1.04$ a.u.

Again, the results of this simulation are in perfect agreement with the theoretical considerations we have developed previously.

c) Results in the model world: The model parameter is now the random-walk FM level $h_{-2}$ and the measurement data are its estimation $\hat{h}_{-2}$. Among the million simulated sequences, 3568 have a model parameter $h_{-2}$ within the interval [0.952, 1.05]: model parameter, $\langle h_{-2} \rangle = 1.0007$ a.u.; fit over classical estimates, $\langle \hat{h}_{-2} \rangle = 0.911$ a.u. ($-8.9\%$); and fit over log-unbiased estimates, $\langle \hat{h}_{-2} \rangle = 0.911$ a.u. ($-0.9\%$).

The estimation of the random-walk noise level is clearly better with the fit over the log-unbiased estimates, but the estimation of $h_{-2}$ over the classical estimates is only slightly underestimated because it is based on all estimates, including those obtained for low $\tau$ values that are little biased (see Fig. 7).

2) White FM + Random Walk FM: We simulated 10000 sequences of a combination of white FM and random-walk FM in such a way that the random-walk FM is the dominating noise over the last decade. For the simulations, we used model parameter: $h_0 = 1$ a.u., $h_{-2} = 1.5 \cdot 10^{-5}$ a.u.; fit over classical estimates: $\langle \hat{h}_0 \rangle = 0.969$ a.u. ($-3.1\%$), $\langle \hat{h}_{-2} \rangle = 1.01 \cdot 10^{-5}$ a.u. ($-33\%$); and fit over log-unbiased esti-
<figure><img src="image_9.png" /><figcaption>Fig. 8. Example of one realization among 10000 simulations of an Allan deviation (ADEV) curve for a combination of white FM and random-walk FM (obtained using SigmaTheta-1.2 [15]). As in Fig. 1, the graph shows ADEV but the fit was performed over the Allan variance (AVAR) estimates. The crosses represent the ADEV estimates, whereas the circles represent the square root of the log-unbiased estimates. For comparison, Theol estimates [14] are represented as asterisks.</figcaption></figure>

mates: $\langle\hat{h}_0\rangle = 0.983$ a.u. ($-1.7\%$), $\langle\hat{h}_{-2}\rangle = 1.42 \cdot 10^{-5}$ a.u. ($-5.5\%$).

In this case, the estimation of the random-walk FM noise level is drastically improved by the use of the log-unbiased estimates because the random walk asymptote is fitted over the very last $\tau$ values, which are significantly affected by the log-log bias (see Fig. 8). For comparison, Theol estimates [14] were added on Fig. 8.

## IV. CONCLUSION

AVAR estimates are not biased and are easily fitted by weighted relative least squares. However, $\ln(\text{AVAR})$ estimates are biased, and this log-log bias affects the fits performed on log-log plots.

A log-unbiased estimator has been defined in such a way that its log is an unbiased estimator of the log of the model parameter. On a log-log AVAR plot, it is then recommended to fit the log-unbiased estimates rather than the classical AVAR estimates.

For this purpose, we developed the SigmaTheta software package [15] which performs fits according to the criterion (4) including log-unbiased estimator. It is free software developed under a GNU license for UNIX and Windows (Microsoft Corp., Redmond, WA) systems, available at http://theta obs-besancon.fr/spip.php?article103.

## REFERENCES

[1] J. Rutman, “Characterization of phase and frequency instabilities in precision frequency sources: Fifteen years of progress,” *Proc. IEEE*, vol. 66, no. 9, pp. 1048–1075, Sep. 1978.

[2] F. Vernotte, E. Lantz, J. Groslambert, and J. Gagnepain, “Oscillator noise analysis: Multivariate measurement,” *IEEE Trans. Instrum. Meas.*, vol. 42, no. 2, pp. 342–350, Apr. 1993.

[3] D. W. Allan, “Statistics of atomic frequency standards,” *Proc. IEEE*, vol. 54, no. 2, pp. 221–230, Feb. 1966.

[4] J. A. Barnes, A. R. Chi, L. S. Cutler, D. J. Healey, D. B. Lesson, T. E. McCunigal, J. A. Mullen, W. L. Smith, R. L. Sydnor, R. Vessot, and G. M. R. Winkler, “Characterization of frequency stability,” *IEEE Trans. Instrum. Meas.*, vol. IM-20, pp. 105–120, May 1971.

[5] P. Lesage and C. Audoin, “Characterization of frequency stability: Uncertainty due to the finite number of measurements,” *IEEE Trans. Instrum. Meas.*, vol. IM-22, pp. 157–161, Jun. 1973. See also corrections published in vol. IM-23, p. 103, Mar. 1974; vol. IM-24, p. 86, Mar. 1975; and vol. IM-25, p. 270, Sep. 1976.

[6] C. Greenhall and W. Riley, “Uncertainty of stability variances based on finite differences,” in *35th Annu. Precise Time and Time Interval Meeting*, 2003, pp. 267–280.

[7] A. Tarantola, *Inverse Problem Theory*. Amsterdam, The Netherlands: Elsevier, 1987.

[8] M. Evans, N. Hastings, and B. Peacock, *Statistical Distributions*, 4th ed., (Wiley Series in Probability and Statistics). New York, NY: Wiley, 2010.

[9] D. V. Lindley, “Fiducial distributions and Bayes theorem,” *J. R. Stat. Soc. B*, vol. 20, no. 1, pp. 102–107, 1958.

[10] M. Abramowitz and I. Stegun, *Handbook of Mathematical Functions*. New York, NY: Dover, 1964.

[11] F. Benford, “The law of anomalous numbers,” *Proc. Am. Philos. Soc.*, vol. 78, no. 4, pp. 551–572, Mar. 1938.

[12] S. Newcomb, “Note on the frequency of use of the different digits in natural numbers,” *Am. J. Math.*, vol. 4, no. 1, pp. 39–40, 1881.

[13] J. Bernardo and A. Smith, *Bayesian Theory*. (Wiley Series in Probability and Statistics). New York, NY: Wiley, 1994.

[14] D. A. Howe and T. N. Tasset, “dTheol: Characterization of very long-term frequency stability,” in *Proc. 18th European Frequency and Time Forum*, 2004, pp. 581–587.

[15] F. Vernotte, P. Y. Bourgeois, and F. Meyer, (2011, Apr. 26) “SigmaTheta software package,” [Online]. Available: http://theta obs-besancon.fr/spip.php?article103, GNU license.