# On the cause of the non-Gaussian distribution of residuals in geomagnetism

A. Khokhlov, G. Hulot

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

## To cite this version:

A. Khokhlov, G. Hulot. On the cause of the non-Gaussian distribution of residuals in geomagnetism. Geophysical Journal International, 2017, 209, pp.1036-1047. 10.1093/gji/ggx071 . insu-03748861

HAL Id: insu-03748861

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

Submitted on 10 Aug 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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# On the cause of the non-Gaussian distribution of residuals in geomagnetism

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

Accepted 2017 February 22. Received 2017 February 15; in original form 2016 December 1

## SUMMARY

To describe errors in the data, Gaussian distributions naturally come to mind. In many practical instances, indeed, Gaussian distributions are appropriate. In the broad field of geomagnetism, however, it has repeatedly been noted that residuals between data and models often display much sharper distributions, sometimes better described by a Laplace distribution. In this study, we make the case that such non-Gaussian behaviours are very likely the result of what is known as mixture of distributions in the statistical literature. Mixtures arise as soon as the data do not follow a common distribution or are not properly normalized, the resulting global distribution being a mix of the various distributions followed by subsets of the data, or even individual datum. We provide examples of the way such mixtures can lead to distributions that are much sharper than Gaussian distributions and discuss the reasons why such mixtures are likely the cause of the non-Gaussian distributions observed in geomagnetism. We also show that when properly selecting subdata sets based on geophysical criteria, statistical mixture can sometimes be avoided and much more Gaussian behaviours recovered. We conclude with some general recommendations and point out that although statistical mixture always tends to sharpen the resulting distribution, it does not necessarily lead to a Laplacian distribution. This needs to be taken into account when dealing with such non-Gaussian distributions.

**Key words:** Archaeomagnetism; Magnetic anomalies: modelling and interpretation; Magnetic field variations through time; Palaeomagnetism; Satellite magnetics; Probability distributions.

## 1 INTRODUCTION

Geomagnetic field modeling consists in converting large sets of data {γᵢ, sᵢ} (where γᵢ denotes a datum, typically a field component, declination, inclination or intensity, measured at some location and time, and sᵢ its assumed uncertainty) into a so-called model **m** = {m<fnref n="1" />, m<fnref n="2" />, ..., m^K} (where m' denotes model parameters). Such models aim to provide the best mathematical description of the geomagnetic field accounting for the observed data when only a limited number of parameters are being used. Such approaches are typically used to recover spherical harmonic representations of the main field (e.g. Olsen *et al.* 2015), the lithospheric field (e.g. Thébault *et al.* 2016) and the ionospheric field (e.g. Chulliat *et al.* 2016) when using satellite data. They also are used to recover similar representations of the historical field when using historical data (e.g. Bloxham *et al.* 1989; Jackson *et al.* 2000; Walker & Jackson 2000) and of the archeomagnetic and time-averaged palaeomagnetic fields when using archeomagnetic, sediment and palaeomagnetic data (e.g. Licht *et al.* 2013; Johnson & Constable 1997). They are further commonly used to recover best representations of the temporal evolution of the

local field (e.g. Thébault & Gallet 2010; Panovska *et al.* 2012; Hellio *et al.* 2014). For such modeling efforts, it is important that appropriate statistical assumptions are being made for the *a priori* distribution of residuals (differences between the data and the model's predictions). The statistical properties of these residuals, however, are not always well characterized. In such circumstances, assuming that residuals follow a Gaussian distribution would make sense, since such a distribution often arises naturally as a consequence of the central limit theorem when errors act in an additive manner (see e.g. Feller 1971). Relying on this assumption, and provided that sᵢ is an adequate measure of the standard deviation of the residual expected for the datum γᵢ, standard maximum likelihood estimations can then be used to infer the model **m̃**, in which case, normalized residuals {γᵢ - ̂γᵢ / sᵢ} (γ̃ᵢ being the datum value predicted by the model **m̃**) would be expected to follow a standard normal distribution. Yet, this turns out to often not be the case, with residuals often displaying a marked trend to display a sharper distribution, sometimes much closer to that of a so-called Laplace distribution (e.g. Jackson *et al.* 2000; Walker & Jackson 2000; Panovska *et al.* 2012, 2015). Various strategies have been put

<footnote n="1">¹IEPT Russian Academy of Sciences, 84/32, Profsoyuznaya 117997 Moscow, Russia</footnote>

<footnote n="2">²Institute of Physics of the Earth, Russian Academy of Sciences, 10, B.Gruzinskaya 123242 Moscow, Russia</footnote>

<footnote n="3">³Institut de Physique du Globe de Paris, Sorbonne Paris Cité, Université Paris Diderot, CNRS, F-75005 Paris, France. E-mail: gh@ipgp.fr</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Result of an equally probable mixture of nine Gaussian distributions with zero means and standard deviations $\sigma = 0.1, \dots, 0.9$. Green solid line is for the resulting pdf $w(x) = \frac{1}{9} \sum_{i=1}^{9} \frac{10}{i\sqrt{2\pi}} \exp[-\frac{1}{2}(\frac{10x}{i})^2]$, dashed line is for a Laplace distribution best-fitting $w(x)$, namely $1.3 \exp(-2.6|x|)$. Left: linear scale and right: decimal semi-logarithmic scale.</figcaption></figure>

forward to take this into account, by acknowledging at the very onset of the modeling procedure that residuals follow a Laplace, rather than a Gaussian, distribution (e.g. Walker & Jackson 2000), by allowing the modeling procedure to empirically seek, rather than a priori assume, the distribution of residuals (e.g. Constable 1988), or by implementing an Iteratively Reweighted Least Squares (IRLS) method and Huber weights (Huber 1981) (see e.g. Olsen 2002; Thébault & Gallet 2010). No studies, however, have yet looked into the reason why such non-Gaussian behaviour arises in the geomagnetic context. The purpose of this paper is precisely to look into this and to show that the main cause can be traced back to the fact that residuals are usually assumed to sample a common distribution, whereas this may not be the case, and that normalized residuals may be incorrectly normalized, resulting in both cases in what is known in the statistical literature as mixture (or randomization) of distributions (Feller 1971; Barndorff-Nielsen et al. 1982). Indeed, as we will illustrate, such mixtures of otherwise normal distributions can produce non-Gaussian distributions with high kurtosis, often very close to Laplacian distributions.

## 2 MIXTURES OF RANDOM DISTRIBUTIONS

The simplest way to introduce the concept of mixtures of random distributions consists in starting from a family of probability densities (pdf) $v(x, y)$ where $x$ is the variable defining the values that a random variable $\alpha$ may take and $y$ is a parameter identifying each pdf of the family. A mixture of such pdfs can then be defined as a new pdf $w(x)$ using the formula:

$$w(x) = \int_{-\infty}^{\infty} v(x, y)u(y)dy \quad (1)$$

where $u(y)$ can be viewed as a weighing factor defining the contribution of each $v(x, y)$ pdf to the $w(x)$ pdf mixture. This mixture, however, can also be viewed as a randomization of the $v(x, y)$ pdfs, if the parameter $y$ is itself viewed as an independent random variable $\beta$ with pdf $u(y)$ (see e.g. Feller 1971).

Fig. 1 provides a first example of mixture of random distributions when considering nine Gaussian distributions with zero means and standard deviations $\sigma = 0.1, \dots, 0.9$. In this case, the random variable $\beta$ is assumed drawn with equal probability from a set of nine values (0.1, \dots, 0.9, the $u(y)$ distribution is therefore assumed discrete), and the random variable $\alpha$ is next drawn from a Gaussian distribution with zero mean and standard deviation $\sigma = y$, where $y$ value is the value taken by $\beta$. Here, therefore, $v(x, y) = \frac{1}{y\sqrt{2\pi}} \exp[-\frac{1}{2}(\frac{x}{y})^2]$. This trivially results in a mixed pdf of the form $w(x) = \frac{1}{9} \sum_{i=1}^{9} \frac{10}{i\sqrt{2\pi}} \exp[-\frac{1}{2}(\frac{10x}{i})^2]$. As can be seen, this mixed pdf is much sharper than that of an individual Gaussian distribution (Fig. 1, left). It in fact is much closer to a Laplace distribution of the form $\frac{a}{2} \exp(-a|x|)$, with $a = 2.6$ in the present case, as can be inferred from a best linear fit to the semi-logarithmic representation of this pdf (Fig. 1, right).

The above example is a simple one. More generally in geomagnetism, one can expect residuals to consist in a collection $\{x_n|n=1, 2, \dots, N\}$ of several populations $X_k = \{x_{i_j}|j=1, 2, \dots, n_k\}$ of residuals, each assumed to be drawn from a Gaussian distribution with some expectation $m_k$ and standard deviation $s_k$. Considering the impact of mixtures of Gaussian distributions with different expectation $m_k$, also modeled as drawn from a random distribution, is of course possible (see Barndorff-Nielsen et al. 1982). Here, however, we ignore this possibility, and assume $m_k = 0$ (as was already the case in the previous example). The reason for this is twofold. First, because such biases and their variability can be expected to be small compared to the variability in the standard deviations $s_k$, and second and foremost, because variability in the biases will usually result in some widening, and not sharpening, of the resulting mixed distribution (as intuition would tell, and as can readily be checked numerically).

In what follows, we thus focus on the consequences of the variability in the standard deviations $s_k$. Then, the total set of residuals $\{x_n|n=1, 2, \dots, N\}$ can be considered as a sample population drawn from a mixture (or randomization) of unbiased Gaussian distributions with random variable $\alpha$, and standard deviations defined by a random variable $\beta$ independently drawn from its own distribution $f_{\beta}(y)$. In that case we may again write $v(x, y) = \frac{1}{y\sqrt{2\pi}} \exp[-\frac{1}{2}(\frac{x}{y})^2]$, which indeed defines an unbiased Gaussian probability for the $\alpha$
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Plots of the probability density function $f_{\xi}(x) = \frac{1}{\sqrt{8\pi}} \Gamma(0, \frac{x^2}{2})$ when one considers uniform mixtures of unbiased Gaussian distributions with standard deviations varying between 0 and 1. Note the divergence of this pdf when $x$ goes to zero. Left: linear scale and right: decimal semi-logarithmic scale.</figcaption></figure>

variable once a value $y$ has been drawn for the $\beta$ variable, and make use of eq. (1) with $u(y) = f_{\beta}(y)$ to infer the pdf $f_{\xi}(x)$ of the resulting mixed distribution:

$$f_{\xi}(x) = \frac{1}{\sqrt{2\pi}} \int_{0}^{\infty} f_{\beta}(y) \exp\left(-\frac{1}{2} \frac{x^2}{y^2}\right) y^{-1} dy. \quad (2)$$

## 3 TYPICAL MIXTURES OF UNBIASED GAUSSIAN DISTRIBUTIONS

Mixtures of unbiased Gaussian distributions can occur in many forms, formally controlled by the distribution chosen for the standard deviation, that is, by the choice of $f_{\beta}(y)$ in eq. (2). Fig. 1 already provided an example of a discrete case of mixture of nine equally probable unbiased Gaussian. In this section, we consider other mixtures of interest.

### 3.1 Uniform mixtures

We begin with the case of uniform mixtures of unbiased Gaussian distributions with standard deviations varying between 0 and $c$, that is,

$$f_{\beta}(y) = \begin{cases} 1/c & 0 \le y \le c \\ 0 & y > c \end{cases} \quad (3)$$

In that case, $f_{\xi}(x)$ can be expressed in terms of an incomplete Gamma function $\Gamma(a, c) = \int_{c}^{\infty} t^{a-1} \exp(-t) dt$ after two variable substitutions in the integral expression of eq. (2). The first substitution is $y = u^{-1}$:

$$\begin{aligned} & \frac{1}{c\sqrt{2\pi}} \int_{0}^{c} \exp\left(-\frac{1}{2} \frac{x^2}{y^2}\right) y^{-1} dy \\ & = \frac{1}{c\sqrt{2\pi}} \int_{1/c}^{\infty} \exp\left(-\frac{1}{2}(xu)^2\right) u^{-1} du \end{aligned}$$

and the second is $t = \frac{1}{2}(xu)^2$:

$$\begin{aligned} & \frac{1}{c\sqrt{2\pi}} \int_{1/c}^{\infty} \exp\left(-\frac{1}{2}(xu)^2\right) u^{-1} du = \frac{1}{2c\sqrt{2\pi}} \int_{x^2/(2c^2)}^{\infty} t^{-1} \exp(-t) dt \\ & = \frac{1}{c\sqrt{8\pi}} \Gamma\left(0, \frac{x^2}{2c^2}\right) \end{aligned}$$

leading to:

$$f_{\xi}(x) = \frac{1}{c\sqrt{8\pi}} \Gamma\left(0, \frac{x^2}{2c^2}\right) \quad (4)$$

In particular, when $c=1$, that is, when one considers uniform mixtures of unbiased Gaussian distributions with standard deviations varying between 0 and 1, this leads to:

$$f_{\xi}(x) = \frac{1}{\sqrt{8\pi}} \Gamma\left(0, \frac{x^2}{2}\right). \quad (5)$$

The corresponding plots are shown in Fig. 2. As can be seen, this leads to a distribution which is extremely sharp, much sharper in its central part than an a Laplace distribution, since $f_{\xi}(x)$ actually goes to infinity when $x$ goes to zero.

Another case of interest is when one considers a uniform mixture of unbiased Gaussian distributions with standard deviations varying between $c-d$ and $c+d$, that is,

$$f_{\beta}(y) = \begin{cases} 1/2d & c-d \le y \le c+d \\ 0 & 0 \le y < c-d \text{ or } y > c+d \end{cases} \quad (6)$$

In that case, as can readily be inferred from the previous case, the resulting distribution will take the form

$$f_{\xi}(x) = \frac{1}{2d\sqrt{8\pi}} \left( \Gamma\left(0, \frac{x^2}{2(c+d)^2}\right) - \Gamma\left(0, \frac{x^2}{2(c-d)^2}\right) \right). \quad (7)$$

Fig. 3 shows such distributions when considering $c = 0.5$, and $d$ progressively increasing between 0.1 and 0.4 (i.e. for uniform mixtures of unbiased Gaussian distributions with standard deviations varying, respectively, between 0.4 and 0.6, 0.3 and 0.7, 0.2 and 0.8, 0.1 and 0.9). As can clearly be seen, as $d$ increases, the distribution progressively changes from a near Gaussian distribution (as one would expect when $d$ is small enough) to a distribution getting much closer to a Laplace distribution. This figure also makes it clear that
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Plots of the probability density functions $f_{\xi}(x) = \frac{1}{2d\sqrt{8\pi}}(\Gamma(0, \frac{x^2}{2(c+d)^2}) - \Gamma(0, \frac{x^2}{2(c-d)^2}))$ when one considers uniform mixtures of unbiased Gaussian distributions with standard deviations varying between $c-d$ and $c+d$, $c = 0.5$ and $d = 0.1, 0.2, 0.3$ and $0.4$ (resulting in a progressively sharpening mixed distribution when $d$ increases). Left: linear scale and right: semi-logarithmic scale.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 4.** Plots of the probability density function $f_{\xi}(x)$ when one considers a triangular mixture as defined by eq. (8). Left: linear scale and right: decimal semi-logarithmic scale.</figcaption></figure>

the sharp peak seen in the previous case (Fig. 2) is the consequence of including Gaussian distributions with vanishingly small standard deviations in the mixed distribution (Note, indeed, that assuming $c = 0.5$ and $d = 0.5$ in eq. (7), is equivalent to considering eq. (5)).

## 3.2 Other mixtures

Any other mixtures of interest can also easily be computed numerically by simply relying on a direct Monte Carlo approach. Consider, for instance, the case of a triangular mixture defined by:

$$f_{\beta}(y) = \begin{cases} 2(1 - |2y - 1|) & 0 \le y \le 1 \\ 0 & y > 1 \end{cases}. \qquad (8)$$

This triangular mixture provides a typical example of a single-mode concave mixture, where Gaussian distributions involved in the mixture have a central maximum probability of having a given standard deviation $\sigma = 0.5$ (i.e. the same central standard deviation as considered in all previous mixture examples), and a probability of having larger or smaller standard deviations decreasing as one moves away from this central value. As is illustrated in Fig. 4, this

is again enough to lead to a behaviour close to that of a Laplace distribution.

# 4 INTERPRETATION OF THE NON-GAUSSIAN BEHAVIOUR OF RESIDUALS IN GEOMAGNETISM

We now turn to the interpretation of a number of non-Gaussian behaviour of residuals identified by previous authors in the general context of geomagnetism.

## 4.1 Historical magnetic observations

We start by considering the case of historical magnetic observations. As has first been pointed out some decades ago (Bloxham et al. 1989), data used to compute spherical harmonic models of the historical main field usually lead to normalized residuals with distributions distinctly sharper than that of a Gaussian distribution, and somewhat closer to that of a Laplace distribution (see e.g. fig. 21 in
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Histogram of errors in 18 940 historical observations of compass declinations (in degrees) taken on board ship during the 17th and 18th centuries. The solid line represents the observations, the dashed line represents a best-fit Laplace distribution. The standard deviation for this distribution is $\sigma = 0.46^{\circ}$ (after fig. 7a of Jackson *et al.* 2000).</figcaption></figure>

Bloxham *et al.* 1989). This sharpening effect was later considered in more details by Jackson *et al.* (2000), who pointed out that a similar effect could also be seen in the distribution of errors in repeated historical observations of declination made on board ships during the 17th and 18th centuries (Fig. 5, reproduced from fig. 7a of Jackson *et al.* 2000). In that case, the errors plotted are not normalized residuals, but angular differences between individual declination measurements taken on a given day and the mean of these measurements. It is not unreasonable to assume that if all observations had been made by the same person under similar conditions, the distribution would have been much more Gaussian. But only few such measurements were usually made on a given day, and conditions of observations varied significantly from one day to the next. In addition, these observations were made on different ships by different observers using different instruments. Fig. 5 thus very likely shows a distribution resulting from a mixture of presumably unbiased Gaussian distributions, similar to those discussed in Section 3 above. Comparing Fig. 5 with Figs 2–4 would in fact suggest that the mixture involves errors with standard deviations varying by at least an order of magnitude, compatible with the fact, also reported by Jackson *et al.* (2000) that errors ‘with a magnitude perceived to be large enough to merit special mention in logs decreased through time, from a couple of degrees in the 17th century to a few minutes late in the 18th century’.

The effect discussed above is related, but not strictly identical to the one originally highlighted by Bloxham *et al.* (1989) (see also Walker & Jackson 2000) when considering normalized residuals with respect to the spherical harmonic models they computed for the historical main field. Such residuals not only result from measurement errors of the type just discussed, but also from errors linked to the fact that these data also contain signals produced by, for example, the ionosphere, the magnetosphere and most prominently, the crustal field. These non-modeled signals contribute to the error budget, and their contribution will therefore also affect the

form of the distribution of the final normalized residuals. Although such contributions can possibly be modeled as independent sources of (approximately unbiased) Gaussian noise, as we will later see, their magnitude depends on where and when (the latter in the case of signals of ionospheric and magnetospheric origin) the data have been acquired. This variability, and that of the measurement errors discussed above, is usually very partially, if ever, taken into account, particularly when modeling the historical field. As a result, if a common Gaussian error is assumed for each type of error and data, both in the modeling process and in the computation of the final normalized residuals (as was done by Bloxham *et al.* 1989), this will again amount to mix data errors originating from Gaussian distributions with varying standard deviations, and result in a sharp near-Laplacian distribution as observed in fig. 21 of Bloxham *et al.* (1989).

## 4.2 Marine magnetic anomalies

Another relevant example of non-Gaussian behaviour of magnetic residuals was provided by Walker & Jackson (2000), who considered marine magnetic anomalies collected over 30 yr by cruises over the world's oceans (held in the Scripps Institute of Oceanography database). As pointed out by these authors, this distribution is again very sharp and looks quite similar to that of a Laplace distribution (their fig. 1a). In that case, we could access a similar database (courtesy of J. Dyment and Y. Choi), which was recently used for the purpose of building the second version of the World Digital Magnetic Anomaly Map (Lesur *et al.* 2016) and an associated global equivalent magnetization map of the oceanic lithosphere (Dyment *et al.* 2015). Fig. 6 shows the resulting histogram. This histogram contains roughly 50 per cent more data than the one originally shown by Walker & Jackson (2000) and looks very similar when plotted in the same linear scale, as one would have expected (left in Fig. 6). Indeed, statistics derived from both distributions are very similar [Note, incidentally, that Walker & Jackson (2000) erroneously assigned a $3\sigma$ value to the standard deviation $\sigma$ they report in their fig. 1a (539.8 nT); correcting for this error leads to a $\sigma$ value very close to the $\sigma = 140.2$ nT we found for the updated distribution shown in Fig. 6]. This distribution is very clearly non-Gaussian and sharp. Plotting it in a semi-logarithmic scale reveals that it in fact is even sharper than a Laplace distribution (which would lead to a triangular shape when plotted in this way, recall Fig. 1). This, again, is likely the result of some complex mixture of distributions with a wide range of standard distributions.

Having access to the full data set and to additional useful ancillary information allowed us to derive evidence that this is indeed the case. Such marine magnetic anomalies are derived from field intensity measurements collected by towing scalar magnetometers at some distance behind ships to avoid non-natural magnetic signals. These intensity measurements are next corrected for the intensity predicted by global field models such as the International Geomagnetic Reference Field (IGRF) model (for details about IGRF models, see e.g. Thébault *et al.* 2015). The resulting so-called magnetic anomalies are then expected to mainly reflect the contribution of the magnetic signal produced by the magnetization of the oceanic crust lying below the ship. This magnetization is well known to be due to the fact that the crust forms at oceanic ridges where it acquires a (mainly) thermoremanent magnetization before moving away from these ridges as a result of seafloor spreading, progressively reaching greater depths and forming linear magnetized structures running parallel on both sides of the ridges (see e.g. Tivey 2007).
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Histogram of the 15 838 070 marine magnetic anomalies recently used to compute the second version of the World Digital Magnetic Anomaly Map (Lesur *et al. 2016*). Left: linear scale and right: decimal semi-logarithmic scale. Anomalies are provided in nanotesla.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 7.** World Digital Magnetic Anomaly Map over the oceans (Lesur *et al. 2016*), adapted from fig. 1a of Dyment *et al. (2015)*. The grey-shaded area southwest of Australia is the one from which the anomalies used to build the histograms of Figs 8 and 9 come.</figcaption></figure>

Magnetization having been acquired at the time of the very initial cooling phase of the oceanic crust, these linear magnetized structures reflect both the seafloor expansion (and subsequent tectonics) and the history of the main magnetic field, which often, but irregularly, reversed in the geological past (see e.g. Hulot *et al. 2010*). Marine magnetic anomalies will therefore vary in magnitude depending on many factors, most importantly the depth of the oceanic basement (recall, indeed, that magnetic anomalies are computed at sea level, and are therefore weaker if sources are further away), the latitude of the region where these anomalies are observed (for similar depths and ages, magnetic anomalies are stronger at high latitudes than near the equator because of the dipolar structure of the main field), as

well as the orientation of the magnetized structures (for details about these subtleties, see e.g. Lesur *et al. 2016*, and references therein). These characteristics can clearly be seen in the WDMAM over the oceans as shown in fig. 1a of Dyment *et al. (2015)*, which we reproduce in Fig. 7. When focusing on specific regions with limited latitude variations, similar orientation of the magnetized structures and a limited range of basement depths, one may thus hope to capture the fundamental distribution of the magnetization responsible for the global distribution of the marine magnetic anomalies shown in Fig. 6. Following this line of reasoning, we focused on the region southwest of Australia, delimited by latitudes varying between 40° and 55° south and longitudes varying between 95°
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Histograms of the marine magnetic anomalies used to compute the second version of the WDMAM when only considering data from the region south of Australia delimited by latitudes between 40° and 55° south and longitudes between by 95° and 120° east, and depth to basement ranging between 3 and 3.5 km (top) and between 3.5 and 4 km (bottom). Left: true data and right: simulated unbiased Gaussian synthetic data with identical standard deviation; all plots in decimal semi-logarithmic scales.</figcaption></figure>

and 120° east (grey-shaded area in Fig. 7). In this region, basement depths vary between a little less than 3 and 5 km for the data available and anomalies are running along a common roughly east to west orientation. Further separating this regional data set into data sets corresponding to 500 m ranges of depth varia-

tions, and considering the most numerous of these subdata sets leads to histograms that are no longer as sharp as the global histogram shown in Fig. 6, and much closer to that of Gaussian distributions. This can be seen in Figs 8 and 9, which show the corresponding histograms together with synthetic histograms
<figure><img src="image_12.png" /><figcaption>**Figure 9.** Same as Fig. 8, but for depth to basement ranging between 4 and 4.5 km (top) and between 4.5 and 5 km (bottom).</figcaption></figure>

produced numerically from an identical number of data drawn from a Gaussian distribution with the same standard deviation as that of the real data (all plotted in semi-logarithmic scales). Similar plots can be found when considering other regions of the world and carefully selecting subregions in the same way. It thus clearly appears that the very sharp distribution seen in the global histogram of marine magnetic anomalies shown in Fig. 6 is indeed the result of some complex mixture of near-Gaussian distributions with a wide range of standard distributions.

## 4.3 Archeomagnetic and sediment data

We now switch to yet another set of examples of non-Gaussian behaviour of residuals, this time found when considering sediment and archeomagnetic data commonly used for main field modeling over the past few millennia and during the Holocene.

We start by considering the spline analysis of Holocene sediment magnetic records recently carried out by Panovska *et al.* (2012). In this study, the authors used a robust spline analysis to find a best-fit spline representation of individual sediment records and, among
<figure><img src="image_13.png" /><figcaption>**Figure 10.** Left: example of a Holocene declination sediment record, showing the raw data (green diamonds) and best-fit spline (solid blue line) and right: histogram of the corresponding residuals (grey boxes) together with a best-fit Laplacian distribution (solid red line); units are in degrees (reproduced from fig. 4 of Panovska et al. 2012).</figcaption></figure>

other things, infer an estimate of the uncertainties affecting each of these records. They found that these uncertainties varied quite significantly from one record to the next, but also noted that for each record, the distribution of residuals with respect to the best-fit spline was usually much better accounted for by a Laplace distribution than by a Gaussian distribution (see figs 4 and 5 in Panovska et al. 2012). However, simple inspection of the misfit of the corresponding declination, inclination or relative palaeointensities series to their respective best-fit splines also clearly points at these misfits being much larger for some periods of time than for others. This can be seen in Fig. 10 (reproduced from fig. 4 of Panovska et al. 2012), where the fit of the spline to the declination data is much better around 1000 AD than around 1000 BC. This is a strong indication that uncertainties not only vary from one record to the next, as noted by Panovska et al. (2012), but also within a given record. These variations inevitably lead to some mixture of distributions when considering the entire set of residuals, and are the likely cause of the close-to-Laplacian distribution shown in Fig. 10 (right).

Using such sediment data to build time-varying spherical harmonic models of the main field also leads to normalized residuals that follow close-to-Laplacian distributions. This can be seen in figures such as fig. 2 of Panovska et al. (2015), who used the sediment data studied by Panovska et al. (2012) and the uncertainties these authors had derived. These sediment data represent 85 per cent of the data used to produce the models (the rest being archeomagnetic data, see below) and therefore contribute most to the distribution of normalized residuals plotted by these authors. Normalizing the residuals by the uncertainty estimates can be expected to rescale residuals with respect to each other. In the present case, however, uncertainties have been assumed uniform throughout each sediment record and renormalization cannot account for the variability of uncertainties within each record. Even renormalized, residuals from sediment data can thus again lead to a near-Laplacian distribution, as is indeed found in fig. 2 of Panovska et al. (2015).

Distributions of residuals with respect to spherical harmonic models are not always as strongly Laplacian as those just discussed. They can sometimes be only slightly sharper than Gaussian distributions. This is typically the case when archeomagnetic data are considered. A nice example is provided in fig. 5 of Korte & Constable (2006), where histograms of residuals of the intensity data used in computing the CALS7K.2 model have been plotted separately for the archeomagnetic intensity data and for the sedimentary intensity data (note that in this example, residuals are not normalized to the assumed uncertainty). Whereas the sediment data again lead to a fairly sharp distribution, the archeointensity data lead to a distribution only mildly sharper than that of a Gaussian. This would suggest that much less mixture of distributions occurs in the case

of archeointensity data than in the case of sediment data, and that the uncertainty with which the main field intensity can be recovered from archeological samples only modestly varies from one sample to the next, despite the complexity of the causes of these uncertainties (see e.g. Genevey et al. 2008; Suttie et al. 2011). Indeed, comparison of fig. 5 (left) of Korte & Constable (2006) with Fig. 3 would suggest that, for most of the data, relative uncertainties do not vary by much more than a factor 2. Similarly, mildly non-Gaussian behaviour of archeomagnetic declination and inclination residuals would indicate a fairly homogeneous data set. Note, however, that if uncertainties are erroneously assigned (such as when ignoring that converting an MAD angle into an $\alpha_{95}$ angle requires a factor close to 3, see e.g. Khokhlov & Hulot 2016), computing normalized residuals rather than raw residuals, can again lead to significantly more Laplacian mixtures of distributions.

## 4.4 Satellite data

We finally turn to the case of contemporary satellite data. To illustrate this case, we rely on data used to produce a model proposed as a candidate model for IGRF 2015 and which is fairly typical of models produced from satellite data (Vigneron et al. 2015). More specifically, we focus on that fraction of the data set which consists of absolute scalar data acquired by two of the Swarm satellites (Satellites Alpha and Bravo) at quasi-dipole latitudes ranging between -55° and +55°, and compute residuals with respect to the so-called VFM model of Vigneron et al. (2015). These scalar data cover a little less than a year (between 2013 November 29 and 2014 September 25) and were further selected following a number of criteria, among which magnetically quiet and night time conditions, to ensure that as little as possible non-modeled external signal is included in the data. This resulted in 42 160 data for the Alpha satellite and 42 175 for the Bravo satellite. These data can be expected to reflect the signal of the field of internal origin the model aims at modeling, any other source of signal being treated as a source of noise acting on top the very low instrumental and satellite noise (less than 0.3 nT, see Léger et al. 2015; Olsen et al. 2015; Fratter et al. 2016). Residuals thus mainly reflect the noise produced by whatever external field the model fails to capture. At the mid latitudes within which these data were acquired, these residuals are typically due to signals related to the ring-current and mid-latitude ionospheric currents, which only produce modest signals, given the quiet night time data selected. Indeed, computing daily standard deviations for the corresponding residuals for each of the two Alpha and Bravo satellites shows that these residuals are quite small, with standard deviations typically ranging between less than 1 nT and a few nanotesla, with a maximum value of a little less than 8 nT,
1-Day std for Alpha and Bravo non polar residuals

<figure><img src="image_14.png" /><figcaption>**Figure 11.** Standard deviations (in nT) computed every day for the mid-latitude residuals of the Swarm scalar data used to compute the VFM model of Vigneron *et al.* (2015). Blue large dots: data from the Swarm Alpha satellite and red dots: data from the Swarm Bravo satellite. Days are counted in Julian days, with 2000 January 1 taken as the reference.</figcaption></figure>

N=42175 scalar residual

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

N=42175 mixture of
simulated scalar residuals

<figure><img src="image_16.png" /><figcaption>**Figure 12.** Left: histogram of the mid-latitude residuals of the Swarm Bravo scalar data used to compute the VFM model of Vigneron *et al.* (2015) and right: histogram of an identical number of simulated data drawn from unbiased Gaussian distributions using the daily standard deviations shown in Fig. 11; all plots in decimal semi-logarithmic scales.</figcaption></figure>

as shown in Fig. 11. This figure, however, also makes it clear that
residuals do vary from one day to the next in a fairly consistent way
for both satellites, indicating that even though quiet magnetic con-
ditions have been selected, the noise produced by the external field
does vary on a daily basis by some factor of order 4. It thus is no sur-
prise that when combining residuals from all days for each of the two
satellites, the resulting histograms reveal a distribution very similar
to the type of mixed distribution predicted in Section 3. Fig. 12
(left) shows the corresponding distribution for satellite Bravo (a
very similar distribution is found for satellite Alpha). As expected,
and as had been observed by field modelers when inspecting simi-
lar satellite data residuals (see e.g. Olsen 2002), this distribution is
quite far from being Gaussian and is very close to Laplacian.

Having the data set at our disposal, we were also able to actually
confirm the fact that this close-to-Laplacian distribution indeed is

the result of a statistical mixture of otherwise essentially unbiased
normal distributions. This was tested in two ways. We first checked
that the distribution observed in Fig. 12 (left) could be reproduced
with the help of a mixture of unbiased normal distributions, using
exactly the same amount of synthetic data with standard deviations
changing every day in the same way as the real residuals (i.e. using
the standard deviations plotted in Fig. 11). This resulted in a syn-
thetic histogram remarkably similar to the one observed (Fig. 12,
compare right-hand plot to the left-hand plot). We next also com-
puted the histogram of the real residuals normalized to the daily
varying standard deviation as defined by Fig. 11 (see Fig. 13, left),
and compared it to the histogram of exactly the same amount of data
produced by a pure unbiased normal distribution (Fig. 13, right).
Again, the match is quite remarkable. Both tests were done inde-
pendently for the data from the Alpha and Bravo satellites, leading
<figure><img src="image_17.png" /><figcaption>**Figure 13.** Left: histogram of the normalized mid-latitude residuals of the Swarm Bravo scalar data used to compute the VFM model of Vigneron *et al.* (2015) and right: histogram of an identical number of simulated data drawn from an unbiased normal Gaussian distributions; all plots in decimal semi-logarithmic scales.</figcaption></figure>

to virtually undistinguishable histograms (Fig. 13 only shows the result for the Bravo satellite). These results thus clearly confirm that for satellite data also, the occurrence of strongly non-Gaussian distributions of residuals can be the result of statistical mixture of otherwise unbiased Gaussian distributions. Importantly, these results again highlight the fact that provided one has some geophysical insight into what may be the cause of the variability of the standard deviations of the residuals (in the present case, the day to day variability of the external field activity), the raw Gaussian nature of these residuals can be recovered.

## 5 CONCLUDING COMMENTS

Gaussian distributions naturally come to mind when it comes to describe errors in the data. The usual justification for such a choice is the well-known central limit theorem, which states that such distributions should arise when errors affecting the data can be expected to act in an additive manner (see e.g. Feller 1971). In many practical instances, indeed, Gaussian distributions properly describe data errors. This, however, is not always the case. In the broad field of geomagnetism, residuals between data and models aiming at providing the best description of these data often display much sharper distributions, sometimes much better described by a Laplace distribution. This has been found to be the case when considering historical magnetic data (e.g. Bloxham *et al.* 1989; Jackson *et al.* 2000; Walker & Jackson 2000), marine magnetic anomalies (e.g. Walker & Jackson 2000), magnetic sediment data (e.g. Korte & Constable 2006; Panovska *et al.* 2012, 2015) and modern satellite data (e.g. Olsen 2002; Olsen *et al.* 2015). In this study, we made the case that such non-Gaussian behaviours are very likely the result of what is known as mixture of distributions in the statistical literature (e.g. Barndorff-Nielsen *et al.* 1982). Such mixtures arise as soon as

the data do not follow a common distribution, the resulting global distribution being then a mixture of the various distributions followed by subsets of the data or even individual datum. We provided theoretical examples of the way such mixtures can lead to distributions that are much sharper than Gaussian distributions (Section 3). We also provided explicit reasons to believe that such mixtures are the underlying cause of the close-to-Laplacian distribution observed when considering historical magnetic data (Section 4.1) and sediment data (Section 4.3). We finally provided direct evidence that statistical mixture is also the major underlying cause of the very sharp distribution of marine magnetic anomalies (Section 4.2) and of the almost as sharp distribution of mid-latitude satellite scalar residuals (Section 4.4). In both these instances, we were further able to explicitly show that when properly selecting subdata sets based on geophysical criteria, much more Gaussian behaviours could be recovered, thereby proving that these distributions indeed result from such statistical mixtures.

Several conclusions of this study are well worth being highlighted. First is the general conclusion that because of the wide range of noise levels (both natural and instrumental) that can affect the data, statistical mixtures are very likely to occur in the context of geomagnetism. Second is the conclusion that relying on erroneous estimates of the *a priori* or *a posteriori* standard deviations (such as when wrongly assuming a common standard deviation for data of uneven quality) can also lead to some statistical mixture of the distribution of the (then wrongly) normalized residuals, even when errors affecting individual data can demonstrably be assumed Gaussian. Finally, and most importantly, is the conclusion that mixtures do not always lead to exact Laplacian distributions. Sometimes, when the underlying distributions involved in the mixture have standard deviations varying by only a small factor (say within a factor 3), as seems to be the case for the archeointensity data discussed in Section 4.3, the distribution will only be slightly
sharper than a Gaussian. In other instances, when the underlying distributions display much stronger variations, as is the case for the marine magnetic anomalies discussed in Section 4.2 and for the high-precision satellite scalar data discussed in Section 4.4, the resulting distribution can be even sharper than a Laplacian. This, together with the fact that other effects may well also contribute to make the distribution non-Gaussian, should serve as a warning that *a priori* assuming a Laplace distribution for data errors may not always be a better alternative to *a priori* assuming a Gaussian distribution. Rather, and whenever possible, one should try to avoid having to deal with statistical mixtures by identifying the geophysical (or other) causes of the occurrence of such mixtures and by recovering adequate measures of the uncertainties affecting each subset of data, in the way we did when considering marine magnetic anomalies and satellite scalar data. Alternatively, and whenever such an approach appears to be intractable (as seems to be the case when considering the sediment data, recall Section 4.3), adaptive approaches, such as the one proposed by Constable (1988) or the IRLS method combined with Huber weights (see e.g. Olsen 2002), are those to be preferred.

## ACKNOWLEDGEMENTS

The authors wish to thank an anonymous reviewer for helpful comments, N.M. Shapiro for discussion and suggestions on an initial draft of this manuscript, J. Dyment and Y. Choi for providing the marine magnetic anomaly data used in Section 4.2 and P. Vigneron for assistance in the analysis of the satellite data discussed in Section 4.4. This work was carried out while AK benefited from an invited Professor position within IPGP. It was partly funded by grant no. 14.Z50.31.0017 of the Russian Ministry of Science and Education, RFBR 15-05-01842, and by the Centre National d'Etudes Spatiales through the 'Travaux préparatoires et exploitation de la mission Swarm' project. This is IPGP contribution no. 3823.

## REFERENCES

* Barndorff-Nielsen, O., Kent, J. & Sørensen, M., 1982. Normal variance-mean mixtures and z distributions, *Int. Stat. Rev.*, **50**, 145–159.

* Bloxham, J., Gubbins, D. & Jackson, A., 1989. Geomagnetic secular variation, *Phil. Trans. R. Soc. Lond. A*, **329**, 415–502.

* Chulliat, A., Vigneron, P. & Hulot, G., 2016. First results from the Swarm Dedicated Ionospheric Field Inversion chain, *Earth Planets Space*, **68**, 104, doi:10.1186/s40623-016-0481-6.

* Constable, C.G., 1988. Parameter estimation in non-Gaussian noise, *Geophys. J. Int.*, **94**(1), 131–142.

* Dyment, J., Choi, Y., Hamoudi, M., Lesur, V. & Thebault, E., 2015. Global equivalent magnetization of the oceanic lithosphere, *Earth planet. Sci. Lett.*, **430**, 54–65.

* Feller, W., 1971. *An Introduction to Probability Theory and its Applications*, 3rd edn, Vol. 2, Wiley.

* Fratter, I., Léger, J.M., Bertrand, F., Jager, T., Hulot, G., Brocco, L. & Vigneron, P., 2016. Swarm Absolute Scalar Magnetometers first in-orbit results, *Acta Astronaut.*, **121**, 76–87.

* Genevey, A., Gallet, Y., Constable, C.G., Korte, M. & Hulot, G., 2008. ArcheoInt: an upgraded compilation of geomagnetic field intensity data for the past ten millennia and its application to the recovery of the past dipole moment, *Geochem. Geophys. Geosyst.*, **9**, Q04038, doi:10.1029/2007GC001881.

* Hellio, G., Gillet, N., Bouligand, C. & Jault, D., 2014. Stochastic modelling of regional archaeomagnetic series, *Geophys. J. Int.*, **199**(2), 931–943.

* Huber, P.J., 1981. *Robust Statistics*, Wiley, 308 pp.

* Hulot, G., Finlay, C.C., Constable, C.G., Olsen, N. & Mandea, M., 2010. The magnetic field of planet Earth, *Space Sci. Rev.*, **152**, 159–222.

* 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.L. & Constable, C.G., 1997. The time-averaged geomagnetic field: global and regional biases for 0–5 Ma, *Geophys. J. Int.*, **131**, 643–666.

* Khokhlov, A. & Hulot, G., 2016. Principal component analysis of paleomagnetic directions: converting a maximum angular deviation (MAD) into an α₉₅ angle, *Geophys. J. Int.*, **204**, 274–291.

* Korte, M. & Constable, C.G., 2006. On the use of calibrated relative paleointensity records to improve millennial scale geomagnetic field models, *Geochem. Geophys. Geosyst.*, **7**, Q09004, doi:10.1029/2006GC001368.

* Léger, J.M. et al., 2015. In-flight performances of the absolute scalar magnetometer vector mode on board the Swarm satellites, *Earth Planets Space*, **67**, 57, doi:10.1186/s40623-015-0231-1.

* Lesur, V., Hamoudi, M., Choi, Y., Dyment, J. & Thébault, E., 2016. Building the second version of the World Digital Magnetic Anomaly Map (WD-MAM), *Earth Planets Space*, **68**, 27, doi:10.1186/s40623-016-0404-6.

* Licht, A., Hulot, G., Gallet, Y. & Thébault, E., 2013. Ensembles of low degree archeomagnetic field models for the past three millennia, *Phys. Earth planet. Inter.*, **224**, 38–67.

* Olsen, N., 2002. A model of the geomagnetic field and its secular variation for epoch 2000 estimated from Ørsted data, *Geophys. J. Int.*, **149**(2), 454–462.

* Olsen, N. et al., 2015. The Swarm Initial Field Model for the 2014 geomagnetic field, *Geophys. Res. Lett.*, **42**, 1092–1098.

* Panovska, S., Finlay, C.C., Donadini, F. & Hirt, A., 2012. Spline analysis of Holocene sediment magnetic records: uncertainty estimates for field modeling, *J. geophys. Res.*, **117**, B02101, doi:10.1029/2011JB008813.

* Panovska, S., Korte, M., Finlay, C.C. & Constable, C.G., 2015. Limitations in paleomagnetic data and modelling techniques and their impact on Holocene geomagnetic field models, *Geophys. J. Int.*, **202**(1), 402–418.

* Suttie, N., Holme, R., Hill, M.J. & Shaw, J., 2011. Consistent treatment of errors in archaeointensity implies rapid decay of the dipole prior to 1840, *Earth planet. Sci. Lett.*, **304**, 13–21.

* Thébault, E. & Gallet, Y., 2010. A bootstrap algorithm for deriving the archeomagnetic field intensity variation curve in the Middle East over the past 4 millennia BC, *Geophys. Res. Lett.*, **37**, L22303, doi:10.1029/2010GL044788.

* Thébault, E. et al., 2015. International Geomagnetic Reference Field: the 12th generation, *Earth Planets Space*, **67**, 79, doi:10.1186/s40623-015-0228-9.

* Thébault, E., Vigneron, P., Langlais, B. & Hulot, G., 2016. A Swarm lithospheric magnetic field model to SH degree 80, *Earth Planets Space*, **68**, 126, doi:10.1186/s40623-016-0510-5.

* Tivey, M.A., 2007. Vine-Matthews-Morley hypothesis, in *Encyclopedia of Geomagnetism and Paleomagnetism*, pp. 980–984, eds Gubbins, D. & Herrero-Bervera, E., Springer.

* Vigneron, P. et al., 2015. A 2015 International Geomagnetic Reference Field (IGRF) candidate model based on swarm's experimental absolute magnetometer vector mode data, *Earth Planets Space*, **67**, 95, doi:10.1186/s40623-015-0265-4.

* Walker, M.R. & Jackson, A., 2000. Robust modelling of the Earth's magnetic field, *Geophys. J. Int.*, **143**(3), 799–808.