Chapter 5.15. Bayesian techniques: an overview

The analysis of X-ray absorption fine-structure (XAFS) data is traditionally based on least-squares fitting methods. However, if there are more parameters than the data alone allow to be determined, this leads to an ill-conditioned system of equations (Stern, 1988; Rehr & Albers, 2000). We instead propose a Bayesian approach (Krappe & Rossner, 1985a,b, 2004). This allows the subspace of the total model-parameter space in which the data decide predominantly upon the outcome of the fit to be determined, whereas in complementary space a priori assumptions essentially fix the result.

We start the discussion with the direct problem, in which the cross section for photoelectron absorption μ_exp(k) on the K edge or an isolated L edge is calculated as a function of the wavenumber k of the absorbed X-ray. The calculation is based on the single-scattering formula for monoatomic, unoriented samples or samples with cubic symmetry, valid for k > k_cut, where the energy k²ℏ²/(2m) is sufficiently larger than the threshold energy for the photoeffect (Stern, 1988; Zabinsky et al., 1995; Bunker, 1983; Tröger et al., 1994), with the length correction δR_j = δR_j∥ + δR_j⊥ (Fornasini et al., 2004; Rossner et al., 2006), where and truncating the sum beyond the Jth term, related to the parameter k_cut. The model parameters in this approach are the following: the correction factor S₀ for many-electron effects, the half-path distances of the J scattering paths, R_j, j = 1, …, J, the projected Debye–Waller (DW) parameters and the anharmonicity parameters, and C_3;j, respectively. The correction factor S₀ does in fact depend slightly on the path j and on k. We define as the average of the actual averaged over j and over k in the relevant k-range, with the remaining k and j dependence being absorbed in the scattering amplitudes f_j(k, R_j). The absorption coefficient for the embedded absorbing atom μ₀(k), the amplitudes f_j(k, R_j), the scattering phases φ_j(k) and the damping parameter λ(k) follow from an (approximate) solution of the n-electron scattering problem and are, for example, calculated by the FEFF code (Zabinsky et al., 1995; Ankudinov, 1996; Rehr & Albers, 2000; Kas et al., 2024).

To solve the inverse problem, that is to determine the model parameters from a given set of measured values μexp(E_l) at energies E_l, we first obtain wavenumbers k_l = [2sm(E_l − E₀)]^1/2/ℏ, where E₀ is the effective threshold energy. Since the latter is known only approximately, it is often treated as another of the model parameters to be determined by the fit. As usual, a smooth background contribution μ_back(k) is subtracted from μ_exp(k), which is obtained from a polynomial extrapolation of the pre-edge μ_exp to the post-edge region as described by Victoreen (1948) and Milledge (1962). Unfortunately, the precise extrapolation recipe influences the final fit parameters somewhat. The μ_exp(k) are measured with an energy-dependent efficiency A(k). As in Krappe & Rossner (2004), we obtain from μ_exp(k) by a polynomial smoothing procedure, described in detail in the appendix to Krappe & Rossner (2004). Similarly, is obtained from μ₀(k). The ratio A(k) = for k > k_edge is then interpreted as the overall efficiency of the experimental setup in the EXAFS energy range k > k_cut.

The FEFF result for μ₀ needs corrections (Krappe & Rossner, 2004). We therefore write μ₀(k) = + , where δμ₀(k) is represented by a cubic spline on an equally spaced grid of support points, the number of which T is to be chosen to make the spline just sufficiently flexible for the purpose for which it is introduced (Krappe & Rossner, 2004). The ordinates δμ_t, t = 1, …, T are also treated as model parameters to be determined in the fit together with all other model parameters.

We have therefore to fit the function for l = 1, …, L, where χ is given for k > k_cut by equation (1). The set of model parameters is

We give the experimental data a Gaussian probability distribution , characterized by the quadratic form It is usually assumed that the matrix F, which is the inverse of the variance matrix, is diagonal: . One also has to associate errors with the FEFF code because the electron multiple-scattering problem can only be treated approximately, for instance by including in the sum in equation (1) only terms which contribute more than 4% of the total, and integrals have to be approximated by finite sums. We again associate with these errors a Gaussian probability that z′ is true for a given x, i.e. , with where the matrix B is the inverse of the variance matrix.

The conditional probability that the outcome of the observation is , once x is given, may be expressed in terms of and by The integral can be evaluated analytically and yields a Gaussian in g(x), , with in terms of the L × L matrix

We expand g(x) around a first-guess value x⁽⁰⁾ for the solution of the inverse problem where the L × N matrix G is defined as Inserting into equation (6) and calling x − x⁽⁰⁾ in the following x to simplify the notation, one obtains a second-order polynomial in x, in terms of the N × N matrix Q = G^TCG and the vector . The matrix Q is the inverse of the variance matrix of the distribution P_cond in terms of the variable x.

In order to find the probability distribution for the parameter values x, once the are given, , we use Bayes' theorem Bayes' theorem therefore solves the inversion problem in probability theory, but at the price of introducing the prior probability P_prior(x), which expresses the knowledge that we have about the model parameters before the experiment is made. Let us assume for the moment that we have an average value x^(prior) and a variance matrix A⁻¹ so that with and let us further restrict the matrix A to be diagonal, and choose x^(prior) = x⁽⁰⁾. Maximizing P_post yields the normal equations

An optimal choice of the diagonal matrix A must obviously take the quality of the data into account. Turchin and Nozik (Turchin & Nozik, 1969; Turchin et al., 1970; Turchin, 1985) assume that there is a probability distribution of α_n which depends on . They first define the conditional probability where equations (10) and (12) have been used to obtain the last equation. The normalization parameter c of this equation only contains terms that do not depend on A, Q or b. The dependence on the matrices A and Q is shown explicitly, including contributions from the normalization factors of P_cond and P_prior.

However, instead of the inverse conditional probability is needed. It is obtained by using Bayes' theorem once more: Very often the function defined in equation (14) is sharply peaked in α-space at a point . One can then choose very broadly without affecting the α dependence of around the peak. Therefore, close to one has . One may use this peak value as the regularization vector in equation (13). With the condition , equation (14) yields N nonlinear equations for the vector of eigenvalues : for n = 1, …, N.

Note that the regularization method sketched above does not require an a priori restriction of the number of model parameters. Instead, it automatically determines that subspace of the whole model-parameter space in which the data determine the outcome of the fit. In the complementary space the a priori values determine the fit. Strong error correlations between two model parameters indicate that the data do not determine them independently.

More extended versions of this article, which includes applications to some typical EXAFS and magnetic EXAFS examples, can be found in Krappe & Rossner (2004) and Krappe et al. (2014).