Strategy¶
The MC_fit strategy for thermobarometric problems is
tied directly to the Nelder-Mead optimization algorithm
([ONeill1971], [MINIM], [Khan2021], [Duesterhoeft2020],
Slide 7).
Given an initial guess of the inversion parameters and a function that quantifies the
misfit between observed and predicted mineralogy, the Nelder-Mead algorithm iteratively
adjusts the inversion parameters to minimize the misfit.
Randomized Nelder-Mead Search¶
A limitation of the Nelder-Mead
algorithm is that it is a local minimizer, which means that it may not find the global
minimum of the misfit function.
To circumvent this limitation, MC_fit does multiple
Nelder-Mead optimizations from randomly generated initial guesses.
This strategy increases the likelihood that MC_fit will find the global minimum
misfit model, and also provides a means of assessing whether the thermobarometric problem is well-posed.
The characteristic of a well-posed problem is that the misfit
function has a pronounced global minimum as a function of the inversion parameters.
If the misfit function is flat with respect to an inversion parameter,
the problem is ill-posed with respect to that parameter, i.e., the
solution for the parameter is not unique.
MC_fit does thermobarometric inversions in two stages:
central model analysis and perturbation analysis.
Fig. 1 Unfiltered MC_fit results of a synthetic garnet + clinopyroxene problem with O2 treated as an unmeasured component and a true solution at {P = 10000 bar, T = 1000 K}. – left panel central model results projected onto the P-T plane. – middle panels the same results projected onto the P-misfit and T-misfit planes. – right panel best central model (gold-filled circle) perturbation analysis results projected onto the P-T plane. The centroid of the inversion parameter ranges, statistically the average of the initial conditions generated for the Nelder-Mead optimizations of the central model analysis, is indicated by a black diamond (MPP). Clustering of central model results toward the upper limit of the specified P-T range is normally an indication that the parameter ranges should be widened. The ranges were not widened here both because the true solution was known and it was desired to make the problem numerically challenging by having the MPP far from the true solution. Because all Nelder-Mead optimizations during the perturbation analysis initiate from the best central model coordinates {P = 9994 bar, T = 999.8 K}, the scatter of the perturbation results is due to the propagation of analytical and thermodynamic uncertainties. This scatter is quantifies the precision of the best central model misfit, \(\epsilon_{misfit}\), and is reported in the legend. The misfit scatter can be visualized by rotating the perturbation analysis results into the P-misfit or T-misfit planes. The misfit of the perturbed models tends to be greater than the misfit of the best central model, when the perturbations result in one or more models with lower misfit, the best overall model is identified (blue-filled circular symbol). When there is a significant difference between the best central and best overall model misfit it may be of interest to investigate the perturbation that leads to the improved fit.¶
Fig. 2 The filtered central model analysis (left panel) for the
synthetic garnet + clinopyroxene problem
illustrated in the previous figure.
The central model results have been filtered to include only models
with misfit values within \(\pm\epsilon_{\text{misfit}}\) of the best central model
misfit, i.e., [0, 2.12e-4].
The filtered central model covariance ellipse quantifies both the uncertainty of the
inversion parameters attributable to ill-posedness of the inverse problem and the
correlations between the inversion parameters.
The covariance ellipse for the perturbed models (right panel) quantifies the
uncertainty of inversion parameters attributable to analytical and
thermodynamic uncertainty.
From Eq 5 and the uncertainties given in the legends of the left and
right panels, the solution to the inverse
problem is T = 999.8 ± 21.2 K, P = 9994 ± 933 bar.
Filled symbols represent models that fit all the analytical data within its uncertainty,
the fit criterion is dependent on the confidence interval (uncertainty_multiplier).
This figure was generated with MC_fit_plot.
Refer to the caption of the previous figure for additional details.¶
Central Model Analysis¶
During the central model analysis, all observed parameters are set to their
central (most probable) values and starting guesses for the unknown
parameters are randomly generated within their specified limits.
From each random starting guess, or
try, a Nelder-Mead optimization is done
to minimize the misfit between the observed and predicted assemblage.
If the optimization is successful, a summary of the result is written to
the user console and the my_project.out file (example).
The best central model is identified (Fig. 1) as
either the try that yields the minimum misfit or the maximum Bayes score, depending on
the Bayes option.
The scatter of the parameter estimates from these tries
(Fig. 1)
generally does not provide a useful means of assessing uncertainties and/or whether the
problem is well-posed because the precision with which the misfit can be quantified
is not known, i.e., typically many of the tries have misfit values that are statistically
indistinguishable from the best central model.
This problem motivates the perturbation analysis, which establishes the misfit function
precision (\(\epsilon_{\text{misfit}}\)) and parameter uncertainties in light of the
uncertainties in the input analytical and thermodynamic data.
As detailed below, once this precision is known, the central model results can be filtered
to include only those tries for which the misfit is statistically indistinguishable from that
of the best central model.
Perturbation Analysis¶
During the perturbation analysis, the input analytical and thermodynamic data are randomly perturbed within their uncertainties. For each perturbation a Nelder-Mead optimization is done starting from the coordinates of the best central model. The scatter of the results (Fig. 1, left panel) measures both the uncertainty of the inversion parameters for the best central model and the precision with which the misfit can be quantified (e.g., \(\sigma_P\), \(\sigma_T\), and \(\sigma_{\text{misfit}}\), as reported in the legend of the left panel of Fig. 1). These parameter uncertainties are independent of the uncertainties attributable to ill-posedness, i.e., if the misfit function is flat with respect to an unknown parameter, the problem is ill-posed with respect to that parameter and the uncertainty on the inferred value of the parameter is effectively infinite, however, the uncertainty analysis captures only the finite uncertainty on the parameter estimates of a specific model. The component of the uncertainty attributable to ill-posedness is obtained from the central model results by filtering the results to include only those tries with misfit within \(\pm\sigma_{\text{misfit}}\) of the best central model misfit (Fig. 2, left panel). The total uncertainty on an estimated parameter \(p\) is then the quadrature sum
where subscripts \(\text{perturbation}\) and \(\text{central}\) denote, respectively, the uncertainties obtained from the perturbation analysis and the filtered central model analysis.
Example, total uncertainty calculation
For the synthetic garnet + clinopyroxene problem, from the legends of the left and right panels of Fig. 2:
\(\{\sigma_{T,\text{central}} = 17.1\,\text{K},\,\sigma_{T,\text{perturbation}} = 12.5\,\text{K}\}\)
\(\{\sigma_{P,\text{central}} = 811.\,\text{bar},\,\sigma_{P,\text{perturbation}} = 462.\,\text{bar}\}\)
and from Eq 5, the total uncertainties on temperature and pressure are
\(\sigma_{T,\text{total}} = \sqrt{\sigma_{T,\text{perturbation}}^2 + \sigma_{T,\text{central}}^2} = \sqrt{12.5^2 + 17.1^2} = 21.2\,\text{K}\)
and
\(\sigma_{P,\text{total}} = \sqrt{\sigma_{P,\text{perturbation}}^2 + \sigma_{P,\text{central}}^2} = \sqrt{462^2 + 811^2} = 933.\,\text{bar}\)
The filtering and processing of MC_fit results illustrated by
Fig. 1 and Fig. 2 was
done with MC_fit_plot.
