Merge branch 'main' of https://github.com/PmagPy/Essentials-JupyterBook

Author: Swanson-Hysell

book/chapters/chapter12.md

@@ -12,7 +12,7 @@ numbering:
 Read {cite:t}`efron1993`; {cite:t}`tauxe1991`.
 :::
 
-Paleomagnetists have depended since the 1950's on the special statistical framework developed by {cite:t}`fisher1953` for the analysis of unit vector data. The power and flexibility of a variety of tools based on Fisher statistics enables quantification of parameters such as the degree of rotation of a crustal block, or whether the geomagnetic field really averages to a geocentric axial dipole independent of polarity. These tools, however, require that the paleomagnetic data belong to a particular parametric distribution -- the Fisher distribution.
+Paleomagnetists have depended since the 1950s on the special statistical framework developed by {cite:t}`fisher1953` for the analysis of unit vector data. The power and flexibility of a variety of tools based on Fisher statistics enables quantification of parameters such as the degree of rotation of a crustal block, or whether the geomagnetic field really averages to a geocentric axial dipole independent of polarity. These tools, however, require that the paleomagnetic data belong to a particular parametric distribution -- the Fisher distribution.
 
 :::{figure} ../figures/chapter12/vgp-di.png
 :name: fig:vgp-di
@@ -26,7 +26,7 @@ In many important situations, the Fisher distribution fails to represent paleoma
 
 Another example of the inadequacy of the Fisher distribution is the fact that the magnetic field exists in two stable polarity states. Because the Fisher distribution allows only uni-modal data, bi-polar data must be separated into separate modes or one mode must be "flipped" to the antipode prior to calculating a mean. Remanence vectors composed of several components tend to form streaked distributions. Structural complications (e.g., folding) can lead to streaked distributions of directional data. And, inclination error arising from flattening of directions tends to form "squashed" directional distributions that are wider in the horizontal plane than in the vertical. These are all commonly observed pathologies in directional data sets that lead to non-Fisherian data distributions.
 
-Thus, non-Fisherian data are a fact of paleomagnetic life. The Fisher-based tests can frequently be inappropriate and could result in flawed interpretations. In [Chapter 11](#chap:fisher) we learned the basics of Fisher statistics and how to test data sets against a Fisher distribution. In this chapter, we will discuss what to do when Fisher statistics fail. We will begin with parametric approaches that treat certain types of non-Fisherian data. We then turn to the use of non-parametric methods such as the bootstrap and jackknife in paleomagnetic applications.
+Thus, non-Fisherian data are a fact of paleomagnetic life. The Fisher-based tests can frequently be inappropriate and could result in flawed interpretations. In [Chapter 11](#chap:fisher), we learned the basics of Fisher statistics and how to test data sets against a Fisher distribution. In this chapter, we will discuss what to do when Fisher statistics fail. We will begin with parametric approaches that treat certain types of non-Fisherian data. We then turn to non-parametric methods for paleomagnetic applications, focusing on the bootstrap.
 
 (sect:nonfisherian-parametric)=
 ## Non-Fisherian Parametric Approaches
@@ -40,7 +40,7 @@ $$
 F = c(\kappa,\beta)^{-1} \exp(\kappa \cos \alpha + \beta \sin^2 \alpha \cos 2\phi)
 $$ (eq:kent)
 
-The mean direction in a Kent distribution is estimated in the same way as for the Fisher distribution (see [Chapter 11](#chap:fisher)). The difference is that when transformed to the mean direction, Kent declinations are not uniformly distributed around the mean. If we calculate eigenparameters for the orientation matrix of the data (see [](#app:eigen)), then the major and minor eigenvectors ($\V_2, \V_3$) lie in a plane orthogonal to the mean direction along the axis with the most and least scatter respectively. In [Equation %s](#eq:kent), $\alpha$ is the angle between a given direction and the true mean direction, and $\phi$ is the angle in the $\V_2, \V_3$ plane with $\phi$ = 0 parallel to $\V_2$. $\kappa$ is a concentration parameter similar to the Fisher $\kappa$ and $\beta$ is the "ovalness" parameter. $c(\kappa,\beta)$ is a complicated function of $\kappa$ and $\beta$. When $\beta$ is zero, the Kent distribution reduces to a Fisher distribution. Details of the calculation of Kent 95% confidence ellipses are given in [](#app:kent).
+The mean direction in a Kent distribution is estimated in the same way as for the Fisher distribution (see [Chapter 11](#chap:fisher)). The difference is that when transformed  into coordinates centered on the mean direction, Kent declinations are not uniformly distributed around the mean. If we calculate eigenparameters for the orientation matrix of the data (see [](#app:eigen)), then the major and minor eigenvectors ($\V_2, \V_3$) lie in a plane orthogonal to the mean direction along the axis with the most and least scatter respectively. In [Equation %s](#eq:kent), $\alpha$ is the angle between a given direction and the true mean direction, and $\phi$ is the angle in the $\V_2, \V_3$ plane with $\phi$ = 0 parallel to $\V_2$. $\kappa$ is a concentration parameter similar to the Fisher $\kappa$ and $\beta$ is the "ovalness" parameter. $c(\kappa,\beta)$ is a complicated function of $\kappa$ and $\beta$. When $\beta$ is zero, the Kent distribution reduces to a Fisher distribution. Details of the calculation of Kent 95% confidence ellipses are given in [](#app:kent).
 
 If we were to collect data from the equatorial region, we might well obtain a set of directions such as those shown in [](#fig:confidence)a. [Note that the center of the diagram is the expected direction -- not down as is more common.] The Fisher $\alpha_{95}$ circle of confidence for this data set is shown in [](#fig:confidence)a. The Kent ellipse ([](#fig:confidence)b) clearly represents the distribution of data better than the Fisher $\alpha_{95}$, being elongate in the same sense as the data themselves.
 
@@ -61,27 +61,34 @@ $$
 F = \frac{1}{4\pi d(k_1,k_2)} \exp( k_1\cos^2 \phi+k_2\sin^2\phi)\sin^2 \alpha
 $$
 
-where $\alpha$ and $\phi$ are as in the Kent distribution, $k_1,k_2$ are concentration parameters ($k_1<k_2<0$) and $d(k_1,k_2)$ is a constant of normalization. Values for $k_1,k_2$ can be estimated by numerical integration and can be converted into 95% confidence ellipses, the details of which are given in [](#app:bing). In a nut shell, the $\V_1$ eigenvector of the orientation matrix (associated with the largest eigenvalue, see [](#app:eigen)) is the principal direction and the semi-axes of the 95% confidence ellipse are proportional to the intermediate and minimum eigenvalues. The Bingham principal direction therefore is not necessarily the same as the Fisher or Kent mean. If we take each vector end-point to be a mass, the Bingham principal direction is the axis about which the moment of inertia of the masses would be least. The Fisher mean is somewhat different, in that it is the vector sum of unit vectors. The Bingham mean (principal direction) is less affected by outliers than the Fisher mean, lying closer to the center of mass of data points.
+where $\alpha$ and $\phi$ are as in the Kent distribution, $k_1,k_2$ are concentration parameters ($k_1<k_2<0$) and $d(k_1,k_2)$ is a constant of normalization. Values for $k_1,k_2$ can be estimated by numerical integration and can be converted into 95% confidence ellipses, the details of which are given in [](#app:bing). In a nutshell, the $\V_1$ eigenvector of the orientation matrix (associated with the largest eigenvalue, see [](#app:eigen)) is the principal direction and the semi-axes of the 95% confidence ellipse are proportional to the intermediate and minimum eigenvalues. The Bingham principal direction therefore is not necessarily the same as the Fisher or Kent mean. If we take each vector end-point to be a mass, the Bingham principal direction is the axis about which the moment of inertia of the masses would be least. The Fisher mean is somewhat different, in that it is the vector sum of unit vectors. The Bingham mean (principal direction) is less affected by outliers than the Fisher mean, lying closer to the center of mass of data points.
 
-The principle drawback of the Bingham distribution is that because the orientation matrix uses the entire data set (normal and reverse) the two modes are assumed to be antipodal and to share the same distribution parameters. The question of whether normal and reverse data sets are antipodal and have the same dispersion is in fact one we may wish to ask! One could separate the two modes prior to calculation of the Bingham ellipse, but then the rationale for using the Bingham distribution is lost. Also, many published descriptions of the Bingham calculation (e.g., {cite:t}`onstott1980`; {cite:t}`borradaile2003`) have errors in them. The source code for calculating Bingham statistics in widely used paleomagnetic data reduction programs is generally not available, and it is unknown whether these programs contain bugs.
+The principal drawback of the Bingham distribution is that because the orientation matrix uses the entire data set (normal and reverse) the two modes are assumed to be antipodal and to share the same distribution parameters. The question of whether normal and reverse data sets are antipodal and have the same dispersion is in fact one we may wish to ask! One could separate the two modes prior to calculation of the Bingham ellipse, but then the rationale for using the Bingham distribution is lost. Also, many published descriptions of the Bingham calculation (e.g., {cite:t}`onstott1980`; {cite:t}`borradaile2003`) have errors in them. The source code for calculating Bingham statistics in widely used paleomagnetic data reduction programs is generally not available, and it is unknown whether these programs contain bugs.
 
 (sect:bingham-legoff)=
 ### The Bingham-LeGoff Approximation
 
 Estimating the parameters for the Bingham ellipse exactly is computationally taxing and all of the available "canned" programs use the look up table of {cite:t}`mardia1977` (see [](#app:bing)). {cite:t}`legoff1992` suggested some approximations which may be valid for concentrated distributions. They also introduced the concept of weighting results according to some reliability criteria. For the general case, however, it seems preferable to use the exact {cite:t}`kent1982` ellipses on uni-modal data sets. These could of course be weighted if such weighting is desired.
 
-:::{figure} ../figures/chapter12/love.png
+:::{figure} ../figures/chapter12/eigenvector.png
 :name: fig:love
-:alt: 3D plot with North-East-Down axes showing two antipodal clusters of red spheres representing a bi-Gaussian distribution of vectors.
+:alt: 3D plot with North-East-Down coordinate axes (dashed lines) showing two antipodal clusters of red spheres representing a bi-Gaussian distribution of vectors, with the principal eigenvector V1 (solid line) passing through both clusters.
 :width: 80%
 
-A bi-gaussian set of vectors suitable for treatment using the method of {cite:t}`love2003`.
+A bi-Gaussian set of vectors suitable for treatment using the method of {cite:t}`love2003`. The dashed lines show the North-East-Down coordinate axes. The solid line is the principal eigenvector $\V_1$ of the orientation matrix, which defines the axis about which the moment of inertia of the data is least.
 :::
 
 (sect:bigaussian)=
 ### The Bi-Gaussian Distribution
 
-Until now we have continued the Fisher assumption of unit vectors. As already mentioned, neglect of the vector strength can lead to bias. {cite:t}`love2003` began the hard work of incorporating intensity information into the parameter estimation problem. Their method can handle bi-modal spherical Gaussian data such as those shown in [](#fig:love). Estimation of the Love parameters are beyond the scope of this book. Moreover, many data sets are not spherically symmetric as already noted and the {cite:t}`love2003` approach must be generalized to elliptical, more "blade-like" data sets than the "cotton balls" currently treatable.
+Until now we have continued the Fisher assumption of unit vectors. As already mentioned, neglect of the vector strength can lead to bias. {cite:t}`love2003` began the hard work of incorporating intensity information into the parameter estimation problem. Their method can handle bi-modal spherical Gaussian data such as those shown in [](#fig:love). Estimation of the Love parameters is beyond the scope of this book. Moreover, many data sets are not spherically symmetric as already noted and the {cite:t}`love2003` approach must be generalized to elliptical, more "blade-like" data sets than the "cotton balls" currently treatable.
+
+(sect:naive-bootstrap)=
+## The Simple (Naïve) Bootstrap
+
+As we have mentioned, real data may be pathological in several respects including bi-modal and elliptically distributed data. None of the methods we have described so far have the test for common mean so critical to paleomagnetic studies nor can they provide confidence ellipses for an off-center mean direction as is likely to occur in records of the geomagnetic field (see [](#fig:vgp-di)b). Finally, data may be overprinted or contain the record of a paleomagnetic transition, resulting in "streaked" or non-antipodal distributions, conditions that make the conventional methods inappropriate. In this section, we will discuss alternative methods for estimating confidence bounds which are sufficiently flexible to accommodate all of these shortcomings, provided the data set is large enough.
+
+In [](#fig:hypeq)a we show a not unusual "not great" paleomagnetic data set. The data are elliptical, bi-modal and one has the suspicion that the normal and reverse modes may be neither antipodal nor share the same concentration or ovalness parameters. Clearly some non-parametric approach would be desirable. The approach for characterizing uncertainties for vectors we will take here is based on a technique known as the statistical *bootstrap*. As we shall see, the bootstrap has the flexibility to allow us to treat awkward data sets like that shown in [](#fig:hypeq)a.
 
 :::{figure} ../figures/chapter12/hypeq.png
 :name: fig:hypeq
@@ -91,13 +98,6 @@ Until now we have continued the Fisher assumption of unit vectors. As already me
 a) Hypothetical non-Fisherian data set. Normal and reversed polarity data that are not symmetrically distributed. Filled (open) circles plot on the lower (upper) hemisphere. b) Equal area projection of 500 bootstrapped means for pseudo-samples drawn from the data shown in a). c) Same as a) but with the bootstrapped confidence ellipses shown.
 :::
 
-(sect:naive-bootstrap)=
-## The Simple (Naïve) Bootstrap
-
-As we have mentioned, real data may be pathological in several respects including bi-modal and elliptically distributed data. None of the methods we have described so far have the test for common mean so critical to paleomagnetic studies nor can they provide confidence ellipses for an off-center mean direction as is likely to occur in records of the geomagnetic field (see [](#fig:vgp-di)b). Finally, data may be overprinted or contain the record of a paleomagnetic transition, resulting in "streaked" or non-antipodal distributions, conditions that make the conventional methods inappropriate. In this section we will discuss alternative methods for estimating confidence bounds which are sufficiently flexible to accommodate all of these short comings, provided the data set is large enough.
-
-In [](#fig:hypeq)a we show a not unusual "not great" paleomagnetic data set. The data are elliptical, bi-modal and one has the suspicion that the normal and reverse modes may be neither antipodal nor share the same concentration or ovalness parameters. Clearly some non-parametric approach would be desirable. The approach for characterizing uncertainties for vectors we will take here is based on a technique known as the statistical *bootstrap*. As we shall see, the bootstrap has the flexibility to allow us to treat awkward data sets like that shown in [](#fig:hypeq)a.
-
 The principles of the bootstrap are given in [](#app:bootstrap). In essence, the parameter of interest (say, the mean vector) is calculated from many resampled data sets, whose data points are selected at random from the original data. The bootstrapped estimates "map out" the likely distribution of the parameter, allowing estimation of confidence regions. Before we extend the bootstrap from the scalar treatment in [](#app:bootstrap) to vectors, it is important to point out that with the bootstrap, it is assumed that the underlying distribution is represented by the data, demanding that the data sets be rather large. Moreover, the bootstrap estimates are only asymptotically valid, meaning that a large number of bootstrap calculations are required for the confidence intervals to be valid. It's a good thing we have fast computers with huge hard-drives.
 
 There are a variety of ways we can use the bootstrap to estimate confidence regions for paleomagnetic data. We will start with the most "Fisher" like approach of taking unit vectors of a single polarity. Then we will accommodate dual polarity data sets and develop analogous tests to those so useful for Fisher distributions.