NOAA / NESDIS / National Geophysical Data Center / WDC for SEG, Boulder, Colorado USA

Example 4.1

1) Use the program fisher to generate a set of 10 points drawn from a Fisher distribution with a kappa of 15. 2) Make an equal area projection with eqarea.

Solution

To generate a Fisher distributed data set with N=10 and kappa=15, save it to a file named ex4.1 and see what is in it, type the following:

% fisher -kns 15 10 66 > ex4.1

Check out what's in ex4.1 using cat.

The -s switch is an integer used as a seed for a random number generator. Different distributions can be made with different values of -s. To make a postscript file of an equal area projection of these data and view it on the screen, type:

% fisher -kns 15 10 66 | eqarea | plotxy

This causes plotxy to make a postscript file named mypost.

Example 4.2

Calculate a mean direction, k, a95 from the distribution generated in Example 4.1 using the program gofish. Repeat for the principal direction using goprinc.

Solution

Type the following:

% gofish < ex4.1

and the computer responds:

318.1 88.9 10 9.1220 10.3 15.8

To find out what these numbers are, check the documentation for gofish.

As a short cut, you could take advantage of UNIX's pipe facility by:

% fisher -kns 15 10 66 | gofish

which does the same thing, without the intermediate file ex4.1.

Now type:

% goprinc < ex4.1

to which the computer reponds:

33.9 89.4 10 0.837

which are the principal directions D, I, N and tau₁ respectively.

Example 4.3

1) Use the program fishdmag to calculate a mean direction from the last seven data points of sample tst1a in the file ex4.3. 2) Use the program pca to calculate the principal component direction from the last 13 data points of sample tst1b. 3) Use program gtcirc to calculate the best-fitting great circle from the last 22 data points of sample tst1c. 4) Finally, use the program lnp to calculate the mean direction and a95 using the two directed lines and the great-circle data for site tst1.

Solution

Type:

% grep tst1a ex4.3 | fishdmag -f 7 13 > ex4.3a % grep tst1a ex4.3 | pca -p 6 17 >> ex4.3a % grep tst1a ex4.3 | gtcirc -g 2 23 >> ex4.3a

The columns are: sample name, f/p/g, first step, last step, a95, MAD_plane, declination, inclination, where the ``f'' stands for Fisher mean, ``p'' stands for principle component, and ``g'' stands for the pole to best-fit plane (great circle). Check what is in the output file: ex4.3a. To calculate an estimate of the Fisher mean of combined lines and planes, type:

% lnp < ex4.3a

and the computer responds:

tst1 2 1 16.7 166.7 -60.1

To find out what these numbers are, check the documentation for lnp.

Example 4.4

Use program fishrot to draw a set of 50 data points from a Fisher distribution with a mean declination of 0, inclination of 20, and a kappa of 30. Calculate the Fisher statistics of the data set with the program gofish. Now awk out the inclinations and estimate the mean and 95% confidence bounds using the program incfish. Repeat this for populations having a mean inclination of 40, 60 and 80. How well does the inclination-only method work for high inclinations?

Solution

The following command will generate the desired distribution:

% fishrot -kndi 30 50 0 20 > ex4.4

Calculate the Fisher statistics using:

% gofish < ex4.4

and the computer responds:

2.3 21.3 50 48.3152 29.1 3.8

Now type the following to select the inclination data (second column) and pipe them to incfish.

% awk '{print $2}' ex4.4 | incfish

and the computer responds:

21.0 24.6 17.3 50 27. 3.7

Check the documentation for incfish (or type incfish -h) to find out what these numbers mean.

Note that the estimated mean inclination using incfish is 21.4 and that calculated using gofish is 21.3. To see what happens at higher inclinations, we can streamline things a bit using the tee utility. The tee utility transcribes the standard input to the standard output and copies it to a specified file name, as follows procedure.

% fishrot -kndi 30 50 0 40 | tee ex4.4a | gofish

to which the computer responds:

358.1 40.1 50 48.4065 30.7 3.7

and silently makes ex4.4a containing the fishrot output. Now type:

% awk '{print $2}' ex4.4a | incfish

and get:

39.3 42.7 35.8 50 32. 3.4

also in good agreement with the gofish calculation.

Repeat for an inclination of 60:

% fishrot -kndi 30 50 0 60 | tee ex4.4b | gofish

to which the computer responds:

357.0 60.1 50 48.4056 30.7 3.7

and

% awk '{print $2}' ex4.4b | incfish

to get:

58.3 61.9 54.8 50 34. 3.6

Finally, try an inclination of 80:

% fishrot -kndi 30 50 0 80 | tee ex4.4c | gofish

to which the computer responds:

351.4 80.0 50 48.4076 30.8 3.7

and

% awk '{print $2}' ex4.4c | incfish

to get the message:

************Result no good***************

This illustrates the point that the method breaks down at high inclinations.

Example 4.5

Use program fishqq to check if a data set produced by fisher is likely to be Fisherian.

Solution

Type the following to draw 25 directions from a Fisher distribution having a kappa of 25:

% fisher -kns 25 25 44 > ex4.5

(To draw a different distribution, you can change the random seed from 44 to any non-zero integer.)

Try plotting the data in equal area projection with the command

% eqarea < ex4.5 | plotxy

(see mypost file).

To make a Q-Q plot against a Fisher distribution, type

% fishqq < ex4.5 | plotxy

This produces a Q-Q plot (mypost).

Example 4.6

Use the program gauss to generate a data set by drawing 200 data points from a normal distribution with a mean of 22, a standard deviation of 5. Plot these data as a histogram with histplot and calculate a mean and standard deviation using stats.

Solution

Type the following:

% gauss -msni 22 5 200 22 > ex4.6

to generate a file with the normal data in it. Now type

% histplot < ex4.6 | plotxy

to cause plotxy to generate mypost.

Finally, type

% stats < ex4.6

to get:

200 22.0720 4414.41 4.89825 22.1921 0.346358 0.679963

To know what these numbers are, type stats -h or check the on-line documentation: stats.

Example 4.7

Use the program qqplot to plot a Q-Q plot of the data in ex4.6, calculate the Kolomogorov-Smirnov statistic D, and the 95% asymptotic significance level for N data points.

Solution

Type:

% qqplot < ex4.6 | plotxy

and view the postscript file mypost.

Example 4.8

Use the program bootstrap to calculate a bootstrap confidence interval for the data generated in Example 4.6.

Solution

Type the following:

% bootstrap -p < ex4.6 | plotxy

To generate a distribution of 1000 bootstrapped means. (This is the default number; a different number of boostraps can be selected using the -b switch). The program spits out plotxy commands for postscript picture mypost.

Example 4.9

Use program bootdi to see if the data set in ex4.9 are 1) likely to be Fisherian, 2) what the approximate 95% confidence ellipses are (try both simple and parametric bootstraps), 3) what the approximate 95% confidence ellipses are using the principal eigenvectors instead of Fisher means.

Solution

Type:

% bootdi < ex4.9

to which the computer responds:

Total N = 49
Mode:	Dec,	Inc,	a95,	N,	kappa,	Fisherian ?
1st:	26.5	39.1	10.7	25	8.	no
2nd:	189.2	-44.3	8.3	24	14.	no
Mode	eta,	dec,	inc,	zeta,	dec,	inc
1st:	6.61	239.44	45.93	11.79	130.96	17.06
2nd:	4.59	152.36	39.38	10.64	259.14	19.38

For a parametric bootstrap, type:

% bootdi -p < ex4.9

to which the computer responds:

Total N = 49
Mode:	Dec,	Inc,	a95,	N,	kappa,	Fisherian ?
1st:	26.5	39.1	10.7	25	8.	no
2nd:	189.2	-44.3	8.3	24	14.	no
Mode	eta,	dec,	inc,	zeta,	dec,	inc
1st:	6.71	240.24	45.67	11.54	131.39	17.52
2nd:	5.30	152.44	39.41	11.25	259.20	19.34

To work on the principal eigenvectors, type:

% bootdi -P < ex4.9

to which the computer responds:

Total N = 49
Mode:	Dec,	Inc,	a95,	N,	kappa,	Fisherian ?
1st:	30.8	40.9	10.7	25	8.	no
2nd:	9.8	44.1	8.3	24	14.	no
Mode	eta,	dec,	inc,	zeta,	dec,	inc
1st:	6.26	222.02	48.52	9.77	125.69	5.57
2nd:	4.70	150.84	38.69	10.99	258.38	20.62

Example 4.10

Use plotdi to make an equal area projection of the data and bootstrap confidence ellipses of the data in ex4.9. Make a plot the bootstrapped eigenvectors.

To plot the data and the bootstrap ellipses, type:

% plotdi < ex4.9 | plotxy

See mypost.

To plot the data and the bootstrap eigenvectors, type:

% plotdi -v < ex4.9 | plotxy

See mypost.

The command bootdi -v < ex4.9 will generate the list of bootstrapped eigenvectors used in plotdi -v.

Example 4.11

Assume that the sampling site for the data in ex4.9 was at latitude 35 N and longitude 120 W. 1) Calculate the direction expected from a dipole field at that location using the dipole formula (tan (I) = 2 tan (lat.)) (see Chapter 1 of Paleomagnetic Principles and Practice). 2) Use cart_hist to compare the data set in ex4.9 with the expected dipole direction. 3) The data were taken in 1997 near sealevel. Use igrf to calculate the expected geomagnetic field direction. Use cart_hist to compare the data with this direction. 4) Use cart_hist to perform a parametric bootstrap reversals test.

Solution

1) The expected dipole direction at the site is D=0,I=54.5.

2) Type the following to find out about cart_hist:

% cart_hist -h

or check the on-line documentation cart_hist.

So to compare a set of directions with a known direction, type:

% cart_hist -d 0 54.5 < ex4.9 | plotxy

and see the file mypost.

3) To determine the IGRF direction at the site (see also Example 1.3), type the following:

% igrf

1997 0 35 -120

(and control-D)

and the computer responds:

14.4 59.5 48925

4) The geomagnetic reference field at the sampling site was therefore D=14.3, I= 59.5, so, to compare with the data, type:

% cart_hist -d 14.4 59.5 <ex4.8 | plotxy

and see mypost.

Example 4.12

Pretend that a reference pole was determined for the same age as the data in ex4.9. with a position of longitude = -140 latitude = 75. The VGPs that went into the calculated had a reported kappa of 30. 1) Assume that the VGPs were Fisher distributed and generate a synthetic data set with the same mean and kappa and N = 100 using fishrot. 2) Convert these data to expected directions at the sampling site for last example. 3) Compare the data in ex4.9 to those just generated. Are they significantly different?

Solution

Since you don't have access to the original data, you can assume that the VGPs were Fisher distributed and create a synthetic distribution with the same mean direction and kappa.

Since the default kappa, N, and seed of fishrot are acceptable, we need only specify the VGP position by:

% fishrot -di -140 75 > ex4.12a

Now we have to put in the site location, rearrange the pole latitude and longitude and pipe it to vgp_di (see Example 1.6).

% awk '{print $2,$1,35,-120}' ex4.12a | vgp_di > ex4.12b

To complete the problem, we must compare the data in ex4.9 with those in ex4.12.b using cart_hist. We also need to flip the reversed data and plot the 95% confidence bounds.

Type:

% cart_hist -crb ex4.12b ex4.9 | plotxy

and look at mypost.

It appears that the data in ex4.9 are significantly different from the dipole field, the present geomagnetic field and the reference directions.

Example 4.13

Use program foldtest to do a parametric foldtest on the data in ex4.13.

Solution

Type the following:

% foldtest -p < ex4.13 | plotxy

which produces the file mypost.