C Paleomagnetic statistics and parameter estimation

Appendix C
Paleomagnetic statistics and parameter estimation

Chapters 5 and 7 discussed various hysteresis parameters and Chapters 11 and 12 developed the major features of paleomagnetic directional statistics. Here we go over some aspects in greater detail.

C.1 Hysteresis Parameters

A typical hysteresis experiment involves determination of a hysteresis loop and frequently also a back-field curve (see Figure C.1). Processing of the data in the PmagPy software package (see, e.g., Example for hysteresis_magic.py) proceeds as follows:

  1. Sometimes the descending and ascending hysteresis loops do not close because of instrument drift (see Figure C.1a). Ordinarily, the experiment should be re-done, but for small differences, we force the loops to close by subtracting the difference, interpolated from the maximum difference at the maximum field (Bmax) to a zero difference at the minimum field (Bmin).
  2. After closing the loops, we calculate the best-fit line to the M,B data for the portion within 70% of Bmax, averaging data from both the ascending and descending loops. A difference in the absolute value of the y-intercepts for the ascending and descending loops indicates a vertical offset of the data, which is adjusted such that the two intercepts are equal. The average slope is the high-field susceptibility (χhf), which is subtracted off. The data after these steps are shown as the dashed line in Figure C.1. The maximum magnetization after adjusting for the χhf is the saturation magnetization Ms.

    Table C.1: Summary of hysteresis parameters.

    Symbol Method Section Figure

    χhf high field susceptibility  5.6 &  C.1
    Ms saturation magnetization  3.2.2  5.6 &  C.1
    Mr saturation remanence  5.6 &  C.1
    Hc or μoHcCoercivity  4.1.3 & &  C.1
    Coercivity of remanence:
    Hcr ΔM method  5.2.1  C.1
    Hcrascending loop intercept method 5.2.1  5.6
    Hcr′′ Back-field method  7.7,  7.20C.1
    Hcr′′′ H12 method  7.7 & 8.4.2  7.20 &  8.9

    Figure C.1: Typical hysteresis experiment. a) Raw data are solid red line. Data are processed (see text) by closing the ascending and descending loops, subtracting the high field slope (χhf) and adjusting such that the y-intercepts are equal (that for the descending loop is labeled Mr). Processed data are dotted blue line. Coercivity (μoHc) is the applied field for which a saturation magnetization (Ms) is reduced to zero. b) Difference between processed ascending and descending loops is the ΔM curve (solid blue line). Back-field IRM data shown normalized by saturation remanence (Mr) – dashed green line. Two methods of estimating coercivity of remanence shown (see text).

  3. Coercivity (μoHc) is the field at which M = 0. We estimate this by finding the values of B between which M switches sign for both the ascending and descending loops (after adjustment), calculate a line and evaluate the B for which M = 0. The coercivity is the average of the two estimates.
  4. We fit a spline to the adjusted ascending and descending loops and resample the loops at even intervals of B (usually 10 mT intervals). The ΔM curve shown in Figure C.1b is the difference between these two interpolated curves, averaging the data for negative and positive B. The saturation remanence Mr is the value of the ΔM curve at B = 0. The coercivity of remanence (μoHcr in Table C.1) is the field for which ΔM is half the value of Mr. This is the “ΔM” method of coercivity of remanence calculation (see Chapter 5).
  5. If there are “back-field” IRM data as in Figure C.1b, the coercivity of remanence can be estimated by finding (through interpolation) the applied field which reduces the saturation remanence (Mr) to zero. This is the “back field” method.

C.2 Directional statistics

Table C.2: Critical values of Ro for a random distribution [Watson, 1956.]

N95%99% N95%99%


C.2.1 Calculation of Watson’s V w

  1. Calculate Ri, and ki where i = 1,2 for the two data sets with N1,N2 samples using Equations 11.6 and  11.8.
  2. Calculate xij (where j = 1,3 for the three axes) using Equation 11.7.
  3. Calculate Xij = Rixij.
  4. Find the weighted means for the two data sets:

Xˆj =     ki �Xij.
  5. Calculate the weighted overall resultant vector Rw by
           ˆ2   ˆ 2   ˆ2 12
Rw =  (X1 + X 2 + X 3) ,
    and the weighted sum Sw by,
Sw =    kiRi.
  6. Finally, Watson’s V w is defined as
    Vw = 2(Sw - Rw ).

C.2.2 Combining lines and planes

  1. Calculate M directed lines (two in our case) and N great circles (one in our case) using principal component analysis (see Chapter 9) or Fisher statistics.
  2. Assume that the primary direction of magnetization for the samples with great circles lies somewhere along the great circle path (i.e., within the plane).
  3. Assume that the set of M directed lines and N unknown directions are drawn from a Fisher distribution.
  4. Iteratively search along the great circle paths for directions that maximize the resultant vector R for the M + N directions.
  5. Having found the set of N directions that lie along their respective great circles, estimate the mean direction using Equation 11.7 and κ as:

    k = -2M--+-N----2-,
    2(M  + N - R )
    The cone of 95% confidence about the mean is given by:

cosα95 = 1-  N----1,[(1-)1∕(N′-1) - 1],
              kR     p
    where N= M + N∕2 and p = .02

Table C.3: Maximum likelihood estimators of k1,k2 in the Bingham distribution for given eigenvalues ω12. Data from Mardia and Zemroch (1977). Upper (lower) number is k1(k2)

ω1 0.02 0.04 0.06 0.08 0.10 0.12 0.14 0.16 0.18 0.20 0.22 0.24 0.26 0.28 0.30 0.32

0.08 -25.6-13.16-9.065-7.035
0.20-25.64-13.18 -9.05-6.974-5.694-4.799-4.118-3.570-3.109-2.709
-1.312-1.267-1.216-1.159-1.096-1.028-0.955-0.876-0.791-0.701-0.604-0.500-0.389-0.269-0.140 0.000
-1.123-1.073-1.017-9.555-0.887-0.814-0.736-0.651-0.561-0.464-0.360-0.249-0.129 0.000
-0.940-0.885-0.824-0.757-0.684-0.606-0.522-0.432-0.335-0.231-0.120 0.000
0.42 -25.5-12.73-8.532-6.388-5.045-4.089-3.349-2.741
-0.589-0.523-0.452-0.374-0.290-0.200-0.104 0.000
-0.418-0.347-0.270-0.186-0.097 0.000
-0.250-0.173-0.090 0.000

C.2.3 Inclination only calculation

We wish to estimate the co-inclination (α = 90 -I) of N Fisher distributed data (αi), the declinations of which are unknown. We define the estimated value of α to be ˆα. McFadden and Reid showed that αˆ is the solution of:

                         ∑                     ∑
N cos ˆα + (sin 2ˆα - cos2ˆα )  cos αi - 2sinαˆcos ˆα  αi = 0,

which can be solved numerically.

They further define two parameters S and C as:

S =     sin(ˆα - αi),
C  =     cos (ˆα- αi).

An unbiassed approximation for the Fisher parameter κ, k is given by:

     N - 1
k =---------.
   2(N - C )

The unbiased estimate of the true inclination is:

ˆ           S-
I = 90 - ˆα + C .

Finally, the α95 is estimated by:

             1 S  2    f
cosα95 = 1 - 2(C-) -  2Ck-,
where f is the critical value taken from the F distribution (see F-distribution tables in statistics textbooks or online) with 1 and (N-1) degrees of freedom.

C.2.4 Kent 95% confidence ellipse

Kent parameters are calculated by rotating unimodal directions x into the data coordinates xby the transformation:
 ′   T
x = Γ  x,

where Γ = (γ123), and the columns of Γ are called the constrained eigenvectors of orientation matrix, T (see Appendix A.3.5). The vector γ1 is parallel to the Fisher mean of the data, whereas γ2 and γ3 (the major and minor axes) diagonalize T as much as possible subject to being constrained by γ1 (see Kent, 1982, but note that his x1 corresponds to x3 in conventional paleomagnetic notation). The following parameters may then be computed:
  ˆμ = N -1  kxk1′
ˆσ2 = N -1∑  (xk2′)2
ˆσ22=  N -1∑ k(x  ′)2.
 3         k  k3

As defined here, ˆμ = R∕N (R is closely approximated by the equation for R in Chapter 11. Also to good approximation, σˆ22 = τ2, and ˆσ32 = τ3, where τi are the eigenvalues of the orientation matrix. The semi-angles ζ95 and η95 subtended by the major and minor axes of the 95% confidence ellipse are given by:
             √ --               √ --
ζ95 = sin-1(σ2 g),  η95 = sin-1(σ3 g),

where g = -2 ln(0.05)(Nˆμ2).

The tensor Γ is, to a good approximation, equivalent to V, the eigenvectors of the orientation matrix. Therefore, the eigenvectors of the orientation matrix V give a good estimate for the directions of the semi-angles by:
D ζ = tan- 1(v22∕v12),Iζ = sin- 1v32,
D η = tan- 1(v23∕v13),Iη = sin -1v33,

where for example the x2 component of the smallest eigenvector V3 is denoted v23.

C.2.5 Bingham 95% confidence parameters

The Bingham distribution is given by:

         1               2         2     2
F =  4πd(k-,k-) exp (k1cos ϕ + k2sin ϕ)sin α,
          1  2
where α and ϕ are as in the Kent distribution, k1,k2 are concentration parameters (k1 < k2 < 0) and d(k1,k2) is a constant of normalization given by:

              ∫   ∫
           -1-  2π  π            2         2     2
d(k1,k2) = 4π  0   0  exp ((k1 cos  ϕ+ k2 sin  ϕ)sin  θ)sin θdθϕ.

To estimate the axes of the Bingham confidence ellipse, we first calculate the eigenparameters of the orientation matrix as for Kent parameters and described in Appendices A.3.5 and  C.2.4. The principle eigenvector V1 of the orientation matrix is associated with the largest eigenvalue τ1. In Bingham (1974), ω1 is the τ3 and ω3 is τ1. In Bingham statistics, the V1 direction is taken as the mean. Beware – it is not always parallel to the Fisher mean of a unimodal set of directions.

The maximum likelihood estimates of k1,k2, the concentration parameters are gotten by first maximizing the log likelihood function:

F = - N log(4π)-  N log d(k1,k2)+ k1ω1 + k2ω2.
These are listed for convenience in Table C.3 as calculated by Mardia and Zemroch (1977). Once these are estimated, the semi-axes of the 95% confidence ellipse around the mean direction V1 are given by:

2     2    2
ϵij = χp(ν)σij,
where χp2(ν) = 5.99 is the χ2 value for significance (p = .05 for 95% confidence) with ν = 2 degrees of freedom and

σ2ij = -------------------.
      2N (ωi - ωj)(ki - kj)

Bingham (1974) set k3 = 0, so the semi axes of the confidence ellipse about the principle direction V1, associated with ω3, are therefore:

ϵ32 = --k-N-(ω---ω-),
        2    3    2
ϵ31 =---------------.
     - k1N (ω3 - ω1)

Because k1 < k2 < 0, the semi-axes are positive numbers. Please note that here we use the corrected version of Tanaka (1999) as opposed to the more oft-quoted but erroneous treatment of Onstott (1980). Note also that the N is required for σ because we have normalized the ωs to sum to unity for consistency with other eigenvalue problems in this book. The N is missing in the treatment of Tanaka (1999) presumably because the eigenvalues sum to N. Finally, note that these values of ϵ are in radians and must be converted to degrees for most applications.

C.3 Paleointensity statistics

Paleointensity statistics have gotten somewhat out of hand of late. There are a bewildering variety of statistics that are used in the literature, with no consensus as to which ones are essential, which ones are helpful and which ones are irrelevant. This appendix will not help the reader in this regard, but merely attempts to assemble the ones we feel are the most useful.

Figure C.2: Illustration of paleointensity parameters. Arai plots: The magnitude of the NRM remaining after each step is plotted versus the pTRM gained at each temperature step. Closed symbols are zero-field first followed by in-field steps (ZI) while open symbols are in-field first followed by zero field (IZ). Triangles are pTRM checks and squares are pTRM tail checks. Horizontal dashed lines are the vector difference sum (VDS) of the NRM steps. Vector endpoint plots: Insets are the x,y (solid symbols) and x,z (open symbols) projections of the (unoriented) natural remanence (zero field steps) as it evolves from the initial state (plus signs) to the demagnetized state. The laboratory field was applied along -Z. Diamonds indicate bounding steps for calculations. a) The fvds is the fraction of the component used of the total VDS. The difference between the pTRM check and the original measurement at each step is δTi. The inset shows the deviation angle (DANG) that a component of NRM makes with the origin. The maximum angle of deviation MAD is calculated from the scatter of the points about the best-fit line (solid green line). b) Data exhibit zig-zag behavior diagnostic for significant difference between blocking and unblocking temperatures. The Zig-zag for slopes compares slopes calculated between ZI and IZ steps (bzi) with those connecting IZ and ZI steps biz). The difference between the pTRM tail check and the original measurement at each step is ΔTi. c) β reflects the scatter (δxy) about the best-fit slope (solid green line). The Zig-zag for directions compares those calculated between ZI and IZ steps (Dzi) with those connecting IZ and ZI steps Diz). [Figures from Ben-Yosef et al., 2008.]

  1. The Deviation of the ANGle (DANG; Tauxe and Staudigel, 2004; see Chapter 9): The angle that the direction of the NRM component used in the slope calculations calculated as a best-fit line (see Appendix A.3.5) makes with the angle of the line anchoring the center of mass (see Appendix A.3.5) to the origin (see insert to Fig. C.2a).
  2. The Maximum Angle of Deviation (MAD; Kirschvink, 1980; see Chapter 9): The scatter about the best-fit line through the NRM steps.
  3. We can calculate the best-fit slope (b) for the data on the NRM-pTRM plot and its standard error σ (York, 1966; Coe et al. 1978). The procedure for calculating the best-fit slope, which is the best estimate for the paleofield, is given as follows:

    a)Take the N data points that span two temperature steps T1 and T2, the best-fit slope b relating the NRM (yi) and the pTRM (xi) data in a least squares sense (taking into account variations in both x and y is given by:
         ∘ -∑----------
b = -   ∑-i(yi---�y)2,
          i(xi - �x)2

    where y is the average of all y values and x is the average of all x values.

    b) The y-intercept (yo) is given by y - bx.

    c) The standard error of the slope σ is:
         ∘ -∑---------------∑-----------------
       2---i(yi --y�)2---2b-i(xi --�x)(yi --�y)
σb =           (N - 2)∑  (xi - �x )2       .

  4. The “scatter” parameter β: the standard error of the slope σ (assuming uncertainty in both the pTRM and NRM data) over the absolute value of the best-fit slope |b| (Coe et al. 1978).
  5. The remanence fraction, f, was defined by Coe et al. (1978) as:

    f = ΔyT ∕yo,
    where ΔyT is the length of the NRM segment used in the slope calculation (see Figure C.2).
  6. The fraction of the total remanence (by vector difference sum), fvds (Tauxe and Staudigel, 2004): While f works well with single component magnetizations as in Fig. C.2d where it reflects the fraction of the total NRM used in the slope calculation, it can be misleading when there are multiple components of remanence as in Fig. C.2a. The values of f for such specimens can be quite high, whereas the fraction of the total NRM is much less. We prefer to use a parameter fvds which is the fraction of the total NRM, estimated by the vector difference sum (VDS; Chapter 9) of the entire zero field demagnetization data. The VDS (see Fig. C.2a) “straightens out” the various components of the NRM by summing up the vector differences at each demagnetization step. fvds is calculated as:

    fvds = ΔyT ∕yvds,

    where yvds is the vector difference sum of the entire NRM (see Figure C.2a and Chapter 9). This parameter becomes small, if the remanence is multi-component, whereas the original f can be blind to multi-component remanences.

  7. The Difference RATio Sum, DRATS: The difference between the original pTRM at a given temperature step (horizontal component of the circles in Fig. C.2) and the pTRM check (horizontal component of the triangles in see Fig. C.2), δi (see Fig. C.2a), can result from experimental noise or from alteration during the experiment. Selkin and Tauxe (2000) normalized the maximum δi value within the region of interest by the length of the hypotenuse of the NRM/pTRM data used in the slope calculation. DRAT is therefore the maximum difference ratio expressed as a percentage. In many cases, it is useful to consider the trend of the pTRM checks as well as their maximum deviations. We follow Tauxe and Staudigel (2004) who used the sum of these differences. We normalize this difference sum by the pTRM acquired by cooling from the maximum temperature step used in the slope calculation to room temperature. This parameter is called the Difference RATio Sum or DRATS. Only pTRM checks at temperatures below the maximum bound are included in the DRATS calculation.
  8. Maximum Difference % MD%: The absolute value of the difference between the original NRM measured at a given temperature step (vertical component of the circles in Fig. C.2) and the second zero field step (known as the pTRM tail check) results from some of the pTRM imparted in the laboratory at Ti having unblocking temperatures that are greater than Ti. These differences (Δi; see Fig. C.2b) are plotted as squares. The Maximum Difference, normalized by the VDS of the NRM and expressed as a percentage is the parameter MD%.
  9. Zig-Zag Z: In certain specimens, the IZZI protocol leads to rather interesting behavior, described in detail by Yu et al. (2004). The solid symbols in Fig. C.2 are the zerofield-infield (ZI) steps and the intervening steps are the infield-zerofield (IZ) steps (open circles). Alternating the two results in a “zigzag” in some specimens. The zigzag can be in either the Arai diagrams (compare slope of solid versus dashed line segments in Fig. C.2b) or in the orthogonal projections or the zero field vectors (compare directions of solid and dashed line segments in Fig. C.2 c). We therefore can define a parameter Z by testing the difference in either the two sets of slopes or the two sets of directions between the IZ steps and the ZI steps.

    To test the significance of the difference between the zero-field IZ directions and those from the ZI zero-field steps, we calculate F-test (Fw) for Watson’s test for common mean (Watson, 1983). The zigzag for directions Zdir is ratio Fw∕F(ν,) where F(ν) is the critical value for F at ν = 2N - 2 degrees of freedom (at the 95% level of confidence).

    For the slopes, we calculate the mean and variance of the slopes for the IZ segments (biziz2) and the ZI segments (bzizi2). The parameter tb is the t test for the two means. The zigzag for the slopes Zslope is the ratio tb∕t(ν) where t(ν) is critical value for t with ν = Niz + Nzi - 2 degrees of freedom (from a statistics table).

    If the difference between the sets of directions and slopes is less than 2 or both Zslope and Zdir are less than unity, then Z = 0. Otherwise Z is the larger of Zdir and Zslope.

  10. The “gap factor” g (Coe et al. 1978) penalizes uneven distribution of data points and is:
    g = 1 - �Δ �y∕Δy  ,

    where Δy is given by :

              i=∑N -1
Δ��y = -1---     Δy2i,
      ΔyT   i=1
    and is the weighted mean of he gaps Δyi between the N data points along the selected segment.
  11. The Coe quality index q combines the standard error of the slope, the NRM fraction and the gap factors by:

    q = βf g.

    As data spacing becomes less uniform, g decreases.

  12. A quick and dirty test for the possibility of anisotropy of TRM is to compare the direction of the pTRM acquired in the laboratory field with the direction of the applied field. The angle between these two at the maximum pTRM temperature used in the slope calculation is defined here as γ. If this exceeds more than a few degrees, it is advisable to perform some sort of test for TRM (or ARM) anisotropy.