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

This needs to be taken into account when dealing with such non-Gaussian distributions


I N T RO D U C T I O N
Geomagnetic field modeling consists in converting large sets of data {γ i , s i } (where γ i denotes a datum, typically a field component, declination, inclination or intensity, measured at some location and time, and s i its assumed uncertainty) into a so-called model m = {m 1 , m 2 , . . . , m K } (where m j 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; 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 i is an adequate measure of the standard deviation of the residual expected for the datum γ i , standard maximum likelihood estimations can then be used to infer the model m, in which case, normalized residuals { γ i − γ i s i } ( γ i 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;Panovska et al. 2012Panovska et al. , 2015.  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. , 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.

M I X T U R E S O F R A N D O M D I S T R I B U T I O N S
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 α 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: 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 β 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 σ = 0.1, . . . , 0.9. In this case, the random variable β is assumed drawn with equal probability from a set of nine values (0.1, . . . , 0.9, the u(y) distribution is therefore assumed discrete), and the random variable α is next drawn from a Gaussian distribution with zero mean and standard deviation σ = y, where y value is the value taken by β.
. This trivially results in a mixed pdf of the form w(x) = 1 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 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, . . . , N} of several populations X k = {x i j | j = 1, 2, . . . , 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, . . . , N} can be considered as a sample population drawn from a mixture (or randomization) of unbiased Gaussian distributions with random variable α, and standard deviations defined by a random variable β independently drawn from its own distribution f β (y). In that case we may again write v(x, y) = 1 y √ 2π exp[− 1 2 ( x y ) 2 ], which indeed defines an unbiased Gaussian probability for the α 1038 A. Khokhlov and G. Hulot  variable once a value y has been drawn for the β variable, and make use of eq. (1) with u(y) = f β (y) to infer the pdf f ξ (x) of the resulting mixed distribution:

T Y P I C A L M I X T U R E S O F U N B I A S E D G AU S S I A N D I S T R I B U T I O N S
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 β (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.

Uniform mixtures
We begin with the case of uniform mixtures of unbiased Gaussian distributions with standard deviations varying between 0 and c, that is, In that case, f ξ (x) can be expressed in terms of an incomplete Gamma function (a, c) = ∞ c t a−1 exp(−t)dt after two variable substitutions in the integral expression of eq. (2). The first substitution is y = u −1 : leading to: 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: 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 ξ (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, In that case, as can readily be inferred from the previous case, the resulting distribution will take the form 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   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)).

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: This triangular mixture provides a typical example of a singlemode concave mixture, where Gaussian distributions involved in the mixture have a central maximum probability of having a given standard deviation σ = 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.

I N T E R P R E TAT I O N O F T H E N O N -G AU S S I A N B E H AV I O U R O F R E S I D UA L S I N G E O M A G N E T I S M
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.

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  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  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).

Marine magnetic anomalies
Another relevant example of non-Gaussian behaviour of magnetic residuals was provided by , 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  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  erroneously assigned a 3σ value to the standard deviation σ they report in their fig. 1a (539.8 nT); correcting for this error leads to a σ value very close to the σ = 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).  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 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. 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 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.

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 bestfit spline representation of individual sediment records and, among 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 α 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.

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 midlatitudes 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, 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.  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. 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 conditions 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 surprise 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 similar 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 synthetic histogram remarkably similar to the one observed (Fig. 12, compare right-hand plot to the left-hand plot). We next also computed 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 independently for the data from the Alpha and Bravo satellites, leading 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.

C O N C L U D I N G C O M M E N T S
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;, marine magnetic anomalies (e.g. , magnetic sediment data (e.g. Korte & Constable 2006;Panovska et al. 2012Panovska et al. , 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.

A C K N O W L E D G E M E N T S
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.