# Testing statistical palaeomagnetic field models against directional data affected by measurement errors

A. Khokhlov, G. Hulot, C. Bouligand

► To cite this version:

A. Khokhlov, G. Hulot, C. Bouligand. Testing statistical palaeomagnetic field models against directional data affected by measurement errors. Geophysical Journal International, 2006, 167, pp.635-648. 10.1111/j.1365-246X.2006.03133.x. insu-03603137

HAL Id: insu-03603137

https://insu.hal.science/insu-03603137

Submitted on 9 Mar 2022

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
# Testing statistical palaeomagnetic field models against directional data affected by measurement errors

A. Khokhlov,<fnref n="1" />,<fnref n="2" /> G. Hulot<fnref n="2" /> and C. Bouligand<fnref n="2" />

Accepted 2006 July 7. Received 2006 June 22; in original form 2005 December 30

## SUMMARY

In a previous paper, Khokhlov et al. introduced a method to test the compatibility of so-called ‘giant Gaussian process’ (GGP) statistical models of the palaeomagnetic field against any palaeosecular variation database. This method did not take measurement errors into account. It therefore lacked practical usefulness. In the present paper, we remedy this and generalize the method to account for measurement errors in a way consistent with both the assumptions underlying the GGP approach and the nature of those errors. The method is implemented to test GGP models against any directional data set affected by Fisherian errors. Simulations show that the method can usefully discriminate which GGP model best explains a given data set. Applying the method to test six published GGP models against a test Bruhnes stable polarity data set extracted from the Quidelleur et al. database, it is found that all but one model (that of Quidelleur & Courtillot) should be rejected. Although this result should be taken with care, and does not necessarily imply that this model is superior to other models (Quidelleur & Courtillot precisely used the Quidelleur et al. database to infer their model), it clearly shows that in practice also, and with the databases currently available, the method can discriminate various candidate GGP models. It also shows that the statistical behaviour of the geomagnetic field at times of stable polarity can indeed be described in a consistent way in terms of a GGP model. This ‘forward’ testing method could ultimately be used to design an ‘inverse’ approach to GGP modelling of the palaeomagnetic field.

**Key words:** archaeomagnetism, geomagnetism, palaeomagnetism, potential field, spherical harmonics, statistical methods.

## INTRODUCTION

Thanks to the large amount of magnetic data provided by satellite missions spherical harmonic models of the main magnetic field produced within the core can readily be computed, providing high-resolution pictures of the way this field has been behaving over the past few decades (Hulot et al. 2002). Additional information can also be recovered from observations carried out by many generations of observers, explorers and navigators worldwide to produce spherical harmonic models describing the main magnetic field over the past four centuries (Jackson et al. 2000). Reconstructing similar, albeit much less accurate, spherical harmonic models of the Earth’s main magnetic field further back in time is also possible (Hongre et al. 1998; Korte & Constable 2005). This, however, requires that indirect measurements be used, only available through human artefacts, lava flows and sediments that have been magnetized in the ancient field. It also requires that a good age control of each sample is available to ensure a satisfying synchronization of the data used in computing spherical harmonic models for a given epoch. Because changes in the non-dipole component of the main field occur on

short timescales for small spatial scales (Hulot & Le Mouël 1994), a poor control indeed limits our ability to produce such spherical harmonic models back in time (so far, only up to 7000 yr in the past, Korte & Constable 2005).

To recover information about the Earth’s main magnetic field further back in time, one must rely on palaeomagnetic data and acknowledge the impossibility of accurately (within a few decades at most) synchronizing data acquired at different locations. This is especially true for the so-called palaeosecular variation (Quidelleur et al. 1994; Johnson & Constable 1996; McElhinny & McFadden 1997) and palaeointensity (Tanaka et al. 1995; Perrin & Schnepp 2004) databases covering the past few million years. Those databases encompass data recovered from volcanic samples and testifying for the instantaneous value of the direction and/or intensity of the field at well-known locations but poorly known times (typically within much more than a millennium). Fortunately, ages can nevertheless be measured with enough accuracy to identify the chron during which each of these samples acquired its magnetization. Assuming that the geodynamo essentially remained in a stationary state at times of stable polarity, this then opens the possibility, recognized

<footnote n="1">¹International Institute of Earthquake Prediction Theory and Mathematical Geophysics 79, b2, Warshavskoe shosse 117556 Moscow, Russia</footnote>

<footnote n="2">²Equipe de Géomagnétisme, Institut de Physique du Globe de Paris (Institut de recherche associé au CNRS et à l'Université Paris 7), 4, Place Jussieu, 75252, Paris, France. E-mail: gh@ipgp.jussieu.fr</footnote>
long ago (see e.g. Merrill *et al.* 1996, for an extensive review of earlier work), that at least some statistical properties could be recovered, characterizing the ancient field at times of stable polarity over the past millions years.

Several approaches have been used to try and infer such statistical properties (again, see e.g. Merrill *et al.* 1996). One approach has drawn particular attention, the so-called giant Gaussian process (GGP) approach first introduced by Constable & Parker (1988), and next generalized by Hulot & Le Mouël (1994) and Kono & Tanaka (1995). In its most general form (see e.g. Hulot & Bouligand 2005), this approach consists in

(1) using the classical spherical harmonic decomposition of the Main Field:

$$ \mathbf{B}(\Re, \Theta, \Psi, t) = -\nabla \left[ a \sum_{l=1}^{\infty} \sum_{m=0}^{l} \left( \frac{a}{\Re} \right)^{l+1} (g_l^m(t) \cos m\Psi + h_l^m(t) \sin m\Psi) P_l^m(\cos \Theta) \right], \quad (1) $$

where $g_l^m(t)$ and $h_l^m(t)$ are the time-varying Gauss coefficients, $a$ denotes the Earth's radius and $\{\Re, \Theta, \Psi\}$ are the standard spherical coordinates (i.e. distance from the Earth's centre, colatitude and longitude), and

(2) assuming that the $g_l^m(t)$ and $h_l^m(t)$ are the coordinates of a vector $\mathbf{k}(t)$ in a multidimensional model space, which behaves as a single realization of a multidimensional stationary random Gaussian process characterized by a statistical mean (or mean model) $E(\mathbf{k}(t))$ and a covariance matrix $\text{Cov}(\mathbf{k}(t), \mathbf{k}(t'))$ (which may very well be non-diagonal).

In addition to being defined with the help of a fully consistent formalism, this GGP approach has the unique advantage of providing a common tool to analyse just as well the historical (Constable & Parker 1988; Hulot & Le Mouël 1994), archaeomagnetic (Hongre et al. 1998) and palaeomagnetic fields (Constable & Parker 1988 and many studies since, see e.g. Kono et al. 2000a; Tauxe & Kent 2004). Even more, it can also be used to carry on similar analysis of the field produced by dynamo numerical simulations (Kono et al. 2000b; McMillan et al. 2001; Bouligand et al. 2005). In the latter case, considerable advantage can be taken of the fact that the Gauss coefficients (and very long time-series) are readily available and can thus be analysed in much detail without having to worry about any observational errors. Those studies have shown that under reasonable circumstances, the field produced by a numerical dynamo is indeed compatible with a GGP description, which can then be used to detect very interesting statistical properties, most notably, spontaneous and forced (because of the boundary conditions imposed at the core–mantle boundary) symmetry breaking properties (Bouligand et al. 2005; Hulot & Bouligand 2005). The GGP approach thus appears to potentially be a tool of choice to characterize the Earth’s main magnetic field at various epochs (present and past), detect possible changes in the regime of the geodynamo, some of which could be due to changes in the boundary conditions, and decide whether a numerical dynamo simulation is ‘Earth-like’ or not.

For this to be possible, we, however, need to have robust and sensitive numerical tools to actually decide which GGP model best fits a given data set. The issue is unfortunately not so trivial when the data to be analysed are the palaeomagnetic data from the palaeosecular variation databases (the largest available palaeomagnetic database for that type of study). In that case indeed, because the (directional) data are non-linearly related to the Gauss coefficients of the main field, no simple method is available to directly produce estimates

of the mean and covariance matrix of the GGP process best describing the data. Furthermore, using approximate approaches can easily lead to biased estimates (Khokhlov *et al.* 2001; Hatakeyama & Kono 2001). The purpose of the present paper, which expands on the earlier work of Khokhlov *et al.* (2001) (hereafter paper I), is to show that in fact a rigorous method can be implemented which (1) is fully consistent with all the assumptions involved in the GGP approach, (2) takes into account the measurement errors provided with each data point (a major improvement with respect to paper I), (3) provides a rigorous quantified assessment of the compatibility of any GGP model with any directional data set and (4) proves very discriminating.

DEALING WITH ERROR-FREE
PALAEDIRECTIONAL DATA

In this section, we first recall the approach proposed in paper I to test any error-free palaeodirectional data set. [Note, however, that we will make use of slightly different notations to avoid confusion among the various quantities we need to define. As a general rule, in the present paper, we will indeed always make use of gothic letters, namely g, s, ℓ and p (possibly with indexes) for the various probability density functions (PDFs) we will need to introduce.]

When dealing with such data, which consist in instantaneous values of the palaeodirection recorded at times separated by much more than the memory of the GGP process, each palaeomagnetic datum can be viewed as a local (both in time and space) independent realization of the GGP process. As a result, the covariance matrix of the process can be assumed to be of the simplified form $\text{Cov}(\mathbf{k}(t), \mathbf{k}(t')) = \text{Cov}(\mathbf{k}, \mathbf{k})\delta(t - t')$ and temporal correlations ignored altogether. Rewriting eq. (1) in the simplified form $\mathbf{B}(\Re, \Theta, \Psi, t) = A(\Re, \Theta, \Psi)\mathbf{k}(t)$, where $A(\Re, \Theta, \Psi)$ is a matrix projecting the multidimensional vector $\mathbf{k}(t)$ into the 3-D vector $\mathbf{B}(\Re, \Theta, \Psi, t)$, a given instantaneous value $\mathbf{B}$ of the field at a given location $(\Re, \Theta, \Psi)$ at the Earth's surface can then itself be viewed as a random drawing from a 3-D Gaussian distribution of random vectors $\mathbf{V}$ characterized by $E(\mathbf{V}) = A(\Re, \Theta, \Psi)E(\mathbf{k})$ and $\text{Cov}(\mathbf{V}, \mathbf{V}) = A(\Re, \Theta, \Psi)\text{Cov}(\mathbf{k}, \mathbf{k})A(\Re, \Theta, \Psi)^T$. Denote g, the corresponding PDF in $\mathbb{R}^3$.

Next, introduce the spherical coordinates $(\mathbf{u}, \rho)$ of the vector $\mathbf{V}$ (where $\mathbf{u} = \mathbf{V}/|\mathbf{V}| \in \Sigma$, the unit sphere in $\mathbb{R}^3$, and $\rho = |\mathbf{V}|$ denotes the distance from the origin $\mathbf{O}$ of the sphere $\Sigma$). Denote s the PDF associated with the direction $\mathbf{u}$ on $\Sigma$. This so-called angular Gaussian distribution (Bingham 1983, in paper I, we also referred to it as the Gaussian Directional distribution) is then defined by:

$$ s(\mathbf{u}) = \int_{0}^{\infty} g(\rho \mathbf{u}) d\rho. \quad (2) $$

In a local Cartesian coordinate system, we may write $E(\mathbf{V}) = \mathbf{m} = (m_1, m_2, m_3)$ and $\text{Cov}(\mathbf{V}, \mathbf{V}) = [\text{cov}(V_i, V_j)]$. Then, let $\mathbf{\Lambda} = [\Lambda^{ij}]$ be the inverse (hence, also symmetric) matrix of $\text{Cov}(\mathbf{V}, \mathbf{V})$. With respect to these local Cartesian coordinates, the PDF of $\mathbf{V}$ is:

$$ g(x_1, x_2, x_3) = \sqrt{\frac{\det \mathbf{\Lambda}}{(2\pi)^3}} \exp \left[ -\frac{1}{2} \sum_{i,j=1}^{3} \Lambda^{ij} (x_i - m_i)(x_j - m_j) \right], \quad (3) $$

or

$$ g(\mathbf{x}) = \sqrt{\frac{\det \mathbf{\Lambda}}{(2\pi)^3}} \exp \left[ -\frac{1}{2} (\mathbf{\Lambda}(\mathbf{x} - \mathbf{m}), \mathbf{x} - \mathbf{m}) \right]. \quad (4) $$

Making use of the $\mathbf{\Lambda}$-inner product $(\mathbf{x}, \mathbf{y})_{\mathbf{\Lambda}} = (\mathbf{\Lambda}\mathbf{x}, \mathbf{y}) = \sum_{i,j=1}^{3} \Lambda^{ij} x_i y_j$ together with the corresponding $\mathbf{\Lambda}$-norm $|\mathbf{x}|_{\mathbf{\Lambda}} =$
$\sqrt{(\mathbf{x}, \mathbf{x})}_{\Lambda}$, and turning back to the spherical coordinates $(\mathbf{u}, \rho)$ in $\mathbb{R}^3$, we may then write:

$$g(\rho\mathbf{u}) = \sqrt{\frac{\det \mathbf{\Lambda}}{(2\pi)^3}} e^{-\frac{1}{2}|\rho\mathbf{u}-\mathbf{m}|_{\Lambda}^2} \rho^2, \qquad (5)$$

and, after $\rho$-integration (eq. 2), finally get

$$\mathbf{s}(\mathbf{u}) = e^{-\frac{1}{2}m^2} \cdot \frac{\sqrt{\det \mathbf{\Lambda}}}{4\pi|\mathbf{u}|_{\Lambda}^3} \left[ z\sqrt{\frac{2}{\pi}} + e^{\frac{1}{2}z^2} (1+z^2) \left[ 1 + \mathrm{Erf}\left(\frac{z}{\sqrt{2}}\right) \right] \right], \qquad (6)$$

where

$$z = \frac{(\mathbf{m}, \mathbf{u})_{\Lambda}}{|\mathbf{u}|_{\Lambda}}, \qquad m = |\mathbf{m}|_{\Lambda}, \qquad (7)$$

correspond to respectively the $\Lambda$-projection of $\mathbf{m}$ on the direction $\mathbf{u}$, and the $\Lambda$-norm of $\mathbf{m}$. Since $\Lambda$ is positive $-m \le z \le m$.

Eq. (6) makes it possible to predict the PDF of the direction of the field at any location at the Earth's surface. Any given GGP model can then easily be tested against the data set corresponding to such a location.

Unfortunately, relatively few data are usually available at a given location. In paper I, we explained how this drawback could be overcome, and sparse data distributed over various sites combined into a uniformized data set to globally test the GGP model.

This consists in converting the directional data set $\{\mathbf{u}_i\}$ under consideration (where the 'i' index varies between 1 and N, the total amount of data) into a univariate data set $\{t_i\}$ with the help of:

$$t_i = P\{\mathbf{u}|\mathbf{s}_i(\mathbf{u}) \ge \mathbf{s}_i(\mathbf{u}_i)\} = \int_{\{\mathbf{u}|\mathbf{s}_i(\mathbf{u}) \ge \mathbf{s}_i(\mathbf{u}_i)\}} \mathbf{s}_i(\mathbf{u}) dU, \qquad (8)$$

where $dU$ is the elementary surface about $\mathbf{u}$ on $\Sigma$, and testing the distribution of the $\{t_i\}$ against a uniform probability density over the segment $[0, 1]$.

This uniformization procedure simply amounts, for each data $\mathbf{u}_i$, to identify on $\Sigma$ the isoprobability line of the PDF $\mathbf{s}_i(\mathbf{u})$ on which the data $\mathbf{u}_i$ lie (note that this PDF depends on the location the data $\mathbf{u}_i$ come from, which is why an 'i' index is being added to $\mathbf{s}_i(\mathbf{u})$), and to assign to $t_i$ the value of the probability that $\mathbf{u}_i$ could have been lying on a higher isoprobability line on $\Sigma$. Then, if the GGP model is compatible with the data set (hence, if the $\mathbf{u}_i$ are compatible with the $\mathbf{s}_i(\mathbf{u})$), the $\{t_i\}$ should be compatible with a uniform distribution on $[0, 1]$.

Whether this is indeed the case can finally be assessed with the help of various standard tests, (in paper I we used the Kolmogorov-Smirnov (KS) and the $\chi^2$ tests). Note that this uniformization procedure is very general and does not rely on any properties of the angular Gaussian distribution. It can be applied in more general situations.

## TAKING MEASUREMENT ERROR INTO ACCOUNT

A main issue that paper I did not address is that of data measurement errors. The tests we just described are appropriate only in the event these errors are negligible. Unfortunately this is far from being the case and we now need to take this into account. In the case of directional data, these measurement errors are classically described in terms of a Fisher distribution characterized by a so-called concentration (or precision) parameter $K$. Within this Fisherian framework, stating that the observed direction $\mathbf{s}$ reflects the actual direction $\mathbf{u}$ on $\Sigma$ with an error characterized by a concentration parameter $K$,

amounts to state that $\mathbf{s}$ is the result of a random drawing from a Fisher distribution centred on $\mathbf{u}$ with PDF:

$$\mathbf{k}_K(\mathbf{u}, \mathbf{s}) = \frac{K}{2\pi(e^K - e^{-K})} e^{K \cdot (\mathbf{u}, \mathbf{s})}. \qquad (9)$$

In spherical coordinates and when centred on $\mathbf{u}_0 = (\theta = 0, \varphi = 0)$, this takes the more usual form:

$$\mathbf{k}_K(\mathbf{u}_0, \mathbf{s}) = \frac{K e^{K \cos \theta}}{2\pi(e^K - e^{-K})}. \qquad (10)$$

This means that in order to test a given GGP model against a given data set $\{\mathbf{u}_i\}$ with associated errors characterized by $\{K_i\}$, we now need to test whether each $\mathbf{u}_i$ can be considered the result of a realization $\mathbf{s}$ of the GGP process (distributed over $\Sigma$ with the PDF $\mathbf{s}_i$), shifted to $\mathbf{u}_i$ as a result of an error with a Fisher distribution $\mathbf{k}_i(\mathbf{s}, \mathbf{u})$ centred on $\mathbf{s}$ and a concentration parameter $K_i$.

The infinitesimal probability that the GGP process first produces a direction within the elementary surface $dS$ about $\mathbf{s}$, and that the error next independently shifts this direction to within the elementary surface $dU$ about $\mathbf{u}_i$ is:

$$\mathbf{s}_i(\mathbf{s}) \, dS \cdot \mathbf{k}_{K_i}(\mathbf{s}, \mathbf{u}) \, dU. \qquad (11)$$

Then, integrating over all possible intermediate $\mathbf{s}$ directions leads to the probability $p_i(\mathbf{u})dU$ of finding $\mathbf{u}_i$ within $dU$ about $\mathbf{u}$, where:

$$p_i(\mathbf{u}) = \int_{\Sigma} \mathbf{s}_i(\mathbf{s}) \mathbf{k}_{K_i}(\mathbf{s}, \mathbf{u}) \, dS, \qquad (12)$$

is the new PDF against which $\mathbf{u}_i$ should be tested. Thus, taking data errors into account only amounts to use $p_i(\mathbf{u})$ in place of $\mathbf{s}_i(\mathbf{u})$ in eq. (8). Tests can then again be used, either on a site-by-site basis or, more interestingly, on a regional and global scale after using the uniformization procedure, which now becomes:

$$t_i = P\{\mathbf{u}|p_i(\mathbf{u}) \ge p_i(\mathbf{u}_i)\} = \int_{\{\mathbf{u}|p_i(\mathbf{u}) \ge p_i(\mathbf{u}_i)\}} p_i(\mathbf{u}) dU. \qquad (13)$$

## PRACTICAL IMPLEMENTATION

In practice, testing a given GGP model against a given directional data set thus involves four successive steps:

(1) for each data $\mathbf{u}_i$, to compute the error-free PDF $\mathbf{s}_i(\mathbf{u})$ predicted by the GGP model at the site where $\mathbf{u}_i$ was collected (the analytical form of which is given by eq. 6);

(2) to compute the error-included PDF $p_i(\mathbf{u})$ through the convolution (12);

(3) to produce the uniformized data $t_i$ with the help of eq. (13) and

(4) to test the uniformized data set $\{t_i\}$ against a uniform distribution over $[0, 1]$.

Unfortunately, no exact analytical solutions of eqs (12) and (13) are known to us. In principle, this is not too much of a problem, since both eqs (12) and (13) can be computed numerically. In practice however, the numerical implementation of formula (12) requires quite some computational time. This drawback can be considered negligible if we simply test a single GGP model against a small
database, or if all data **u**<sub>i</sub> coming from the same site share the same data error (in which case **p**<sub>i</sub>(**u**) only needs to be computed once for each of the typically thirty something sites). However, if **p**<sub>i</sub>(**u**) needs to be computed for each data of the data set (i.e. typically a thousand times), this computational time can become very long (on the order of a few days on a current PC). For that kind of more realistic situations (and for the 'forward' testing method we propose here to possibly become of any use for future much more computer-intensive 'inverse' searches of best models), faster algorithms are clearly desirable. This prompted us to further look for an approximate, accurate enough, analytic solution to the convolution (12).

As is described in the appendix, one such approximation can indeed be found. For any given data **u**<sub>i</sub>, this consists in approximating **p**<sub>i</sub>(**u**) by the angular Gaussian distribution defined by eqs (6) and (7), where **m** is then the mean field value predicted by the GGP model to be tested at the site associated with the data **u**<sub>i</sub> (i.e. predicted by **g**<sub>i</sub>(**x**) as defined by eq. (4) for the site associated with the data **u**<sub>i</sub>), and **Λ** is (Cov(**V**, **V**) + |**m**|<sup>2</sup> |**I**|<sup>-1</sup>, where Cov(**V**, **V**) is the covariance matrix associated with **g**<sub>i</sub>(**x**) and K (=**K**<sub>i</sub>) is the concentration parameter characterizing the Fisherian error associated with **u**<sub>i</sub>.

In what follows, whenever we use this approximation of **p**<sub>i</sub>(**u**) to implement the tests, we will refer to the 'approximate' method, and whenever a direct numerical implementation of the exact convolution (12) is used, we will refer to the 'exact' method. It turns out that, as we shall see, the approximate method is most often just as good as the exact method, and about 100 times faster to run.

## MODELS, DATA SET, AND STATISTICAL TOOLS USED IN THIS STUDY

To illustrate the power of the method we propose, several published GGP models have been tested. We will refer to these as:

(i) CP model, which is the preferred model of Constable & Parker (1988);

(ii) QC model, C1 (preferred) model of Quidelleur & Courtillot (1996);

(iii) CJ model, which is the CJ98 model proposed by Constable & Johnson (1999);

(iv) JC model, which is the CJ98.nz model also proposed by Constable & Johnson (1999);

(v) TK model, which is the TK03.GAD model recently proposed by Tauxe & Kent (2004) and further discussed in Tauxe (2005);

(vi) HK model, which is the final model for the normal polarity of Hatakeyama & Kono (2002).

Models CP, QC, CJ, JC and TK share many characteristics. They are defined by simple axisymmetric mean models for which only E(g<sub>1</sub><sup>0</sup>) and E(g<sub>2</sub><sup>0</sup>) can take non-zero values, and purely diagonal covariance matrices with Cov(g<sub>n</sub><sup>m</sup>, g<sub>n</sub><sup>m</sup>) = Cov(h<sub>n</sub><sup>m</sup>, h<sub>n</sub><sup>m</sup>) = (σ<sub>n</sub><sup>m</sup>)<sup>2</sup>, except in the case of model JC, which assumes different values for Cov(g<sub>1</sub><sup>2</sup>, g<sub>2</sub><sup>1</sup>) = (σ(g<sub>2</sub><sup>1</sup>))<sup>2</sup> and Cov(h<sub>1</sub><sup>2</sup>, h<sub>2</sub><sup>1</sup>) = (σ(h<sub>2</sub><sup>1</sup>))<sup>2</sup>. In all but one case, for n ≥ 3, the σ<sub>n</sub><sup>m</sup> are further assumed to be independent of m and defined by σ<sub>n</sub><sup>m</sup> = σ<sub>n</sub> where σ<sub>n</sub> = α(c/a)<sup>n</sup> / ((n + 1)(2n + 1))<sup>1/2</sup>. The only exception is model TK which distinguishes σ<sub>n</sub><sup>m</sup> = σ<sub>n</sub> for (n − m) even from σ<sub>n</sub><sup>m</sup> = βσ<sub>n</sub> for (n − m) odd. Table 1 gives the values of the relevant parameters for each of those five models. Finally, model HK differs from the other models because of a more elaborate mean field, defined up to degree and order 4. However, its covariance matrices are otherwise defined in much the same way. A

<table><caption>**Table 1.** Parameters defining models CP, QC, CJ, JC and TK. Units are in μT, except for β only defined for model TK and which is dimensionless (see main text for further details). Parameters defining model HK are to be found in Table 2 of Hatakeyama & Kono (2002).</caption><thead><tr><th></th><th>CP</th><th>QC</th><th>CJ</th><th>JC</th><th>TK</th></tr></thead><tbody><tr><td>E(g<sub>1</sub><sup>0</sup>)</td><td>-30.0</td><td>-30.0</td><td>-30.0</td><td>-30.0</td><td>-18.0</td></tr><tr><td>E(g<sub>2</sub><sup>0</sup>)</td><td>-1.8</td><td>-1.2</td><td>-1.5</td><td>-1.5</td><td>0.0</td></tr><tr><td>σ<sub>1</sub><sup>0</sup></td><td>3.0</td><td>3.0</td><td>11.72</td><td>11.72</td><td>6.36</td></tr><tr><td>σ<sub>1</sub><sup>1</sup></td><td>3.0</td><td>3.0</td><td>1.67</td><td>1.67</td><td>1.67</td></tr><tr><td>σ<sub>2</sub><sup>0</sup></td><td>2.14</td><td>1.3</td><td>1.16</td><td>1.16</td><td>0.58</td></tr><tr><td>σ(g<sub>1</sub><sup>2</sup>)</td><td>2.14</td><td>4.3</td><td>4.06</td><td>1.16</td><td>2.20</td></tr><tr><td>σ(h<sub>1</sub><sup>2</sup>)</td><td>2.14</td><td>4.3</td><td>4.06</td><td>8.12</td><td>2.20</td></tr><tr><td>σ<sub>2</sub><sup>2</sup></td><td>2.14</td><td>1.3</td><td>1.16</td><td>1.16</td><td>0.58</td></tr><tr><td>α</td><td>27.7</td><td>27.7</td><td>15.0</td><td>15.0</td><td>7.5</td></tr><tr><td>β</td><td></td><td></td><td></td><td></td><td>3.8</td></tr></tbody></table>

full description of this model can be found in table 2 of Hatakeyama & Kono (2002).

All those models have been constructed in the hope that they would properly describe the statistical properties of the palaeomagnetic field at times of normal polarity over the past 5 Myr. In particular, they have been constructed with the help of databases including a large number (if not a majority) of directional data corresponding to the current Bruhnes chron (e.g. Quidelleur et al. 1994; Johnson & Constable 1996). They are, therefore, good candidates for an example test against a well-controlled data set covering the Bruhnes chron and corresponding to volcanic directional data acquired at various sites distributed worldwide. Ideally, such a test data set would have to be built by extracting data from the most recent and permanently updated IAGA palaeomagnetic reference database for palaeofield direction available at the National Geophysical Data Center (http://www.ngdc.noaa.gov/seg/geomag/paleo.shtml), and by relying on a set of fine-tuned criteria agreeable to the community. This, however, is not a trivial matter since, as noted by one of the reviewers, there currently is no general agreement among investigators on what comprises a satisfactory database for this kind of study. Such a substantial endeavour, which we plan to carry on at a later stage, is therefore clearly beyond the scope of the present paper, which only intends to introduce, test, and illustrate a new methodology. For such a purpose, we felt that a simpler, more readily available data set (and one that at least already went through some type of selection processes relevant to the present study), would be appropriate enough.

We, therefore, decided to use a test data set extracted from the Quidelleur et al. (1994) database, originally used by Quidelleur & Courtillot (1996) to build their QC model, last updated in January 1998, and currently available at http://www.ipgp.jussieu.fr/rech/paleomag/var-secu/. We note however, even before proceeding further, that because of this choice, a close fit of the QC model to our test data set can not be anticipated, so that test results reported here will necessarily not test the QC model as stringently as the other models. This test data set consists of a total of 990 independent estimates of the local direction of the palaeomagnetic field at 36 sites (Fig. 1, Table 2). Each site is located with the help of its latitude and longitude. At each site, each estimate is based on the direction of the resultant vector **R** of n (≥3) volcanic samples (unit vectors) and is given in the form of a declination D and an inclination I together with both n and the norm **R** = |**R**| of the resultant vector. Although **R** is usually not published as such in the original papers, it can be accurately recomputed from the published material. This then makes it possible to present all the data in
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Locations of the 36 sites of the test data set used in the present study.</figcaption></figure>

a homogeneous and complete form in the database. Indeed, and as is well known (e.g. McElhinny & McFadden 2000), an estimate of the local direction of the palaeomagnetic field is completely characterized (with is Fisherian error) as soon as the four parameters *D*, *I*, *n* and *R* are given.

For each estimate, the corresponding vector $\mathbf{u}_i$ could then be located on the unit sphere $\Sigma$, and an estimate of the concentration parameter $K_i$ defining the Fisherian error computed [with the help of $K = R(\frac{n-1}{n-R})$; note that this is not the concentration parameter $k = \frac{n-1}{n-R}$ classically published with the data, which characterizes the dispersion of the n individual samples about their average direction, but the concentration parameter characterizing the error on $\mathbf{u}_i$]. Starting from $\mathbf{u}_i$ and $K_i$, and for each GGP model to be tested, a uniformized value $t_i$ could finally be computed with the help of eq. (13) to produce the final uniformized data set $\{t_i\}$.

To test this final data set against a uniform distribution over [0, 1], we relied on two different well-known tests: the Kolmogorov-Smirnov test (KS-test) and the Anderson-Darling test (AD-test). These two tests rely on the fact that, if a given uniformized data set $\{t_i\}$ (with $i = 1$ to $N$) is indeed compatible with a uniform distribution over [0, 1], its empirical cumulative distribution function (CDF) $F_N(x)$ should fluctuate within predictable limits about the theoretical CDF value $F(x) = x$. The null hypothesis should then be rejected if the empirical CDF $F_N(x)$ either departs too much from $x$, or remains to close to it (the latter revealing the presence of suspicious regularities within the $\{t_i\}$). The KS-test and AD-test differ in the measure chosen to assess how distant $F_N(x)$ is from $x$ over [0, 1]:

(i) the KS-test (as defined by e.g. Press et al. 1996) uses the maximum value $M_N$ of $|F_N(x) - x|$ over [0, 1] and is, therefore, most sensitive to departures of the $\{t_i\}$ from a uniform distribution towards the middle of the segment [0, 1];

(ii) the AD-test (as defined by e.g. Marsaglia & Marsaglia 2004) uses the integral quantity $I_N = N \int_0^1 \frac{F_N(x)-x}{x(1-x)} dx$. Because of the weight $[x(1-x)]^{-1}$, it is much more sensitive to the behaviour of $\{t_i\}$ at both extremes of the segment [0, 1].

The KS and AD tests are both excellent tests for small samples and can also be used for large samples. Neither require any additional parameter and both can, therefore, be considered as completely objective. Also, they complement each other very well. In practice, for each of the uniformized data set $\{t_i\}$ we had to test; we, therefore, computed the values of $M_N$ and $I_N$, together with (from known software, e.g. Press et al. 1996; Marsaglia & Marsaglia 2004) the probabilities $P(M_N)$ and $P(I_N)$ for the null hypothesis to have possibly produced such large, or even larger, values for respectively $M_N$ and $I_N$.

<table><caption>**Table 2.** Site number, latitude, longitude and number of directional estimates available for each of the 36 sites of the test data set used in the present study.</caption><thead><tr><th>Site</th><th>Lat (deg)</th><th>Lon (deg)</th><th>N</th></tr></thead><tbody><tr><td>1</td><td>19.50</td><td>205.00</td><td>112</td></tr><tr><td>2</td><td>39.00</td><td>29.00</td><td>16</td></tr><tr><td>3</td><td>-35.50</td><td>173.60</td><td>21</td></tr><tr><td>4</td><td>43.40</td><td>2.80</td><td>31</td></tr><tr><td>5</td><td>38.70</td><td>332.80</td><td>21</td></tr><tr><td>6</td><td>4.50</td><td>9.50</td><td>14</td></tr><tr><td>7</td><td>3.50</td><td>9.00</td><td>37</td></tr><tr><td>8</td><td>-38.80</td><td>77.50</td><td>14</td></tr><tr><td>9</td><td>-21.10</td><td>55.50</td><td>17</td></tr><tr><td>10</td><td>-46.50</td><td>52.20</td><td>40</td></tr><tr><td>11</td><td>-46.50</td><td>51.70</td><td>34</td></tr><tr><td>12</td><td>-46.90</td><td>37.80</td><td>16</td></tr><tr><td>13</td><td>19.00</td><td>-99.00</td><td>36</td></tr><tr><td>14</td><td>-27.10</td><td>250.80</td><td>53</td></tr><tr><td>15</td><td>35.10</td><td>139.00</td><td>8</td></tr><tr><td>16</td><td>46.00</td><td>3.00</td><td>38</td></tr><tr><td>17</td><td>50.00</td><td>7.00</td><td>44</td></tr><tr><td>18</td><td>50.00</td><td>7.00</td><td>44</td></tr><tr><td>19</td><td>19.00</td><td>261.00</td><td>39</td></tr><tr><td>20</td><td>37.70</td><td>241.00</td><td>33</td></tr><tr><td>21</td><td>35.00</td><td>139.00</td><td>35</td></tr><tr><td>22</td><td>35.00</td><td>139.00</td><td>23</td></tr><tr><td>23</td><td>-1.00</td><td>271.00</td><td>16</td></tr><tr><td>24</td><td>-38.00</td><td>77.00</td><td>25</td></tr><tr><td>25</td><td>-21.00</td><td>55.00</td><td>27</td></tr><tr><td>26</td><td>-36.90</td><td>174.80</td><td>17</td></tr><tr><td>27</td><td>71.10</td><td>8.20</td><td>10</td></tr><tr><td>28</td><td>38.20</td><td>140.50</td><td>11</td></tr><tr><td>29</td><td>38.50</td><td>14.90</td><td>33</td></tr><tr><td>30</td><td>41.70</td><td>238.50</td><td>19</td></tr><tr><td>31</td><td>57.20</td><td>189.70</td><td>8</td></tr><tr><td>32</td><td>60.00</td><td>194.00</td><td>13</td></tr><tr><td>33</td><td>19.30</td><td>-101.00</td><td>9</td></tr><tr><td>34</td><td>20.10</td><td>-155.50</td><td>10</td></tr><tr><td>35</td><td>38.20</td><td>15.30</td><td>45</td></tr><tr><td>36</td><td>28.80</td><td>342.30</td><td>21</td></tr></tbody></table>

## VALIDATION AND DISCRIMINATING POWER OF THE METHOD

Before actually implementing the testing method we propose, we ran a number of instructive validating tests. For that purpose, various artificial data sets have been produced in the following way. Starting from a given GGP model, we first produced random error-free directional data values. Exactly the same number of data values were produced for each of the 36 sites, as available in the (real) test data set (thus amounting to a total of 990 independent artificial values). This could easily be done by just randomly drawing 990
independent sets of Gauss coefficients ($g_i^m, h_i^m$) from a multidimensional Gaussian generator with mean values and covariances chosen as specified by the GGP model under consideration, and computing the predicted values at each site. A random directional error was next added to each of those 990 error-free directional data values. Because, as shown in the appendix, a Fisher distribution is extremely close to an angular Gaussian distribution (for realistic Fisherian errors such as those we deal with here), this was actually done by introducing a random Gaussian error (an easier procedure to control than that involved in producing random Fisherian errors). Using this procedure and starting from various specified GGP models, a number of artificial directional data sets affected by specified errors could be produced, which we used to validate our method.

We first validated the 'exact' method and checked that for artificial data produced from any GGP models and perturbed by some known error, this method would always conclude to the compatibility of the data set with the starting GGP model, provided the data are assumed to be affected by the correct amount of error. This first series of tests not only allowed us to validate the method. It also allowed us to check our software and the level of numerical accuracy needed to produce accurate results. It turns out that from a numerical point of view, the most sensitive step of the entire method is the numerical implementation of the uniformization procedure (i.e. the computation of eq. 13). This indeed requires the numerical computation of $p_i(u)$ on a grid over the unit sphere $\Sigma$. The most demanding computations, in terms of grid size, happen to be those corresponding to $p_i(u)$ predicting the largest dispersions of the $u_i$ on the unit sphere $\Sigma$. In that case, when the grid is not tight enough, numerical errors happen to mainly affect the computation of the $t_i$ lying close to 0 or 1. This in turn mainly affects the results of the AD-test which produces too large values of $I_N$ and, therefore, too small values of $P(I_N)$. As a result, the AD-test tends to more easily produce negative results than it should.

For the purpose of the present paper, we used a grid corresponding to a decent compromise, tight enough for the KS-test to always lead to an accurate result, but not too tight to make it possible to run all tests in a reasonable amount of time. This had the drawback that, for the reasons just outlined, the AD-test produced slightly erroneous results when testing the CJ and JC models. For those two models, we therefore concluded that only the KS-test could be considered reliable. For all other models (CP, QC, TK and HK), both the KS-test and the AD-test were otherwise found to be reliable. Those conclusions have been reached when testing artificial data affected by errors characterized by $\alpha_{95}$ of up to $12.5^\circ$. (Note also that to produce those artificial data sets and run all the tests reported in this study, all GGP models have been used only up to degree and order 7, since, as we shall later see, higher degrees contribute negligibly to the final results for all the GGP models we tested).

We next proceeded to compare the results of the 'approximate' method to those of the 'exact' method, using exactly the same artificial data sets in both cases. Again and as anticipated (see appendix), significant disagreements between the two methods were found to arise only when $p_i(u)$ predicts large dispersions of the $u_i$ on the unit sphere $\Sigma$. In practice, the approximate method was found to lead to accurate enough results for all the models we tested (testing artificial data with all levels of errors up to $\alpha_{95} = 12.5^\circ$), except when testing models CJ and JC against data affected by errors with $\alpha_{95}$ greater than $7.5^\circ$. Table 3 summarizes the domain of validity of each type of test as established from all those preliminary tests.

We then went on to test the ability of the method to properly discriminate the source GGP model and the level of error affecting

<table><caption>**Table 3.** Domain of validity of the tests used in the present study. The test can be considered reliable if the error affecting the data is less than the $\alpha_{95}$ value shown in the table (units in degrees).</caption><thead><tr><th></th><th>CP</th><th>QC</th><th>CJ</th><th>JC</th><th>TK</th><th>HK</th></tr></thead><tbody><tr><td>KS-test with exact method</td><td>12.5</td><td>12.5</td><td>12.5</td><td>12.5</td><td>12.5</td><td>12.5</td></tr><tr><td>KS-test with app. method</td><td>12.5</td><td>12.5</td><td>7.5</td><td>7.5</td><td>12.5</td><td>12.5</td></tr><tr><td>AD-test with exact method</td><td>12.5</td><td>12.5</td><td>N.A.</td><td>N.A.</td><td>12.5</td><td>12.5</td></tr><tr><td>AD-test with app. method</td><td>12.5</td><td>12.5</td><td>N.A.</td><td>N.A.</td><td>12.5</td><td>12.5</td></tr></tbody></table>

the data. Because, as we shall soon see, the real data turn out to be compatible only with model QC (from the point of view of the tests described in the present paper), and since the average $\alpha_{95}$ for the 990 data in the real database is of $4.7^\circ$, this test was carried out on an artificial data set produced from model QC and affected by errors of $\alpha_{95} = 4.7^\circ$. This artificial data set was used to test all GGP models, assuming various possible levels of errors ($0^\circ, 2.5^\circ, 5^\circ, 7.5^\circ, 10^\circ, 12.5^\circ$). An additional test was also run to check that the method would properly conclude that this data set is compatible with the starting model QC when assuming an error of $\alpha_{95} = 4.7^\circ$. Since both the exact and approximate methods provide reliable results for errors of that magnitude, both methods have been used with the KS-test and AD-test (except for the CJ and JC models for which only the KS test was used, for the reasons outlined above). Table 4 reports the results of those tests. (Note that Table 4 also provides a good illustration of the accuracy of the results of the approximate method, when compared to the results of the exact method).

This table shows that, as expected, the artificial data set is compatible with the starting model (QC model) provided the correct level of error affecting the data is being assumed ($\alpha_{95} = 4.7^\circ$). It also shows that by contrast, making a wrong *a priori* assumption with respect to the errors affecting the data can lead to a negative result. A minimum error of $\alpha_{95} = 2.5^\circ$ and a maximum error of $\alpha_{95} = 7.5^\circ$ is indeed required for the test with model QC to produce positive results (with no rejection at a 95 per cent level of confidence by either of the KS or AD tests). Finally, Table 4 shows that the method also appears to be able to discriminate the various GGP models. The data set produced by QC with $\alpha_{95} = 4.7^\circ$ cannot be attributed to a CP model, unless (but marginally so) one wrongly assumes $\alpha_{95} = 10^\circ$. It can no more be attributed to a CJ, JC or HK model, unless one again wrongly assumes $\alpha_{95} = 12.5^\circ$. And it can hardly be attributed to a TK model, whatever the value assumed for $\alpha_{95}$. Altogether those results thus show that the tests we introduced are sensitive to both the choice of the model to be tested and the level of error assumed to affect the data. In particular they show that if we have a reasonable knowledge of the error affecting the data, the method applied to the type of data set we deal with here (i.e. 990 data distributed over the 36 sites of Fig. 1) makes it possible to fully discriminate the models. Indeed Table 4 shows that if we know *a priori* that the error affecting the data is somewhere between $2.5^\circ$ and $7.5^\circ$, the method can detect that all of the CP, CJ, JC, HK and TK models are incompatible with the data set produced from QC model with $\alpha_{95} = 4.7^\circ$, at more than a 95 per cent level of confidence. Those encouraging results show that the method proposed in the present paper has enough discriminating power to decide which, if any, of the six GGP models can best explain the test data set we have assembled.

## TESTING THE REAL DATA

Having validated the method, and assessed its limits and discriminating power, we then went on to test the various GGP models against the real test data set.
<table><caption>**Table 4.** Testing all models against an artificial data set produced from model QC and affected by errors of $\alpha_{95} = 4.7^{\circ}$. Results are shown for all models (first column), assuming various possible levels of errors (second column). For each of the KS and AD tests, both the approximate and exact methods have been used. Results are reported in the form of the $M_N$ and $P(M_N)$ values for the KS test, and the $I_N$ and $P(I_N)$ values for the AD test. Non-zero values of the $P(M_N)$ and $P(I_N)$ have been highlighted in bold. The AD test for the CJ and JC models being unreliable, their results have not been reported. Note also that when no error is being assumed (i.e. for a 0° error), the approximate method is not relevant and the calculation is directly carried out with the help of eq. (8) in place of eq. (13).</caption><thead><tr><th rowspan="3">model</th><th rowspan="3">error</th><th colspan="4">KS</th><th colspan="4">AD</th></tr><tr><th colspan="2">approx.</th><th colspan="2">exact</th><th colspan="2">approx.</th><th colspan="2">exact</th></tr><tr><th>M<sub>N</sub></th><th>P(M<sub>N</sub>)</th><th>M<sub>N</sub></th><th>P(M<sub>N</sub>)</th><th>I<sub>N</sub></th><th>P(I<sub>N</sub>)</th><th>I<sub>N</sub></th><th>P(I<sub>N</sub>)</th></tr></thead><tbody><tr><td rowspan="10">CP</td><td>0.0</td><td></td><td></td><td>0.0855</td><td>0.00</td><td></td><td></td><td>22.3</td><td>0.00</td></tr><tr><td>2.5</td><td>0.0804</td><td>0.00</td><td>0.0807</td><td>0.00</td><td>19.5</td><td>0.00</td><td>19.7</td><td>0.00</td></tr><tr><td>4.7</td><td>0.0692</td><td>0.00</td><td>0.0705</td><td>0.00</td><td>13.7</td><td>0.00</td><td>14.3</td><td>0.00</td></tr><tr><td>5.0</td><td>0.0666</td><td>0.00</td><td>0.0687</td><td>0.00</td><td>12.8</td><td>0.00</td><td>13.5</td><td>0.00</td></tr><tr><td>7.5</td><td>0.0521</td><td><strong>0.01</strong></td><td>0.0539</td><td><strong>0.01</strong></td><td>5.88</td><td>0.00</td><td>6.90</td><td>0.00</td></tr><tr><td>10.0</td><td>0.0374</td><td><strong>0.12</strong></td><td>0.0371</td><td><strong>0.13</strong></td><td>2.88</td><td><strong>0.03</strong></td><td>3.73</td><td><strong>0.01</strong></td></tr><tr><td>12.5</td><td>0.0747</td><td>0.00</td><td>0.0728</td><td>0.00</td><td>6.73</td><td>0.00</td><td>6.76</td><td>0.00</td></tr><tr><td>0.0</td><td></td><td></td><td>0.0415</td><td><strong>0.06</strong></td><td></td><td></td><td>2.57</td><td><strong>0.05</strong></td></tr><tr><td>2.5</td><td>0.0360</td><td><strong>0.15</strong></td><td>0.0358</td><td><strong>0.15</strong></td><td>1.90</td><td><strong>0.10</strong></td><td>1.87</td><td><strong>0.11</strong></td></tr><tr><td>4.7</td><td>0.0243</td><td><strong>0.60</strong></td><td>0.0247</td><td><strong>0.58</strong></td><td>0.766</td><td><strong>0.51</strong></td><td>0.734</td><td><strong>0.53</strong></td></tr><tr><td rowspan="6">QC</td><td>5.0</td><td>0.0230</td><td><strong>0.67</strong></td><td>0.0221</td><td><strong>0.71</strong></td><td>0.636</td><td><strong>0.61</strong></td><td>0.616</td><td><strong>0.63</strong></td></tr><tr><td>7.5</td><td>0.0231</td><td><strong>0.66</strong></td><td>0.0274</td><td><strong>0.44</strong></td><td>0.590</td><td><strong>0.66</strong></td><td>0.869</td><td><strong>0.43</strong></td></tr><tr><td>10.0</td><td>0.0467</td><td><strong>0.03</strong></td><td>0.0547</td><td><strong>0.01</strong></td><td>12.4</td><td>0.00</td><td>13.8</td><td>0.00</td></tr><tr><td>12.5</td><td>0.0658</td><td>0.00</td><td>0.0779</td><td>0.00</td><td>25.2</td><td>0.00</td><td>27.5</td><td>0.00</td></tr><tr><td>0.0</td><td></td><td></td><td>0.169</td><td>0.00</td><td></td><td></td><td></td><td></td></tr><tr><td>2.5</td><td>0.163</td><td>0.00</td><td>0.161</td><td>0.00</td><td></td><td></td><td></td><td></td></tr><tr><td>4.7</td><td>0.149</td><td>0.00</td><td>0.143</td><td>0.00</td><td></td><td></td><td></td><td></td></tr><tr><td rowspan="6">CJ</td><td>5.0</td><td>0.146</td><td>0.00</td><td>0.139</td><td>0.00</td><td></td><td></td><td></td><td></td></tr><tr><td>7.5</td><td>0.121</td><td>0.00</td><td>0.108</td><td>0.00</td><td></td><td></td><td></td><td></td></tr><tr><td>10.0</td><td>0.0892</td><td>0.00</td><td>0.0709</td><td>0.00</td><td></td><td></td><td></td><td></td></tr><tr><td>12.5</td><td>0.0593</td><td>0.00</td><td>0.0371</td><td><strong>0.13</strong></td><td></td><td></td><td></td><td></td></tr><tr><td>0.0</td><td></td><td></td><td>0.203</td><td>0.00</td><td></td><td></td><td></td><td></td></tr><tr><td>2.5</td><td>0.193</td><td>0.00</td><td>0.190</td><td>0.00</td><td></td><td></td><td></td><td></td></tr><tr><td>4.7</td><td>0.171</td><td>0.00</td><td>0.163</td><td>0.00</td><td></td><td></td><td></td><td></td></tr><tr><td rowspan="6">JC</td><td>5.0</td><td>0.167</td><td>0.00</td><td>0.157</td><td>0.00</td><td></td><td></td><td></td><td></td></tr><tr><td>7.5</td><td>0.128</td><td>0.00</td><td>0.115</td><td>0.00</td><td></td><td></td><td></td><td></td></tr><tr><td>10.0</td><td>0.0893</td><td>0.00</td><td>0.0657</td><td>0.00</td><td></td><td></td><td></td><td></td></tr><tr><td>12.5</td><td>0.0455</td><td><strong>0.03</strong></td><td>0.0259</td><td><strong>0.51</strong></td><td></td><td></td><td></td><td></td></tr><tr><td>0.0</td><td></td><td></td><td>0.0709</td><td>0.00</td><td></td><td></td><td>9.12</td><td>0.00</td></tr><tr><td>2.5</td><td>0.0673</td><td>0.00</td><td>0.0650</td><td>0.00</td><td>8.27</td><td>0.00</td><td>7.93</td><td>0.00</td></tr><tr><td>4.7</td><td>0.0586</td><td>0.00</td><td>0.0532</td><td><strong>0.01</strong></td><td>6.62</td><td>0.00</td><td>5.62</td><td>0.00</td></tr><tr><td rowspan="5">TK</td><td>5.0</td><td>0.0568</td><td>0.00</td><td>0.0508</td><td><strong>0.01</strong></td><td>6.39</td><td>0.00</td><td>5.30</td><td>0.00</td></tr><tr><td>7.5</td><td>0.0525</td><td><strong>0.01</strong></td><td>0.0464</td><td><strong>0.03</strong></td><td>5.11</td><td>0.00</td><td>3.49</td><td><strong>0.02</strong></td></tr><tr><td>10.0</td><td>0.0560</td><td>0.00</td><td>0.0441</td><td><strong>0.04</strong></td><td>6.42</td><td>0.00</td><td>5.28</td><td>0.00</td></tr><tr><td>12.5</td><td>0.0683</td><td>0.00</td><td>0.0761</td><td>0.00</td><td>12.2</td><td>0.00</td><td>12.8</td><td>0.00</td></tr><tr><td>0.0</td><td></td><td></td><td>0.156</td><td>0.00</td><td></td><td></td><td>68.8</td><td>0.00</td></tr><tr><td rowspan="6">HK</td><td>2.5</td><td>0.148</td><td>0.00</td><td>0.149</td><td>0.00</td><td>62.8</td><td>0.00</td><td>62.5</td><td>0.00</td></tr><tr><td>4.7</td><td>0.133</td><td>0.00</td><td>0.131</td><td>0.00</td><td>49.3</td><td>0.00</td><td>48.8</td><td>0.00</td></tr><tr><td>5.0</td><td>0.130</td><td>0.00</td><td>0.129</td><td>0.00</td><td>47.1</td><td>0.00</td><td>46.6</td><td>0.00</td></tr><tr><td>7.5</td><td>0.102</td><td>0.00</td><td>0.0980</td><td>0.00</td><td>27.8</td><td>0.00</td><td>27.6</td><td>0.00</td></tr><tr><td>10.0</td><td>0.0646</td><td>0.00</td><td>0.0638</td><td>0.00</td><td>11.3</td><td>0.00</td><td>11.6</td><td>0.00</td></tr><tr><td>12.5</td><td>0.0365</td><td><strong>0.14</strong></td><td>0.0404</td><td><strong>0.08</strong></td><td>2.22</td><td><strong>0.07</strong></td><td>2.88</td><td><strong>0.03</strong></td></tr></tbody></table>

## Global tests

In a first series of tests, we did not rely on the actual errors provided with the data. Rather, for each model tested, we successively assumed values of $\alpha_{95} = 0^{\circ}, 5^{\circ}, 7.5^{\circ}, 10^{\circ}$, and $12.5^{\circ}$. This made it possible to again run both the ‘exact’ and ‘approximate’ methods in parallel for further validation of the ‘approximate’ method. It also provided us with a quick way of assessing which GGP model could possibly be considered compatible with the test data set. Table 5 shows the results of those tests.

Interestingly, we first note some common features with the results of the previous test (recall Table 4). Model QC is rejected at more than a 95 per cent confidence level for $\alpha_{95} = 0^{\circ}, 10^{\circ}$, and $12.5^{\circ}$, but not for $\alpha_{95} = 5^{\circ}$ and $7.5^{\circ}$. Again we see that model CP shows a tendency to be compatible with the data if one assumes $\alpha_{95} = 10^{\circ}$, though this possibility is now clearly less likely. Finally we see that only by assuming $\alpha_{95} = 12.5^{\circ}$ can models CJ, JC, HK possibly be consistent with the data, while model TK is always rejected at a 99 per cent confidence level. Based on these tests, we therefore reach the conclusion that if we were to ignore the actual level of error
<table><caption>**Table 5.** Testing all models against the real test data set, assuming various error levels (this error being assumed to be the same for each data of the test data set). Results are presented as in Table 4.</caption><thead><tr><th rowspan="3">model</th><th rowspan="3">error</th><th colspan="4">KS</th><th colspan="4">AD</th></tr><tr><th colspan="2">approx.</th><th colspan="2">exact</th><th colspan="2">approx.</th><th colspan="2">exact</th></tr><tr><th>M<sub>N</sub></th><th>P(M<sub>N</sub>)</th><th>M<sub>N</sub></th><th>P(M<sub>N</sub>)</th><th>I<sub>N</sub></th><th>P(I<sub>N</sub>)</th><th>I<sub>N</sub></th><th>P(I<sub>N</sub>)</th></tr></thead><tbody><tr><td rowspan="6">CP</td><td>0.0</td><td></td><td></td><td>0.0978</td><td>0.00</td><td></td><td></td><td>22.3</td><td>0.00</td></tr><tr><td>5.0</td><td>0.0780</td><td>0.00</td><td>0.0787</td><td>0.00</td><td>13.2</td><td>0.00</td><td>13.9</td><td>0.00</td></tr><tr><td>7.5</td><td>0.0574</td><td>0.00</td><td>0.0605</td><td>0.00</td><td>6.83</td><td>0.00</td><td>7.90</td><td>0.00</td></tr><tr><td>10.0</td><td>0.0406</td><td><strong>0.07</strong></td><td>0.0430</td><td><strong>0.05</strong></td><td>4.49</td><td><strong>0.01</strong></td><td>5.83</td><td>0.00</td></tr><tr><td>12.5</td><td>0.0757</td><td>0.00</td><td>0.0752</td><td>0.00</td><td>9.16</td><td>0.00</td><td>10.5</td><td>0.00</td></tr><tr><td>0.0</td><td></td><td></td><td>0.0481</td><td><strong>0.02</strong></td><td></td><td></td><td>4.41</td><td><strong>0.01</strong></td></tr><tr><td rowspan="6">QC</td><td>5.0</td><td>0.0365</td><td><strong>0.14</strong></td><td>0.0382</td><td><strong>0.11</strong></td><td>1.60</td><td><strong>0.16</strong></td><td>1.48</td><td><strong>0.18</strong></td></tr><tr><td>7.5</td><td>0.0247</td><td><strong>0.58</strong></td><td>0.0289</td><td><strong>0.37</strong></td><td>0.648</td><td><strong>0.60</strong></td><td>1.11</td><td><strong>0.30</strong></td></tr><tr><td>10.0</td><td>0.0435</td><td><strong>0.05</strong></td><td>0.0555</td><td>0.00</td><td>2.97</td><td><strong>0.03</strong></td><td>5.00</td><td>0.00</td></tr><tr><td>12.5</td><td>0.0715</td><td>0.00</td><td>0.0861</td><td>0.00</td><td>10.5</td><td>0.00</td><td>15.1</td><td>0.00</td></tr><tr><td>0.0</td><td></td><td></td><td>0.154</td><td>0.00</td><td></td><td></td><td></td><td></td></tr><tr><td>5.0</td><td>0.131</td><td>0.00</td><td>0.124</td><td>0.00</td><td></td><td></td><td></td><td></td></tr><tr><td rowspan="6">CJ</td><td>7.5</td><td>0.107</td><td>0.00</td><td>0.0957</td><td>0.00</td><td></td><td></td><td></td><td></td></tr><tr><td>10.0</td><td>0.0749</td><td>0.00</td><td>0.0574</td><td>0.00</td><td></td><td></td><td></td><td></td></tr><tr><td>12.5</td><td>0.0601</td><td>0.00</td><td>0.0419</td><td><strong>0.06</strong></td><td></td><td></td><td></td><td></td></tr><tr><td>0.0</td><td></td><td></td><td>0.168</td><td>0.00</td><td></td><td></td><td></td><td></td></tr><tr><td>5.0</td><td>0.136</td><td>0.00</td><td>0.129</td><td>0.00</td><td></td><td></td><td></td><td></td></tr><tr><td>7.5</td><td>0.105</td><td>0.00</td><td>0.0898</td><td>0.00</td><td></td><td></td><td></td><td></td></tr><tr><td rowspan="6">JC</td><td>10.0</td><td>0.0737</td><td>0.00</td><td>0.0470</td><td><strong>0.02</strong></td><td></td><td></td><td></td><td></td></tr><tr><td>12.5</td><td>0.0530</td><td><strong>0.01</strong></td><td>0.0225</td><td><strong>0.69</strong></td><td></td><td></td><td></td><td></td></tr><tr><td>0.0</td><td></td><td></td><td>0.0642</td><td>0.00</td><td></td><td></td><td>8.53</td><td>0.00</td></tr><tr><td>5.0</td><td>0.0566</td><td>0.00</td><td>0.0529</td><td><strong>0.01</strong></td><td>6.34</td><td>0.00</td><td>5.17</td><td>0.00</td></tr><tr><td>7.5</td><td>0.0630</td><td>0.00</td><td>0.0567</td><td>0.00</td><td>5.64</td><td>0.00</td><td>3.98</td><td><strong>0.01</strong></td></tr><tr><td>10.0</td><td>0.0706</td><td>0.00</td><td>0.0602</td><td>0.00</td><td>7.65</td><td>0.00</td><td>6.56</td><td>0.00</td></tr><tr><td rowspan="6">TK</td><td>12.5</td><td>0.0844</td><td>0.00</td><td>0.0701</td><td>0.00</td><td>14.0</td><td>0.00</td><td>15.0</td><td>0.00</td></tr><tr><td>0.0</td><td></td><td></td><td>0.156</td><td>0.00</td><td></td><td></td><td>82.3</td><td>0.00</td></tr><tr><td>5.0</td><td>0.134</td><td>0.00</td><td>0.133</td><td>0.00</td><td>56.5</td><td>0.00</td><td>54.4</td><td>0.00</td></tr><tr><td>7.5</td><td>0.107</td><td>0.00</td><td>0.103</td><td>0.00</td><td>34.2</td><td>0.00</td><td>31.8</td><td>0.00</td></tr><tr><td>10.0</td><td>0.0689</td><td>0.00</td><td>0.0703</td><td>0.00</td><td>14.7</td><td>0.00</td><td>13.1</td><td>0.00</td></tr><tr><td>12.5</td><td>0.0409</td><td><strong>0.07</strong></td><td>0.0442</td><td><strong>0.04</strong></td><td>2.89</td><td><strong>0.03</strong></td><td>2.95</td><td><strong>0.03</strong></td></tr></tbody></table>

affecting the data in the test data set, all models would be rejected at more than a 95 per cent confidence level except model QC when assuming $\alpha_{95} = 5^{\circ}$ or $7.5^{\circ}$, or model CJ and JC when assuming $\alpha_{95} = 12.5^{\circ}$.

The key point, now, is that we do have some knowledge of the errors affecting the data in the test data set. We already noted that on average, this error is of $\alpha_{95} = 4.7^{\circ}$. This value is clearly incompatible with the $\alpha_{95} = 12.5^{\circ}$ value needed for models CJ and JC to possibly be compatible with the test data set. However, it seems compatible with the results for model QC.

To finally assess the issue, we eventually ran the test for all six models, no longer assuming the same *a priori* error for all the data of the test data set, but the true error associated (i.e. published) with each of the 990 data. Since, as already noted, this requires considerably more computational time, only the approximate method was used for this comparative assessment. (Note that all previous tests have anyway shown that for errors of the order of $5^{\circ}$, the approximate method always leads to accurate enough results.) Table 6 gives the results of this final test, and Fig. 2 shows the CDFs of the corresponding uniformized data. Those results clearly show that only model QC is not rejected by the test data set. All other models can be rejected at more than a 99 per cent confidence level, TK being the second, CP the third, HK the fourth, CJ the fifth and JC the sixth best models in terms of both the $M_N$ and $I_N$ measures (note that

<table><caption>**Table 6.** Testing all models against the real test data set, using the true (published) error for each data. Those tests have been run only using the ‘approximate’ method except in the case of model QC for which the ‘exact’ method was also used (last line).</caption><thead><tr><th rowspan="3">model</th><th colspan="2">KS</th><th colspan="2">AD</th></tr><tr><th>M<sub>N</sub></th><th>P(M<sub>N</sub>)</th><th>I<sub>N</sub></th><th>P(I<sub>N</sub>)</th></tr><tr><th colspan="4">'Approximate' method</th></tr></thead><tbody><tr><td>CP</td><td>0.0773</td><td>0.00</td><td>13.1</td><td>0.00</td></tr><tr><td>QC</td><td>0.0319</td><td><strong>0.26</strong></td><td>1.48</td><td><strong>0.18</strong></td></tr><tr><td>CJ</td><td>0.131</td><td>0.00</td><td></td><td></td></tr><tr><td>JC</td><td>0.138</td><td>0.00</td><td></td><td></td></tr><tr><td>TK</td><td>0.0552</td><td>0.00</td><td>6.38</td><td>0.00</td></tr><tr><td>HK</td><td>0.125</td><td>0.00</td><td>54.3</td><td>0.00</td></tr><tr><td colspan="5">'Exact' method</td></tr><tr><td>QC</td><td>0.0322</td><td><strong>0.25</strong></td><td>1.46</td><td><strong>0.19</strong></td></tr></tbody></table>

when both quantities are available in Table 6, $M_N$ and $I_N$ lead to the same relative ranking).

As an ultimate check and only for model QC, we finally duplicated this last test using the ‘exact’ method. The result of this test (also shown in Table 6) confirms the previous results. Moreover, the difference found between those two results for
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Cumulative distribution functions (CDFs) of the N = 990 uniformized data {$t_i$} for each of the six models, when testing them against the real test data set, assuming the true errors and using the 'approximate' method. Those CDFs lead to the $M_N$, $P(M_N)$, $I_N$, $P(I_N)$ shown in Table 6. In the case of model QC, the CDF obtained when using the 'exact' method leads to a CDF almost indistinguishable from the one shown here. Color code: CP(red), QC(green), CJ(dark blue), JC(pink), TK(light blue), HK(black).</figcaption></figure>

model QC clearly appears to be much smaller than the differences among the results of the approximate method for all models. We, therefore, feel quite confident that these results can be trusted.

It is interesting to take advantage of model QC, the only model not rejected by the test data set so far, to also quickly address the issue of the number of Gauss coefficients that should be taken into account to carry on the tests. As stated earlier, we used all Gauss coefficients up to degree and order 7. However, all models, including QC, assume non-zero variances for all degrees up to infinity. However, the higher the degree n, the weaker the contribution of the coefficients to the local PDFs $p_i(u)$, and thus to the final result. Table 7 illustrates this point and shows that the calculation is indeed already converged when taking into account all degrees up to 7. Note that this also shows that the test (applied to the test data set used in the present study) is, therefore, hardly sensitive to Gauss coefficients of the GGP models with degrees above 7.

## Local tests

It is also important to keep in mind that whereas the results presented so far show that models CP, CJ, JC TK and HK should be rejected, they do not imply that model QC is fully compatible with the data. Those results only correspond to a global test, for which all the 990 data $u_i$ have been uniformized into a single global distribution of {$t_i$} over [0, 1], and it may be that model QC, not rejected by this global test, could be rejected by a more local test. Table 8 shows results of such tests. The same testing procedure has been used as before (with both the exact and approximate method, using both the KS and AD tests and assuming the true error), except for the fact

<table><caption>**Table 7.** Convergence test. Model QC is tested here against the real test data set, using the true errors and the 'approximate' method. The test is repeated eight times, taking into account more and more of the Gauss coefficients, up to successively degree 1, 2, 3, 4, 5, 6, 7, and finally 8.</caption><thead><tr><th rowspan="2">degree</th><th colspan="2">KS</th><th colspan="2">AD</th></tr><tr><th><i>M</i><sub>N</sub></th><th><i>P</i>(<i>M</i><sub>N</sub>)</th><th><i>I</i><sub>N</sub></th><th><i>P</i>(<i>I</i><sub>N</sub>)</th></tr></thead><tbody><tr><td>1</td><td>0.362</td><td>0.00</td><td>645.</td><td>0.00</td></tr><tr><td>2</td><td>0.0808</td><td>0.00</td><td>21.3</td><td>0.00</td></tr><tr><td>3</td><td>0.0467</td><td><b>0.03</b></td><td>4.36</td><td><b>0.01</b></td></tr><tr><td>4</td><td>0.0356</td><td><b>0.16</b></td><td>2.13</td><td><b>0.08</b></td></tr><tr><td>5</td><td>0.0330</td><td><b>0.23</b></td><td>1.64</td><td><b>0.15</b></td></tr><tr><td>6</td><td>0.0323</td><td><b>0.25</b></td><td>1.52</td><td><b>0.17</b></td></tr><tr><td>7</td><td>0.0319</td><td><b>0.26</b></td><td>1.48</td><td><b>0.18</b></td></tr><tr><td>8</td><td>0.0319</td><td><b>0.26</b></td><td>1.47</td><td><b>0.18</b></td></tr></tbody></table>

that for each test, only the data from one of the 36 sites of the test data set have been used at a time, to produce the (then much less numerous) uniformized {$t_i$}.

Those results show that the KS-test reveals four sites (sites 17, 18, 30 and 35) rejecting model QC at a 99 per cent confidence level because of too low values of $P(M_N)$. This amounts to about 10 per cent of the sites, an uncomfortably high number. By contrast, the AD-test only reveals two such sites (sites 3 and 18). However, it confirms the strong level of disagreement of site 18 with model QC. This site clearly needs to be double-checked. Indeed, a closer look at sites 17 and 18, reveals some interesting features. Those sites are geographically very close to one another, at latitude 50.0° and longitude 7.0°. Site 17 (data selected from Böhnel et al. 1987), is in the quaternary East-Eifel volcanic field, whereas site 18 (data selected from Böhnel et al. 1982), is in the quaternary West-Eifel volcanic field. For those two sites, we show in Fig. 3 how the selected data $u_i$ plot on the unit sphere $\Sigma$. On the same plots, we also show the isovalues of the $p_i(u)$, which have been computed assuming $\alpha_{95} = 2.76^\circ$ for site 17 and $\alpha_{95} = 4.68^\circ$ for site 18. Those values correspond to the average errors for all data from respectively sites 17 and 18. (Using such averages is done here mainly for illustrative purpose, the plotted isovalues being then simultaneously meaningful for all the data shown on the same plot). Also shown in Fig. 3 are the CDFs for the {$t_i$} corresponding to this calculation, and the analogous CDFs corresponding to the exact calculation (i.e. when the true, rather than the average, error is being used for each individual data). The reasons sites 17 and 18 reject model QC are clear: for site 17, the data tend to display too little dispersion compared to the expected distribution (producing too many low values for the {$t_i$}), whereas for site 18, the data show a strong tendency to lie East of the expected distribution (producing too many high values for the {$t_i$}).

The interesting observation here is that beyond the disagreement of each site with model QC, we see that both sites also disagree with each other. Since they share virtually the same location, this would suggest that two different stable polarity data sets acquired at the same location but at different times during the Bruhnes chron, could disagree. This is quite an embarrassing finding in view of the stationarity assumption underlying the GGP modelling approach. Fortunately, a closer look at the Böhnel et al. (1982, 1987) papers provides us with an explanation for this disagreement. As far as site 17 is concerned, it simply appears that the 44 data selected in the Quidelleur et al. (1994) database are not fully independent. In seven instances, several data points appear to have basically recorded
<table><caption>**Table 8.** Testing model QC against the real test data set on a site by site basis. The true errors and both the ‘approximate’ and ‘exact’ methods are being used. For each site, results are provided for both the KS and AD tests.</caption><thead><tr><th rowspan="3">Site</th><th colspan="4">KS</th><th colspan="4">AD</th></tr><tr><th colspan="2">Approx.</th><th colspan="2">Exact</th><th colspan="2">Approx.</th><th colspan="2">Exact</th></tr><tr><th>M<sub>N</sub></th><th>P(M<sub>N</sub>)</th><th>M<sub>N</sub></th><th>P(M<sub>N</sub>)</th><th>I<sub>N</sub></th><th>P(I<sub>N</sub>)</th><th>I<sub>N</sub></th><th>P(I<sub>N</sub>)</th></tr></thead><tbody><tr><td>1</td><td>0.106</td><td>0.15</td><td>0.105</td><td>0.16</td><td>1.81</td><td>0.12</td><td>1.73</td><td>0.13</td></tr><tr><td>2</td><td>0.149</td><td>0.84</td><td>0.163</td><td>0.75</td><td>0.626</td><td>0.62</td><td>0.598</td><td>0.65</td></tr><tr><td>3</td><td>0.284</td><td>0.05</td><td>0.286</td><td>0.05</td><td>3.63</td><td>0.01</td><td>3.63</td><td>0.01</td></tr><tr><td>4</td><td>0.106</td><td>0.86</td><td>0.105</td><td>0.86</td><td>0.485</td><td>0.76</td><td>0.478</td><td>0.77</td></tr><tr><td>5</td><td>0.122</td><td>0.90</td><td>0.115</td><td>0.93</td><td>0.576</td><td>0.67</td><td>0.542</td><td>0.70</td></tr><tr><td>6</td><td>0.162</td><td>0.82</td><td>0.162</td><td>0.82</td><td>1.01</td><td>0.35</td><td>0.991</td><td>0.36</td></tr><tr><td>7</td><td>0.0745</td><td>0.98</td><td>0.0762</td><td>0.98</td><td>0.401</td><td>0.85</td><td>0.413</td><td>0.83</td></tr><tr><td>8</td><td>0.288</td><td>0.16</td><td>0.287</td><td>0.16</td><td>2.96</td><td>0.03</td><td>2.92</td><td>0.03</td></tr><tr><td>9</td><td>0.265</td><td>0.15</td><td>0.264</td><td>0.16</td><td>1.28</td><td>0.24</td><td>1.27</td><td>0.24</td></tr><tr><td>10</td><td>0.208</td><td>0.05</td><td>0.211</td><td>0.05</td><td>3.46</td><td>0.02</td><td>3.33</td><td>0.02</td></tr><tr><td>11</td><td>0.167</td><td>0.27</td><td>0.165</td><td>0.29</td><td>1.46</td><td>0.19</td><td>1.46</td><td>0.19</td></tr><tr><td>12</td><td>0.199</td><td>0.50</td><td>0.198</td><td>0.51</td><td>0.469</td><td>0.78</td><td>0.458</td><td>0.79</td></tr><tr><td>13</td><td>0.156</td><td>0.32</td><td>0.158</td><td>0.30</td><td>2.08</td><td>0.08</td><td>2.23</td><td>0.07</td></tr><tr><td>14</td><td>0.0867</td><td>0.80</td><td>0.0878</td><td>0.79</td><td>0.686</td><td>0.57</td><td>0.625</td><td>0.62</td></tr><tr><td>15</td><td>0.304</td><td>0.38</td><td>0.325</td><td>0.30</td><td>2.71</td><td>0.04</td><td>2.88</td><td>0.03</td></tr><tr><td>16</td><td>0.116</td><td>0.66</td><td>0.108</td><td>0.74</td><td>0.876</td><td>0.43</td><td>0.789</td><td>0.49</td></tr><tr><td>17</td><td>0.234</td><td>0.01</td><td>0.233</td><td>0.01</td><td>2.81</td><td>0.03</td><td>2.83</td><td>0.03</td></tr><tr><td>18</td><td>0.287</td><td>0.00</td><td>0.288</td><td>0.00</td><td>8.67</td><td>0.00</td><td>8.76</td><td>0.00</td></tr><tr><td>19</td><td>0.185</td><td>0.12</td><td>0.195</td><td>0.09</td><td>2.31</td><td>0.06</td><td>2.49</td><td>0.05</td></tr><tr><td>20</td><td>0.136</td><td>0.55</td><td>0.136</td><td>0.54</td><td>0.618</td><td>0.63</td><td>0.614</td><td>0.63</td></tr><tr><td>21</td><td>0.212</td><td>0.07</td><td>0.201</td><td>0.10</td><td>1.88</td><td>0.11</td><td>1.67</td><td>0.14</td></tr><tr><td>22</td><td>0.200</td><td>0.28</td><td>0.198</td><td>0.30</td><td>1.77</td><td>0.12</td><td>1.85</td><td>0.11</td></tr><tr><td>23</td><td>0.199</td><td>0.51</td><td>0.199</td><td>0.50</td><td>1.08</td><td>0.32</td><td>1.08</td><td>0.32</td></tr><tr><td>24</td><td>0.147</td><td>0.62</td><td>0.158</td><td>0.52</td><td>1.79</td><td>0.12</td><td>1.91</td><td>0.10</td></tr><tr><td>25</td><td>0.198</td><td>0.21</td><td>0.199</td><td>0.21</td><td>2.41</td><td>0.06</td><td>2.48</td><td>0.05</td></tr><tr><td>26</td><td>0.0941</td><td>1.00</td><td>0.0916</td><td>1.00</td><td>0.301</td><td>0.94</td><td>0.315</td><td>0.93</td></tr><tr><td>27</td><td>0.222</td><td>0.65</td><td>0.239</td><td>0.56</td><td>0.796</td><td>0.48</td><td>0.846</td><td>0.45</td></tr><tr><td>28</td><td>0.368</td><td>0.08</td><td>0.373</td><td>0.07</td><td>2.64</td><td>0.04</td><td>2.79</td><td>0.04</td></tr><tr><td>29</td><td>0.140</td><td>0.50</td><td>0.141</td><td>0.50</td><td>0.908</td><td>0.41</td><td>0.994</td><td>0.36</td></tr><tr><td>30</td><td>0.364</td><td>0.01</td><td>0.357</td><td>0.01</td><td>2.52</td><td>0.05</td><td>2.43</td><td>0.05</td></tr><tr><td>31</td><td>0.327</td><td>0.30</td><td>0.321</td><td>0.32</td><td>1.64</td><td>0.15</td><td>1.64</td><td>0.15</td></tr><tr><td>32</td><td>0.179</td><td>0.76</td><td>0.179</td><td>0.76</td><td>1.45</td><td>0.19</td><td>1.40</td><td>0.20</td></tr><tr><td>33</td><td>0.209</td><td>0.78</td><td>0.210</td><td>0.77</td><td>0.925</td><td>0.40</td><td>0.927</td><td>0.40</td></tr><tr><td>34</td><td>0.273</td><td>0.38</td><td>0.274</td><td>0.38</td><td>1.47</td><td>0.19</td><td>1.45</td><td>0.19</td></tr><tr><td>35</td><td>0.248</td><td>0.01</td><td>0.240</td><td>0.01</td><td>3.42</td><td>0.02</td><td>3.23</td><td>0.02</td></tr><tr><td>36</td><td>0.119</td><td>0.91</td><td>0.112</td><td>0.94</td><td>0.480</td><td>0.77</td><td>0.455</td><td>0.79</td></tr></tbody></table>

the same field (each time reflecting a single volcanic event). Those events correspond to low values of $t_i$ and bias the distribution of the {$t_i$} towards low values. Site 18 is a slightly different story. In that case, Böhnel et al. (1987) make a good case that a subset of the data they published must have recorded a relatively short-lived excursion or event during the Bruhnes chron, biasing the data set towards the East. Although Quidelleur et al. (1994) rejected most of those transitional data (based on the selection criteria that the data must correspond to a period of stable polarity) some of the data biased towards that transitional direction have clearly been included within the database. Thus the issue with sites 17 and 18 appears to be more one of ill sampling than one of overall disagreement with model QC. Interestingly, we further note that if we actually bin the two sites into a single one, as is also shown in Fig. 3, the local test no

longer rejects model QC. This then leaves only two sites (30 and 35) out of 36 truly rejecting model QC at a 99 per cent confidence level because of too low values of $P(M_N)$, and only one site (site 3) rejecting QC because of a too low value of $P(I_N)$. These are no longer unexpected proportions and model QC thus indeed appears to be compatible with the Quidelleur et al. (1994) Bruhnes stable polarity data set.

CONCLUSION

In the present paper we introduced the first quantitative method capable of assessing the compatibility of a GGP model with a given palaeodirectional database in a way consistent with both the
<figure><img src="image_4.png" /></figure>

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

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

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

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

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

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

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

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

**Figure 3.** Local tests for model QC at sites 17 (a), 18 (b), and both sites combined (c). Left column shows the way the directional data plots on the unit sphere Σ which is Lambert-projected, the centre point (pointing downwards through the page) corresponding to the direction produced by a pure axial dipole at the site; north is towards the top of the page, East towards the right. Angular distances with respect to the centre of the plot are shown every 5° (green circles). Isovalues of the *p*<sub>*i*</sub>(*u*) (computed with α<sub>95</sub> = 2.76° for site 17 (a), α<sub>95</sub> = 4.68° for site 18 (b), and α<sub>95</sub> = 3.72° for both sites combined (c)) are also shown as red ellipses. Those correspond to the bounds within which, respectively, 10 per cent, 30 per cent, 50 per cent, 70 per cent, 90 per cent of the data points should plot. Central column shows the CDF of the corresponding uniformized data {t<sub>i</sub>}.
Right column shows the analogous CDF when the true errors are used.

statistical assumptions underlying the GGP approach, and the nature
of the measurement error affecting palaeodirectional data.

The method was successfully implemented and tested, in its exact
form, but also in an approximate form, relying on some useful ap-
proximations which we introduced to significantly reduce the com-
putational burden of the tests. The limits of the approximate method
have been assessed. Those limits are reached when both the GGP
model and the measurement errors introduce strong dispersions of
the directional data on the unit sphere. For realistic measurement
errors, this approximate method, however, turns out to be just as
good as the exact method. This brought many encouraging results.

First, it was shown that the method is usefully discriminating. A
GGP model is compatible with an artificial data set affected by a
given (known) level of measurement error only if the model tested
is the same as the model used to produce the artificial data and if the

correct level of error affecting the data is assumed when carrying the
test. However, if one assumes either too small or too large an error
affecting the data, or tests a GGP model different from the starting
model, this conclusion may no longer be reached. The method we
propose thus appears to be relevant for deciding which GGP model,
if any, best describes the statistical behaviour of the past geomagnetic
field.

We applied the method to test six published GGP models against
palaeodirectional Bruhnes stable polarity data from the Quidelleur
et al. (1994) database. All but one model (model QC) had to be re-
jected. This shows that in practice also, our method can successfully
discriminate GGP models attempting to describe a given palaeomag-
netic data set.

The method we propose further offers the possibility of assessing
a GGP model at any scale, local, regional or global. This makes it
possible to better scrutinize a GGP model which first passed tests at a global scale. Any disagreement between the prediction of the GGP model and the data set at a local level can then be used to double-check not only the GGP model, but also the data set itself. Applying this checking procedure to model QC allowed us to identify such problems with two sites of the Quidelleur et al. (1994) database and to confirm the compatibility of this model with the rest of the test data set.

This success of at least one GGP model is motivating. It shows, for the first time, that the statistical behaviour of the geomagnetic field at times of stable polarity can indeed be described in a consistent way in terms of a GGP model. However, the specific success of model QC should not be overemphasized. This model was inferred in an empirical way by Quidelleur & Courtillot (1996) from the Quidelleur et al. (1994) database from which our test data set was also extracted. All other models were inferred (also in empirical ways) from different databases, and even though all those databases share many common data, they may very well not be fully compatible with each other and with our own data set. Part of the failure of those GGP models at being compatible with our test data set might originate from this. With this respect, it is quite clear that finally confirming the success of the QC model (and the failure of other models) would still require a more thorough study, involving a more recent, extensive and independently assessed data set, extracted from databases such as the IAGA palaeomagnetic reference database we mentioned early on.

As a matter of fact, it is interesting to note that such a data set could also include palaeointensity data extracted from databases such as the IAGA palaeointensity reference database (PINT, Perrin & Schnepp (2004)). Indeed, although the present paper mainly described the way to deal with palaeodirectional data, palaeointensity data could easily be taken into account [if such data provide the full 3-D palaeofield at a given location, the way those data plot in 3-D can be compared with the local 3-D Gaussian statistics predicted by the GGP model (i.e. g) by using just the same kind of uniformization procedure as the one we described, see paper I; and taking measurement error into account would just consist in making use of a modified, but still Gaussian, 3-D statistics of the type $g_a$ as defined in the appendix, where $g_e$ is then directly the Gaussian error assumed to affect the data]. In fact such tests might eventually show that model QC also is incompatible with such a more complete data set. However, if such is the case, the method we propose here could be used to seek yet another model. Several strategies could then be pursued. One could involve using the same kind of empirical procedure as those used by Quidelleur & Courtillot (1996) to first guess which parameter of the GGP model should be changed, and using our method to assess the improvement brought. Another strategy, more advanced and much more far fetched at this point, could involve designing a more systematic 'inverse' approach based on the present 'forward' testing method. In either case, however, several issues would still have to be faced. In particular, we mentioned the fact that the method is not sensitive to parameters corresponding to degrees above 7, which reflects the weak contribution of those high degrees to the geomagnetic field observed at the Earth's surface. However, we did not investigate in detail the sensitivity of our test to the various parameters defining the mean field ($E(\mathbf{k})$) and the covariance matrix (Cov($\mathbf{k}, \mathbf{k}$)) of a GGP model. Assessing this would help us better understand how complex a GGP model, and in particular its mean field, really needs to be to explain the data. For the time being, and as far as the present test study can suggest, the Bruhnes stable polarity data does not seem to call for more complex a model than model QC: a mean field with only

$g_1^0$ and a weak $g_2^0$ contribution, and a simple diagonal covariance matrix only involving an enhanced $\sigma_2^1$ contribution along the lines first suggested by Kono & Tanaka (1995) and Hulot & Gallet (1996).

We plan to make our software available upon request (please be in touch with the corresponding author: gh@ipgp.jussieu.fr).

## ACKNOWLEDGMENTS

We thank Richard Holme and two anonymous reviewers for their constructive comments. This is IPGP contribution 2147.

## REFERENCES

Bingham, C., 1983. A series expansion for Angular Gaussian Distribution. Appendix C., in *Statistics on Spheres*, ed. Watson, G., Wiley-Interscience Pub., New York.

Böhnel, H., Kohnen, H., Negendank, J. & Schmincke, H.-U., 1982. Palaeomagnetism of Quaternary volcanics of the East-Eifel, Germany, *J. Geophys.*, **51**, 29–37.

Böhnel, H., Reismann, N., Jäger, G., Haverkamp, U., Negendank, J.F.W. & Schmincke, H.-U., 1987. Paleomagnetic investigation of Quaternary West Eifel volcanics (Germany): indication for increased volcanic activity during geomagnetic excursion/event?, *J. Geophys.*, **62**, 50–61.

Bouligand, C., Hulot, G., Khokhlov, A. & Glatzmaier, G.A., 2005. Statistical paleomagnetic field modelling and dynamo numerical simulation, *Geophys. J. Int.*, **161**, 603–626.

Constable, C.G. & Parker, R.L., 1988. Statistics of the geomagnetic secular variation for the past 5 Myr, *J. geophys. Res.*, **93**, 11 569–11 581.

Constable, C.G. & Johnson, C.L., 1999. Anisotropic paleosecular variation models: implications for geomagnetic field observables, *Earth planet. Sci. Lett.*, **115**, 35–51.

Hatakeyama, T. & Kono, M., 2001. Shift of the mean magnetic field values; effect of scatter due to secular variation and errors, *Earth Planets Space*, **53**, 31–44.

Hatakeyama, T. & Kono, M., 2002. Geomagnetic field model for the last 5 Myr; time-averaged field and secular variation, Paleosecular variation and reversals of the Earth's magnetic field, *Phys. Earth planet. Int.*, **133**, 181–215.

Hongre, L., Hulot, G. & Khokhlov, A., 1998. An analysis of the geomagnetic field over the past 2000 years, *Phys. Earth planet. Int.*, **106**, 311–315.

Hulot, G. & Gallet, Y., 1996. On the interpretation of virtual geomagnetic pole (VGP) scatter curves, *Phys. Earth planet. Int.*, **95**, 37–53.

Hulot, G. & Le Mouël, J.-L., 1994. A statistical approach to the Earth's main magnetic field, *Phys. Earth planet. Int.*, **82**, 167–183.

Hulot, G. & Bouligand, C., 2005. Statistical paleomagnetic field modelling and symmetry considerations, *Geophys. J. Int.*, **161**, 591–602.

Hulot, G., Eymin, C., Langlais, B., Mandea, M. & Olsen, N., 2002. Smallscale structure of the geodynamo inferred from Oersted and Magsat satellite data, *Nature*, **416**, 620–623.

Jackson, A., Jonkers, A.R.T. & Walker, M.R., 2000. Four centuries of geomagnetic secular variation from historical records, *Phil. Trans. R. Soc. Lond.*, A, **358**, 957–990.

Johnson, C. & Constable, C., 1996. Paleosecular variation recorded by lava flows over the last 5 Myr, *Phil. Trans. R. Soc. Lond.*, **354**, 89–141.

Khokhlov, A., Hulot, G. & Carlut, J., 2001. Towards a self-consistent approach to paleomagnetic field modelling, *J. Int. Geoph.*, **145**, 157–171.

Kono, M. & Tanaka, H., 1995. Mapping the Gauss coefficients to the pole and the models of paleosecular variation, *J. Geomag. Geoelectr.*, **47**, 115–130.

Kono, M., Tanaka, H. & Tsukanawa, H., 2000a. Spherical harmonic analysis of paleomagnetic data: the case of linear mapping, *J. geophys. Res.*, **105**, 5817–5833.

Kono, M., Sakuraba, A. & Ishida, M., 2000b. Dynamo simulations and paleosecular variation models, *Phil. Trans. R. Soc. Lond.*, A, **358**, 1123–1139.
Korte, M. & Constable, C.G., 2005. Continuous geomagnetic field models for the past 7 millennia: 2. CALS7K *Geochem. Geophys. Geosys.*, **6**, 2004GC000801.

Love, J.J. & Constable, C.G., 2003. Gaussian statistics for paleomagnetic vectors, *Geophys. J. Int.*, **152**, 515–565.

Marsaglia, G. & Marsaglia, J.C.W., 2004. Evaluating the Anderson-Darling distribution, *Journal Statistical Software*, **9**, Issue 5.

McElhinny, M. & McFadden, P., 1997. Paleosecular variation over the past 5Myr based on a new generalized database, *Geophys. J. Int.*, **131**, 240–252.

McElhinny, M. & McFadden, P., 2000. Paleomagnetism, Continents and oceans, *International Geophysics Series*, **73**, Academic Press, San Diego, CA, USA.

McMillan, D.G., Constable, C.G., Parker, R.L. & Glatzmaier, G.A., 2001. A statistical analysis of magnetic fields from some geodynamo simulations, *Geochem. Geophys. Geosys.*, **2**, 2000GC000130.

Merrill, R., McElhinny, M. & McFadden, P., 1996. *The Magnetic Field of the Earth*, Academic Press, London.

Perrin, M. & Schneppe, E., 2004. IAGA paleointensity database: distribution and quality of the data set, *Phys. Earth planet. Int.*, **147**, 255–267.

Press, C., Teukolsky, S., Vetterling, W. & Flannery, B., 1996. Numerical Recipes in C. 2nd edn, Cambridge, see also http://library.lanl.gov/numerical/index.html.

Quidelleur, X. & Courtillot, V., 1996. On low-degree spherical harmonic models of paleosecular variation, *Phys. Earth planet. Int.*, **95**, 55–77.

Quidelleur, X., Valet, J.-P., Courtillot, V. & Hulot, G., 1994. Long-term geometry of the geomagnetic field for the last five million years: an updated secular variation database, *Geophys. Res. Lett.*, **15**, 1639–1642.

Tanaka, H., Kono, M. & Uchimura, H., 1995. Some global features of paleointensity in geological time, *Geophys. J. Int.*, **120**, 97–102.

Tauxe, L., 2005. Inclination flattening and the geocentric axial dipole, *Earth planet. Sci. Lett.*, **233**, 247–261.

Tauxe, L. & Kent, D.V., 2004. A simplified statistical model for the geomagnetic field and the detection of shallow bias in paleomagnetic inclinations: was the ancient magnetic field dipolar?, in *Timescales of the Paleomagnetic field*, Vol. 145, pp. 101–115. eds Channell, J.E.T., Kent, D.V., Lowrie, W. & Meert, J., Am. Geophys. Un. Monongr.

# APPENDIX: INTRODUCING USEFUL ANALYTIC APPROXIMATIONS

## Fisher and Angular Gaussian distributions

Let us first note that a Fisher distribution $t_K(\mathbf{u}, \mathbf{s})$ centred on $\mathbf{u} \in \Sigma$ with concentration parameter $K$ is very close (when $K$ is large enough) to the angular Gaussian distribution derived from a 3-D Gaussian distribution characterized by $\mathbf{m}$ and $\Lambda^{ij} = \sigma^{-2}\delta_{ij}$, provided $\mathbf{u} = \mathbf{m}/|\mathbf{m}|$ and $\sigma = |\mathbf{m}|K^{-1/2}$. This fairly intuitive result (see also Love & Constable 2003), is in the same line as the better known result that a Fisher distribution is very close to a 2-D Gaussian distribution within the plane tangent to $\Sigma$ at $\mathbf{u}$ (see e.g. Merrill et al. 1996).

In spherical coordinates, and when assuming the Fisher distribution centred on $\mathbf{u}_0 = (\theta = 0, \varphi = 0)$, from eq. (7) we infer $z =$

$m \cos \theta$ and $m = K^{1/2}$. Then, the (axisymmetric) angular Gaussian distribution (with $\theta$ as natural argument) takes the form:

$$s(\theta) = \frac{K^{1/2} e^{-K/2} \cos \theta}{2\sqrt{2\pi} 3^{3/2}} + \frac{1}{4\pi} e^{-K(\sin \theta)^2} (1 + K \cos^2 \theta) \left[ 1 + \text{Erf}\left( \frac{K^{1/2} \cos \theta}{\sqrt{2}} \right) \right] \quad (A1)$$

To leading order when $K$ becomes large, this becomes:

$$s(\theta) \approx \frac{K (\cos \theta)^2}{2\pi} e^{-\frac{1}{2} K (\sin \theta)^2}, \quad (A2)$$

which, to leading order when $\theta$ is small enough, further becomes:

$$s(\theta) \approx \frac{K}{2\pi} e^{-\frac{1}{2} K \theta^2}. \quad (A3)$$

This indeed is the approximate form $t_K(\mathbf{u}_0, \mathbf{s})$ (as defined by 10) takes when $K$ is large enough, and $\theta$ small enough.

To numerically show that for useful values of the concentration parameter $K$, $s(\theta)$ is indeed very close to the $t_K(\mathbf{u}_0, \mathbf{s})$ Fisher distribution, we computed the $P$-probability angles $\alpha(P)$ (defined as the maximum angle the direction may depart from the mean direction $\mathbf{u}_0$ with probability $P$) predicted by $t_K(\mathbf{u}_0, \mathbf{s})$:

$$\alpha_{[P]} = \arccos \left[ \frac{1}{K} \log(e^K + P e^{-K} - P e^K) \right], \quad (A4)$$

and checked that for those angles, the angular Gaussian distribution $s(\theta)$ given by eq. (A1), leads to a correct prediction of $P$. For various relevant values of $K$ and $P$, Table A1 shows values of the departure $\varepsilon(P)$ of the probability predicted by the angular Gaussian distribution from the expected probability $P$:

$$\varepsilon(P) = 2\pi \int_{\{\theta \le \alpha_{[P]}\}} s(\theta) \sin \theta \, d\theta - P. \quad (A5)$$

As can be seen, $\varepsilon(P)$ never exceeds 0.01 and gets all the smaller that $K$ increases.

This then suggests that for our purposes, data error could be taken into account with the help of an equivalent angular Gaussian distribution in place of the more conventional Fisher distribution.

## Expanding the directional error into a 3-D error

The near-equivalence we just established between the Fisher distribution and an angular Gaussian distribution shows that a good approximation of the convolution (12) defining $p_i(\mathbf{u})$, could be obtained by computing the convolution between two angular Gaussian distributions. Unfortunately there is again no simple analytic form

<table><caption>**Table A1.** Departure $\varepsilon(P)$ (as defined by A5) of the probability $P$ predicted by the angular Gaussian distribution from that of the Fisher distribution it approximates, for various values of $P$ and of the concentration parameter $K$ (and corresponding $\alpha_{95}$ values).</caption><thead><tr><th>K</th><th>α<sub>{[P]}</sub></th><th>ε(0.1)</th><th>ε(0.2)</th><th>ε(0.3)</th><th>ε(0.4)</th><th>ε(0.5)</th><th>ε(0.6)</th><th>ε(0.7)</th><th>ε(0.8)</th><th>ε(0.9)</th></tr></thead><tbody><tr><td>25</td><td>28.3°</td><td>0.004</td><td>0.006</td><td>0.009</td><td>0.009</td><td>0.009</td><td>0.009</td><td>0.006</td><td>0.003</td><td>0.000</td></tr><tr><td>100</td><td>14°</td><td>0.001</td><td>0.001</td><td>0.002</td><td>0.002</td><td>0.002</td><td>0.001</td><td>0.001</td><td>0.000</td><td>0.000</td></tr><tr><td>225</td><td>9.4°</td><td>0.000</td><td>0.001</td><td>0.001</td><td>0.001</td><td>0.001</td><td>0.001</td><td>0.001</td><td>0.000</td><td>0.000</td></tr><tr><td>400</td><td>7.0°</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.001</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td></tr><tr><td>625</td><td>5.6°</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td></tr><tr><td>2500</td><td>2.8°</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td></tr></tbody></table>
for the result of such a convolution. However, this result suggests
that yet another approximation could usefully be made (in this ap-
pendix, we will drop the ‘i’ indexes for simplicity, since only a single
data **u**<sub>*i*</sub> is to be considered at a time).

To see this, first recall that **s**(**u**) is in fact the projection on Σ of
the 3-D Gaussian distribution **g**(**x**) predicted by the GGP process at
the site associated with the data to be tested. This then suggests that
we:

(i) seek a 3-D Gaussian distribution g_e(x) expending in 3-D the
directional error,

(ii) compute the convolution of g(**x**) with g_e(**x**) in R³, which
would then analytically lead to yet another 3-D Gaussian distribution
g_a(**x**) and,

(iii) analytically compute the angular Gaussian distribution
s_a(u) associated with this new 3-D Gaussian distribution (thanks
to eq. (6)).

Provided g_e(x) is properly chosen, we may then hope that the final
angular Gaussian distribution s_a(u) provides an accurate enough
approximation of p(u) as defined by the convolution (12).

In the limit case the GGP process predicts a g(**x**) distribution
characterized by E(**V**) = **m** and Cov(**V**, **V**) = 0, while the directional
error is characterized by the concentration parameter **K**, we already
know that choosing E(**V**_e) = 0 and Cov(**V**_e, **V**_e) = |**m**|²/**K** **I** for g_e(**x**)
would do the trick. In that case indeed, s_e(u) reduces to the single
direction **m**/|**m**| so that p(u) given by eq. (12) reduces to a Fisher
distribution centred on **m**/|**m**| with a concentration parameter **K**.
In the same time, g_a(x) would be characterized by E(V_a) = **m** and
Cov(V_a, V_a) = |**m**|²/**K** **I** which leads to s_a(u). This, we just saw, is very
close to a Fisher distribution centred on **m**/|**m**| with a concentration
parameter **K**.

In the more general case when, for g(x), Cov(V, V) is non-zero
but remains small enough, we can next anticipate that sticking
to the same 3-D expansion g_e(x) would also do the trick. Then
g_a(x) would be characterized by E(V_a) = m and Cov(V_a, V_a) =
Cov(V, V) + |m|²/K I, and hopefully, s_a(u) would again be close enough
to p(u).

<figure><img src="image_13.png" /><figcaption>**Figure A1.** Applying the same 3-D Gaussian error g<sub>e</sub>(x) (as symbolized by the two small-dash circles) to two different realizations V<sub>1</sub> and V<sub>2</sub> of a general 3-D Gaussian distribution g(x) (symbolized by the dashed ellipses), leads to different angular Gaussian errors on the unit sphere Σ if |V<sub>1</sub>| is larger than |V<sub>2</sub>| (as symbolized by the thick arcs on Σ).</figcaption></figure>

Note that this again can only be an approximation, and that in fact,
s_a(u) is not even be the exact result of the convolution of s(u) with
s_e(u) (the near-Fisherian angular Gaussian distribution associated
with g_e(x)). This is qualitatively illustrated in Fig. A1, which shows
that the same 3-D Gaussian error g_e(x) applied to two different
realizations V_1 and V_2 of g(x) with |V_1| larger than |V_2|, would
lead to an angular Gaussian error on Σ smaller for u_1 = V_1/|V_1|
than for u_2 = V_2/|V_2|. Fig. A1, however, also shows that provided
|m| is large enough compared to the dispersion implied by g(x) and
g_e(x), the distortion introduced should remain small. For the GGP
models we have tested, this fortunately happens to be the case most
of the time.