The Probability Density Function (PDF)

The information relevant to the kind of detrital thermochronology discussed in this paper is not so much the actual ages, but their probability distribution. Underlying any set of detrital ages is a Probability Density Function (PDF), describing the probability of occurence of any detrital age t:

$$Pr(a<t<b) = \int^b_a pdf(t) dt        \forall  a < b$$ (1)

In practice, the PDF can never be precisely known, because that would require an exhaustive sampling of the detrital population. Therefore, we must work with estimates of the PDF based on a finite sample of detrital ages (typically tens to hundreds of ages). Besides limited data, measurement uncertainty is a second factor reducing the precision of density estimates. The most popular estimators of probability density are histograms and ``kernel density plots'' [e.g., Silverman, 1986]. Both of these methods apply some degree of ``smoothing'' to the data, either by binning them into a histogram, or by assigning a Gaussian uncertainty distribution to each measurement:

$$\widehat{pdf}(t) = \frac{1}{n}\sum_{i=1}^{n}N (t\vert\mu = t_i, \sigma = \alpha \hat{\sigma}(t_i))$$ (2)

With N(t $\mid\mu,\sigma$) the normal distribution of t with mean $\mu$ and standard deviation $\sigma$, and t$_i$ and $\hat{\sigma}$(t$_i$) the measured ages and their respective 1-$\sigma$ uncertainties. To illustrate the different approaches to detrital thermochronological density estimation, consider the degenerate case of a ``diving board'' hypsometry: all detrital grains are derived from a single elevation, corresponding to a single ``true'' age t$_{true}$. The PDF of the true ages is a delta-function (spike at t$_{true}$, zero probability elsewhere). For further simplification, all grains have identical, Gaussian measurement uncertainties. In the following, t$_{true}$ = 10Ma so all grains are 10Ma old and have Gaussian measurement uncertainties of 1Ma.

Suppose we have access to an infinite number of measurement from this detrital population. The PDF of these age measurements can then be determined by a histogram with infinitessimal binwidth or a kernel density estimate with infinitessimal $\alpha$. Note that the PDF of the measurements is not the same as the PDF of the underlying ages (Figure 3.a). Unless we deconvolve the measurement uncertainties, the measurement distribution will always be a ``smoothed'' version of the ``true'' age distribution. In our toy example, the age distribution is a delta-function at 10 Ma, whereas the measurement distribution is a Gaussian distribution with mean 10 Ma and standard deviation 1 Ma (Figure 3.a).

Figure 3: (a) dashed line: expected (errorless) age distribution of a ``diving board'' hypsometry, black line: measurement distribution including 1Ma (normally distributed) measurement uncertainty; (b) circles are three measurements drawn from the measurement distribution. Using the nomenclature of Ruhl and Hodges [2005], SPDF$_z$ is the hypsometric curve and SPDF$_t$ is the Gaussian kernel density estimator of the measured age distribution; (c) as the number of measurements increases (50 in this figure), SPDF$_t$ converges to a Normal distribution with standard deviation $\sigma = \sqrt{2}$. SPDF$_{t*m}$, which has been used as the point of comparison between the hypsometric predictions and the measurements, is a Normal distribution with standard deviation $\sigma$ = 1.

Given a set of age data, the Gaussian kernel density estimator stacks a bell curve on top of each measurement (Equation 2 and Figure 3.b). Repeating this for a large number of measurements drawn from our ``diving board'' hypsometry yields a Gaussian distribution with mean 10 Ma and standard deviation $\sigma$ = $\sqrt{(1 Ma)^2 + (\alpha \times 1 Ma)^2}$. Thus, $\widehat{pdf}$ is ``double-smoothed'': once by the measurement uncertainties, and a second time by the construction of the kernel density estimator. The amount of additional smoothing depends on the parameter $\alpha$. Although it can be shown that $\alpha$ = 0.6 is an optimal value [ Silverman, 1986; Brandon, 1996], previous studies by Brewer et al. [2003], Ruhl and Hodges [2005] and Stock et al. [2006] have used $\alpha$ = 1, and so does Figure 3.b. Ruhl and Hodges [2005] gave this curve the name ``Synoptic Probability Density Function'' (SPDF). These authors distinguish between three kinds of SPDF. SPDF$_{z}$ is the true underlying age distribution, in the hypothetical case of errorless measurements (dashed lines in Figure 3). In our toy example, SPDF$_z$ is a delta function. SPDF$_{t}$ is the kernel density estimator generated by Equation 2 (gray lines in Figure 3). Finally, SPDF$_{t*m}$ effectively is the PDF of the measurements (black lines in Figure 3). Because SPDF$_{t*m}$ is only smoothed once, whereas SPDF$_{t}$ is smoothed twice, SPDF$_{t}$ is a biased estimator of SPDF$_{t*m}$ (Figure 3.c).

One of the requirements for the application of Gaussian kernel density estimation is that the measurement uncertainties are normally distributed. This may be a reasonable assumption for $^{40}$Ar/$^{39}$Ar thermochronology, but not necessarily for fission tracks, which are governed by a Poisson process. However, by using the logistic transform, a set of fission track data can be recast in terms of a new parameter z [Brandon, 1996], which is estimated by

$$\hat{z} = ln(\lambda \zeta g \rho_D) + ln(\frac{N_s + 0.5}{N_i + 0.5})$$ (3)

where $\lambda$ is the decay constant of $^{238}$U (=1.55125$\times$10$^{-10}$a$^{-1}$), $\zeta$ a (zeta) calibration factor measured on an AFT age standard, g a geometric factor (=0.5), N$_s$ the number of spontaneous fission tracks, N$_i$ the number of induced tracks in a mica detector and $\rho_D$ the induced track density of a glass standard that was irradiated along with the sample. N$_s$ and N$_i$ are Poisson variables, but $\hat{z}$ is normally distributed with standard error

$$\hat{\sigma}(z) = \sqrt{\frac{1}{N_s+0.5} + \frac{1}{N_i+0.5}}$$ (4)