<a id="ch-estimation"></a>

<a id="estimation"></a>

# Estimation

This chapter states what the least-squares fit minimises, how it stops, and how each agreement statistic and esd is computed.

<a id="objective-and-weights"></a>

## Objective and weights

<a id="equation-est-obj"></a>
$$
S(\theta) \;=\; \sum_i w_i \bigl(y_{\mathrm{obs},i}
- y_{\mathrm{calc},i}(\theta)\bigr)^2,
\qquad w_i = 1/\sigma_i^2,
$$

*Source:* [`rietx.optimize.least_squares`](https://github.com/yue-here/rietx/blob/v1.10.0/src/rietx/optimize/least_squares.py#L1)

minimised over the residual rows of [(1.5)](forward-model.html.md#equation-fm-rows) [[Rietveld, 1969](references.html.md#id2)].
Weights come from the data file’s esd column whenever it is present;
Poisson $\sigma = \sqrt{\max(y, 1)}$ is only the fallback for bare counts.

The Jacobian is assembled column by column, taking exact work wherever it
exists: linear background columns, analytic peak-chain columns (everything
flowing through per-peak position, width, mixing and intensity), analytic
site-DOF columns for coordinates and anisotropic ADPs over the frozen operator
subsets, analytic axial columns for the FCJ apertures, and plain forward
differences as the fallback.
Under a differentiable backend the same residual is traced and columns
come from forward-mode autodiff, whose cost scales with the parameter
count [[Nocedal and Wright, 2006](references.html.md#id97)]; every backend is held to per-column agreement
with the analytic Jacobian.

<a id="agreement-statistics"></a>

## Agreement statistics

Defined per Toby [[Toby, 2006](references.html.md#id73)]:

<a id="equation-est-indices"></a>
$$
R_{wp} = \sqrt{\frac{\sum w (y_o - y_c)^2}{\sum w y_o^2}}, \qquad
R_{\exp} = \sqrt{\frac{N - P}{\sum w y_o^2}}, \qquad
\chi^2_{\mathrm{red}} = \frac{\sum w (y_o - y_c)^2}{N - P},
$$

*Source:* [`rietx.optimize.statistics`](https://github.com/yue-here/rietx/blob/v1.10.0/src/rietx/optimize/statistics.py#L1)

Here $y_o$ and $y_c$ are the observed and calculated intensities at a channel,
and $w$ is its weight. $N$ is the number of fitted channels and $P$ the number
of refined parameters, so $N - P$ is the degrees of freedom. Also defined are
$\mathrm{GoF} = \sqrt{\chi^2_{\mathrm{red}}} = R_{wp}/R_{\exp}$ and
the background-subtracted $R_{wp}$ variant Toby recommends when the
background carries much of the raw intensity. The Durbin-Watson statistic
on weighted residuals [[Hill and Flack, 1987](references.html.md#id74)] flags serial correlation
($d \approx 2$ ⇒ uncorrelated).

Across the shipped corrections, $\Delta R_{wp}$ recurs as a poor judge of a
physical improvement ([Intensity corrections](corrections.html.md#ch-corrections), [Reading a paper against its own numbers](method.html.md#ch-method)). These indices
measure agreement, and correctness is a separate question.

<a id="structure-agreement-indices"></a>

## Structure agreement indices

Those indices compare profiles. The two that compare the *structure* are
unweighted sums over reflections [[McCusker *et al.*, 1999](references.html.md#id9)] (their equations 13 and 14):

<a id="equation-est-structure-r"></a>
$$
R_B = \frac{\sum_{hkl} \bigl| I_o - I_c \bigr|}{\sum_{hkl} I_o},
\qquad
R_F = \frac{\sum_{hkl} \bigl| |F_o| - |F_c| \bigr|}{\sum_{hkl} |F_o|},
\qquad I_{hkl} = m\,|F_{hkl}|^2,
$$

*Source:* [`rietx.optimize.statistics.structure_r_factors`](https://github.com/yue-here/rietx/blob/v1.10.0/src/rietx/optimize/statistics.py#L1008)

with $m$ the reflection multiplicity, so $|F| = \sqrt{I/m}$ and $R_F$ is the
index a single-crystal $R$ is comparable with. Powder data measure neither sum
directly. $I_o$ is the observed profile partitioned in proportion to $I_c$, the
same partition [Intensities](intensities.html.md#ch-intensities) performs for Le Bail extraction, evaluated
once on the converged structural model. Two consequences follow from the
definition itself.

Both indices are biased towards the model being tested, because a wrong model
receives the intensity it predicted. The paper introduces them for monitoring a
structure’s improvement, and not for judging one in isolation. Both are also
unweighted, so a reflection the weighted fit barely constrains counts as much as
one that dominates it. No weighted variant is computed here, and a trace phase’s
$R_B$ does not compare with the major phase’s.

<a id="the-effective-observation-count"></a>

## The effective observation count

$N$ in [(8.2)](#equation-est-indices) counts profile steps, and the steps across one peak
are repeated measurements of one number. Only the integrated intensities of
individual reflections are unique observations [[McCusker *et al.*, 1999](references.html.md#id9)], and
overlap reduces even those: two reflections at one $2\theta$ are one
observation, and two that partly overlap lie between one and two, because the
profile shape still says how to split them. Altomare et al.
[[Altomare *et al.*, 1995](references.html.md#id13)] make that a count,

<a id="equation-est-mind"></a>
$$
M_{\mathrm{ind}} = \sum_k \frac{I'_k}{I_k},
\qquad
I'_k = I_k - \int_{\chi_k} |F_k|^2\, G(\Delta 2\theta_k)\, \mathrm{d}(2\theta),
$$

*Source:* [`rietx.optimize.statistics.effective_observations`](https://github.com/yue-here/rietx/blob/v1.10.0/src/rietx/optimize/statistics.py#L1214)

where $G$ is the symmetric profile of [Peak profiles](profiles.html.md#ch-profiles) and $\chi_k$ is the
part of reflection $k$’s own interval ($\pm$2.0 FWHM) on
which some overlapping reflection stands higher. An isolated line contributes 1
and the weaker of an exactly coincident pair contributes 0, so the pair is one
observation. The guideline built on it asks for at least
3.0 and preferably 5.0
effective observations per structural parameter. Those are the atomic ones,
since peak positions rather than intensities pay for the cell, profile and
background terms.

The estimate is not a theorem, and both papers say so: the approach “may not
have a rigorous basis”. Its own $\alpha$ is a case in point. The paper’s check
at $\alpha = 4$ lands 6.5 % lower on average, so the value tabulated at
$\alpha = 2$, reported here for comparability with it, runs a little
generous.

<a id="esds-and-the-berar-lelann-inflation"></a>

## Esds and the Bérar-Lelann inflation

rietx reports two esds for each refined parameter.
`RefinedParameter.stderr` is the least-squares esd multiplied by the
Bérar-Lelann factor $f$ [[Bérar and Lelann, 1991](references.html.md#id72)], which corrects for neighbouring
residuals that share a sign. `RefinedParameter.stderr_uncorrected` is the same
esd without $f$. It is the number GSAS-II and FullProf print, so compare theirs
with it. `Statistics.esd_inflation` reports $f$, to within 1 % on a
single-pattern fit.

<a id="equation-est-cov"></a>
$$
\begin{gathered}
\mathrm{Cov} \;=\; \chi^2_{\mathrm{red}} \cdot (J^\top J)^{-1},
\qquad \mathrm{esd}^{\mathrm{uncorr}}_i = \sqrt{\mathrm{Cov}_{ii}},
\\
\mathrm{esd}_i = f \cdot \mathrm{esd}^{\mathrm{uncorr}}_i,
\qquad f = \sqrt{S'' / S},
\end{gathered}
$$

*Source:* [`rietx.optimize.least_squares.covariance_estimates`](https://github.com/yue-here/rietx/blob/v1.10.0/src/rietx/optimize/least_squares.py#L2918)

where $J$ is the Jacobian of the weighted residuals, $\mathrm{esd}_i$ is
`stderr` and $\mathrm{esd}^{\mathrm{uncorr}}_i$ is `stderr_uncorrected`.
[(8.6)](#equation-est-berar-lelann) defines $S$ and $S''$.

<a id="why-the-uncorrected-esd-is-too-small"></a>

### Why the uncorrected esd is too small

The uncorrected esd treats the residual at each channel as independent of its
neighbours. A Rietveld residual rarely is. Where the model misses a peak’s
shape, the misfit keeps one sign across many consecutive channels. Those
channels repeat one piece of information, and the uncorrected esd counts each
of them as new. So it comes out too small.

It also depends on the step width. Take a scan whose misfit is systematic, so
that $\chi^2_{\mathrm{red}}$ is well above 1. Measure the same specimen again
with the same total counts in steps half as wide. Each channel then holds half
the counts and half the misfit. The weighted misfit summed over a peak is
unchanged, and the number of channels doubles. So $\chi^2_{\mathrm{red}}$
halves, $J^\top J$ is unchanged, and every uncorrected esd shrinks by
$\sqrt 2$. The specimen and the information in the scan are the same. Where the
residual is counting noise alone, $\chi^2_{\mathrm{red}} \approx 1$ at either
step width and the uncorrected esd does not move.

<a id="the-factor"></a>

### The factor

$f$ is the homogeneous correction of Bérar and Lelann’s § IV. On the weighted
residuals $a_i$, with $S = \sum a_i^2$, each point carries the probability
$z_i$ that it is correlated with its predecessor. It enters the sum
quadratically with weight $1 - z_i^2$ and linearly with weight $z_i$:

<a id="equation-est-berar-lelann"></a>
$$
\begin{gathered}
S'' = \sum_i (1 - z_i^2)\, a_i^2 + \sum_{\mathrm{runs}}
      \Bigl(\sum_{i \in \mathrm{run}} z_i a_i\Bigr)^2,
\\
z_i = \frac{\sqrt{2\pi x_i}}{2 + \sqrt{2\pi x_i}},\;
x_i = a_i^2 + a_{i-1}^2 \quad (a_i a_{i-1} > 0),
\end{gathered}
$$

*Source:* [`rietx.optimize.statistics.berar_lelann_factor`](https://github.com/yue-here/rietx/blob/v1.10.0/src/rietx/optimize/statistics.py#L420)

and $z_i = 0$ where the sign changes. A run is a stretch of consecutive
$z_i > 0$. Correlated neighbours do not carry independent information, so they
are added coherently, in proportion to how likely they are to be correlated.
The $z_i$ is the paper’s stated derivation: the $\chi^2$ density for two
degrees of freedom over the sum of the densities for one and two,
$f_2/(f_1 + f_2)$ at $x_i$. That is the $2\pi$ of the paper’s own Fortran and
of Andreev’s restatement [[Andreev, 1994](references.html.md#id82)]. The typeset equation (12) prints
2 instead.

Two properties follow. White residuals do not give 1. Independent unit normals
give $E[S'']/E[S] = 1.269$, a factor of 1.13, so a value near 1.13 is no
evidence of correlation. And $z_i$ reads $a_i^2$ itself, so the factor is not
invariant to the residuals’ scale. It grows with $\chi^2_{\mathrm{red}}$ at a
fixed correlation. The paper intends this, since normalised differences grow
with counting time. Before 1.6 the package summed every same-sign run with
weight 1 instead. That is the paper’s rejected § III product test at level
zero, and it reads 1.51 on white residuals.

Two fits shipped with the package, measured on 2026-10-21, with the esd of
the lattice parameter $a$:

| Fit | $f$ | `stderr` of $a$ | `stderr_uncorrected` of $a$ |
|---|---|---|---|
| Fluorapatite, lab Cu Kα (`examples/guides/cell.py`) | 2.17 | 8.3×10<sup>−5</sup> Å | 3.8×10<sup>−5</sup> Å |
| NAC, APS 11-BM synchrotron (`examples/nac_11bm.py`) | 8.02 | 3.9×10<sup>−5</sup> Å | 4.9×10<sup>−6</sup> Å |

So on a fit with nothing left to model, `stderr` is 13 % above
`stderr_uncorrected`. On lab data it is about twice as large, and on a
fine-stepped synchrotron scan about eight times.

Neither esd includes systematic error. On NIST SRM 660c, a LaB<sub>6</sub>
standard measured on a lab instrument, the refined $a$ is 4.156882(21) Å
against 4.156780 Å recomputed from the certificate. That is 4.9 `stderr`
away, with $f$ applied.

The correlation matrix does not carry $f$. It is the true Pearson matrix, so a
genuinely degenerate pair reports $|\rho| \approx 1$ and the 0.98
high-correlation guard means what it says. The esd of a derived quantity, such
as a bond length or a weight fraction, carries $f$ through [(8.7)](#equation-est-derived).
Divide it by $f$ for the uncorrected number. The exact $f$ is `stderr` over
`stderr_uncorrected` on any row of the same fit. `Statistics.esd_inflation`
measures $f$ again at the values the result reports, while the esds carry the
$f$ of the solve’s last residual. The two agreed to 0.75 % or better on the
single-pattern fits measured, and differed by 4 % on a joint fit with
histogram weights 1:4. Values are quoted with two-significant-figure su’s per
the IUCr convention [[Schwarzenbach *et al.*, 1989](references.html.md#id83)].

<a id="what-other-programs-report"></a>

### What other programs report

| Program | Its esd | The rietx number to compare with |
|---|---|---|
| GSAS-II 5.8.4 | $\chi^2_{\mathrm{red}} \cdot (J^\top J)^{-1}$, with no correlation factor | `stderr_uncorrected` |
| FullProf 8.20 | the least-squares esd, with $f$ printed beside it as SCOR | `stderr_uncorrected`, or its esd × SCOR against `stderr` |
| TOPAS (Technical Reference, Version 9) | from the least-squares matrix, or by bootstrap, with no correlation factor described | `stderr_uncorrected` |
| XND 1.42 | multiplied by Bérar’s factor | `stderr` |

GSAS-II multiplies the inverse normal matrix by GOF² (`GSASIIstrMain.py`,
line 397, in 5.8.4, read on 2026-10-20), and GOF² is $\chi^2_{\mathrm{red}}$.
It computes the Durbin-Watson statistic and applies no correlation factor.

FullProf computes Bérar and Lelann’s correction and leaves its own esds as they
are. Its `.sum` file says why, and gives the factor beside two others. From
FullProf 8.20’s refinement of the PbSO<sub>4</sub> round-robin pattern:

```text
 => Standard deviations concern the PRECISION of parameters and
    represent ACCURACY only if there is no systematic errors
    A better estimate of the accuracy of structural parameters
    is obtained multiplying sigmas by the parameter SCOR
 -> SCOR =   2.2790 { Berar'S FORMULA}
 -> SCOR =   2.6633 {Pawley'S FORMULA}
```

The `.out` file repeats it as “Standard deviations have to be multiplied by:
2.2790”. A FullProf esd times SCOR is the counterpart of rietx’s `stderr`, as
far as the two factors agree. rietx’s [(8.6)](#equation-est-berar-lelann) applied to the
final residuals of that FullProf fit gives 2.37, 4 % above FullProf’s 2.28.
The cause of that difference has not been traced.

The TOPAS Technical Reference, Version 9 (October 2026), computes esds
at the end of a refinement when `do_errors` is set (§ 5.9), by singular value
decomposition, and offers `bootstrap_errors` as the alternative (§ 5.10). It
lists the Durbin-Watson statistic among its criteria of fit (§ 5.19). It does
not mention Bérar and Lelann or any correction for serial correlation.

XND 1.42 applies Bérar’s correction, so its esds compare with `stderr`. The
silicon comparison on the validation page relies on that.

<a id="esds-of-derived-quantities"></a>

## Esds of derived quantities

A bond length, an angle or a weight fraction is a function of the refined
parameters rather than one of them, and its esd is the quadratic form

<a id="equation-est-derived"></a>
$$
\sigma_f^2 \;=\; g^\top \mathrm{Cov}\, g,
\qquad g_i = \frac{\partial f}{\partial \theta_i},
$$

*Source:* [`rietx.model.geometry`](https://github.com/yue-here/rietx/blob/v1.10.0/src/rietx/model/geometry.py#L1)

over the whole covariance of [(8.5)](#equation-est-cov): “the whole correlation matrix, not
just the diagonal elements, should be included in the calculation”
[[McCusker *et al.*, 1999](references.html.md#id9)]. Dropping the off-diagonal terms looks conservative and is
not. Measured across the 88 interatomic distances of an 11-BM NAC structural
refinement, the diagonal-only number runs from 0.86 to 1.41 times the full one,
so it is as often too small as too large. Both are reported, so the difference
is visible. It exists only where the coordinates refine,
since a quantity depending on one free parameter has no off-diagonal term to
drop.

A derived esd is absent rather than zero whenever that form cannot be evaluated
honestly. Four cases qualify:

* no covariance at all (an evaluate-only pass, a replayed history node);
* no free parameter the quantity depends on;
* a $g^\top \mathrm{Cov}\, g$ that reaches zero by cancellation, detected as a
  variance below 1e-12 of its own form’s absolute
  terms. A symmetry-fixed 90° angle has exactly zero variance while its partials
  against $x$, $y$ and $z$ do not vanish, so the quadratic form lands on
  roundoff;
* an angle within 2.6e-03° of 0° or 180°, where the
  angle is a stationary point of the coordinates and the linearisation the
  propagation rests on does not hold at all.

The quantity itself is exact in all four cases. Only its uncertainty is
withheld.

The partials are the restraint derivative chain of [(9.17)](parameterisation.html.md#equation-par-restraint)
evaluated at $\sigma = w = 1$. A geometry row is a restraint row, so the two
cannot drift apart, and the weight scale of [(9.18)](parameterisation.html.md#equation-par-restraint-weight) is kept
out of them for the same reason.

<a id="a-derived-esd-describes-one-basin"></a>

## A derived esd describes one basin

Every esd of [(8.5)](#equation-est-cov) and [(8.7)](#equation-est-derived) is built from the curvature of
χ² at the converged point, so it describes the basin the fit stopped in. A
trace phase’s weight fraction is where that stops being enough. The phase’s
scale trades against its width: broadened far enough, its peaks become a hump
the background shares, and the scale can then grow at almost no cost in χ².
That ridge can hold several basins at one χ², each with ordinary curvature and
a tight esd.

The width profile finds them. It pins one width parameter $p$ of the phase on a
grid, refits every other parameter at each value, and admits a point when

<a id="equation-est-profile"></a>
$$
\chi^2(p) - \chi^2_{\min} \;\le\; c \cdot \chi^2_{\mathrm{red}} \cdot
\chi'^2 / \chi^2,
$$

*Source:* [`rietx.strategy.fraction_profile.profile_fraction`](https://github.com/yue-here/rietx/blob/v1.10.0/src/rietx/strategy/fraction_profile.py#L150)

with $\chi^2$ the data’s own weighted sum, $\chi^2_{\min}$ the lowest found
(the fit’s included), and the reduced χ² and Bérar-Lelann factor those of the
fit. Each pinned refit is the constrained hypothesis of Hamilton’s test with a
single constraint [[Hamilton, 1965](references.html.md#id77)], whose statistic is asymptotically
$\chi^2_1$, and $c$ = 3.84 is its 95 % point. The two
factors put the cut on the
scale of the esds themselves, so on a single quadratic basin the admissible
fractions reproduce $W \pm 1.96\,\sigma_W$. Every admissible point’s fraction
lies inside the 95 % profile interval for $W$, so the range the profile reports
is an inner bound, and a finer grid can only widen it. The profile reports how
far its farthest admissible fraction lies from the fit’s, in units of
$1.96\,\sigma_W$, and no threshold turns that distance into a verdict.

The pinned parameter is the width rather than the scale. At a pinned scale the
refit can reach the hump basin only by broadening the phase below the noise,
and there the stage holds the phase’s structure (a phase the data cannot see
is held rather than bounded), so the scan would stop short of the basin it is
looking for. At a pinned width the rest of the problem is close to linear in
the scale, and each refit has one answer.

<a id="an-interval-for-a-parameter-on-its-bound"></a>

## An interval for a parameter on its bound

A parameter the fit leaves on a bound keeps its value and its esd from
[(8.5)](#equation-est-cov), and is flagged. That esd is no interval. An estimate on the
boundary of the parameter space has no symmetric sampling distribution. When
the true value lies on the boundary, the likelihood-ratio statistic follows the
mixture $\tfrac12\chi^2_0 + \tfrac12\chi^2_1$ [[Self and Liang, 1987](references.html.md#id78)]. Half the time the
unconstrained optimum falls outside the range, and the constrained fit then
sits on the bound with a statistic of zero.

The profile interval [[Venzon and Moolgavkar, 1988](references.html.md#id80)] pins the parameter $p$ at a value,
refits every other free parameter, and keeps each value for which

<a id="equation-est-interval"></a>
$$
\chi^2(p) - \chi^2(\hat p) \;\le\; c \cdot \chi^2_{\mathrm{red}} \cdot S''/S,
\qquad
c = \begin{cases}
  F_1^{-1}(\gamma) & \hat p \text{ interior,} \\
  F_1^{-1}(2\gamma - 1) & \hat p \text{ on its bound,}
\end{cases}
$$

*Source:* [`rietx.strategy.parameter_interval.interval_quantile`](https://github.com/yue-here/rietx/blob/v1.10.0/src/rietx/strategy/parameter_interval.py#L86)

with $F_1$ the $\chi^2_1$ distribution function, $\gamma$ the confidence level,
$\hat p$ the fitted value, and $\chi^2_{\mathrm{red}}$ and $S''/S$ those of the
fit. Each pinned refit is Hamilton’s constrained hypothesis with one constraint
[[Hamilton, 1965](references.html.md#id77)], as in [(8.8)](#equation-est-profile). In the interior the statistic is
asymptotically $\chi^2_1$, and $c$ = 3.84 at 95 %. On the bound the mixture
gives $P(T > c) = \tfrac12 P(\chi^2_1 > c)$, so $c$ = 2.71 at 95 %, and the
interval is one-sided with the bound as one end. The two scaling factors put a
two-sided interval on a quadratic profile at the value ± 1.96 esd. Each end is
bracketed by steps of 1, 2, 4, … esd and solved by Brent’s method on the
refitted $\chi^2$.

A bound the data press against is not a stationary point of $\chi^2$. There
$\chi^2$ rises from the bound with a slope as well as a curvature, and the
one-sided end lies inside $\sqrt{c}$ esd. On the BT-1 neutron fit in the test
suite the Lorentzian $Y$ refines to 0 ± 0.0181. Its 95 % end lies at 0.0166,
0.56 of the $\sqrt{2.71}$ esd a stationary bound would give.

<a id="a-moment-angle-at-a-direction-the-symmetry-fixes"></a>

## A moment angle at a direction the symmetry fixes

A moment refines as a modulus and one or two angles. At a direction the site
symmetry fixes, the pattern is even in each angle, so the angle’s Jacobian
column vanishes there while its curvature does not. The covariance of
[(8.5)](#equation-est-cov) scales every column to unit norm before the pseudo-inverse. Beside
such a direction the scaled column points wherever the solver happened to stop,
and so do the esds of everything correlated with it. On a simulated two-site
$Pnma$ pattern with both moments along $\mathbf{a}$, the two moduli read 0.031
and 0.159 μ_B with the angle columns in, and 0.0015 and 0.0011 μ_B with them
left out.

So after each solve the fit tests its free angles $S$. These are the $k$ free
columns that move only moment angles, where some free column also moves each of
those sites’ moduli. A column is read by what it moves, so a named variable
driving two sites’ azimuths is one angle of $S$. The fit moves them, with every
other parameter fixed, until each angle column is at most
1e-06 of the largest of its sites’ modulus columns.
The move is Gauss-Newton on the stacked columns, with first and second
derivatives by central differences, and it ends at $\theta_N$. With
$J_o$ the solve’s Jacobian at the stop $\hat\theta$ without the columns of $S$,
the columns stay in the covariance when

<a id="equation-est-stationary-angle"></a>
$$
\begin{gathered}
\frac{\min_\beta \lVert r(\theta_N) + J_o \beta \rVert^2
      - \min_\beta \lVert r(\hat\theta) + J_o \beta \rVert^2}
     {\chi^2_{\mathrm{red}} \cdot S''/S} \;>\; c_k,
\\
\sum_{j=1}^{k} \binom{k}{j}\, 2^{-k} \bigl(1 - F_j(c_k)\bigr) = 1 - \gamma,
\end{gathered}
$$

*Source:* [`rietx.optimize.least_squares.stationary_angle_test`](https://github.com/yue-here/rietx/blob/v1.10.0/src/rietx/optimize/least_squares.py#L801)

with $r$ the weighted data residual, $F_j$ the $\chi^2_j$ distribution
function and $\gamma$ = 0.95. The two minimisations
refit every other parameter by one linear step at each end. Measured on 96
simulated fits, that step gave the same verdict as refitting in full. The
pattern is even in each angle, so the direction is the boundary of the
parameter space in the squared distance from it. The statistic then follows
the $\bar\chi^2$ mixture of $k$ one-sided directions [[Shapiro, 1985](references.html.md#id79)], and
for one angle that is the $\tfrac12\chi^2_0 + \tfrac12\chi^2_1$ of
[[Self and Liang, 1987](references.html.md#id78)], so $c_1$ is the bound quantile of [(8.9)](#equation-est-interval). At 95 %
$c_k$ is 2.71, 4.23 and 6.50 for one, two and four angles.

*Source:* [`rietx.optimize.statistics.chi_bar_squared_quantile`](https://github.com/yue-here/rietx/blob/v1.10.0/src/rietx/optimize/statistics.py#L483)

A set within its cut leaves the covariance. Its angles keep their values and
have no esd, and the fit says so with `MOMENT_ANGLE_STATIONARY`. A set that
fails drops one angle and is tested again, at most $k$ times. The angle dropped
is the one with the largest second-order change
$\tfrac12 \lVert \partial^2 r / \partial \alpha^2 \rVert\, \delta^2$, or the
one with the largest column where the move did not reach the floor. A moment
tilted 0.05 rad from $\mathbf{a}$ scores about 2500 with the four angles
together, so its polar angle keeps its column. The two azimuths beside it pass
on their own, and leaving them out gives the moduli 0.0015 and 0.0011 μ_B
again.

An angle left out has no esd, but the data still determine it through the
curvature. So a single fit gives each such angle the two-sided interval of
[(8.9)](#equation-est-interval) and stores it on the angle’s row. Its steps start from the
distance either side at which $\chi^2$, with every other parameter held, rises
by $\chi^2_{\mathrm{red}} \cdot S''/S$. On the one-site set both angles’
intervals include the direction, after 24 refits each.

*Source:* [`rietx.strategy.parameter_interval.attach_stationary_intervals`](https://github.com/yue-here/rietx/blob/v1.10.0/src/rietx/strategy/parameter_interval.py#L365)

<a id="a-scale-and-a-b-the-range-cannot-separate"></a>

## A scale and a B the range cannot separate

A uniform displacement parameter $B$ multiplies every reflection of a phase by
$\exp(-2Bs^2)$, beside its scale $S$ [[Prince, 2004](references.html.md#id19)], so

<a id="equation-est-scale-b"></a>
$$
\ln I_{hkl} \;=\; \ln S \;-\; 2B\,s_{hkl}^2 \;+\; \ln K_{hkl},
$$

*Source:* [`rietx.refine._scale_b_separation`](https://github.com/yue-here/rietx/blob/v1.10.0/src/rietx/refine.py#L750)

where $K_{hkl}$ holds every term neither parameter touches. If all of a
phase’s reflections in the fitted range sit at one $s^2$, the data measure only
$\ln S - 2Bs^2$. The two Jacobian columns are then one direction. The
covariance of [(8.5)](#equation-est-cov) does not report that direction as unmeasured. It
equilibrates the normal matrix and passes it to a pseudo-inverse, which
discards every eigenvalue below 1e-15 of the largest. A discarded
direction comes back with zero variance. The scale therefore carries a tight
esd, conditional on a $B$ the data never saw. On a 25–50° Cu Kα scan of a
four-phase mixture, bcc Fe has one reflection in range, (110). It walked to
$B = -150$ Å² and was reported at 0.000 ± 0.000 wt%.

One $s^2$ is the smallest case of a wider one. A displacement parameter
changes peak heights and never peak shapes, so each displacement column lies in
the span of the phase’s reflection profiles in range, and so does the scale’s.
A phase with more displacement parameters beside its scale than it has
reflections in range has a combination of them that imitates the scale.
Fluorite on 25–33° Cu Kα has (111) and (200) against $B_\mathrm{Ca}$ and
$B_\mathrm{F}$. Unbounded, a fit walked them to +99 and −111 Å² and the phase to
0.0 wt% (weighed 1.36).

So at the values each stage starts from, the stage measures a separation:

<a id="equation-est-scale-b-separation"></a>
$$
r \;=\; \min_{b \in \operatorname{span}\{b_j\}}
\left\lVert \hat a - b \right\rVert
\;<\; 2\sqrt{\mathrm{rcond}},
$$

*Source:* [`rietx.refine.SCALE_B_SEPARATION_FLOOR`](https://github.com/yue-here/rietx/blob/v1.10.0/src/rietx/refine.py#L685)

Here $\hat a$ is the phase’s weighted component $y_p/\sigma$ at unit length,
and the $b_j$ are its changes under a step in each free displacement column.
For one column, $r$ is the sine of the angle between the two. Where $r$ falls
under the floor, 6.3e-08, the stage holds the phase’s
displacement parameters. The floor is the cut itself, not a tuned number. The
scale’s column and its projection onto the span, both at unit length, form the
block $[[1, c], [c, 1]]$ with $c = \cos\varphi$. Its small eigenvalue,
$1 - |c| \approx r^2/2$, is discarded against $\lambda_{\max} \approx 2$
exactly when $r < 2\sqrt{\mathrm{rcond}}$. A joint fit has one scale per
histogram and shares $B$, so it reads the smallest angle between the span of
the scales and the span of the $b_j$, stacked over histograms.

With $B$ held, the phase’s weight fraction is conditional on it. At fixed
intensity, an error $\delta B$ moves the fraction by the factor
$\exp(2\,\delta B\,\bar{s}^2)$, which `SCALE_B_INSEPARABLE` quotes per Å². A
pair above the floor is not held, and its esd is honest, however large.

<a id="a-combination-of-three-or-more"></a>

### A combination of three or more

The same cut discards a direction made of more than two columns, and then no
pairwise correlation shows it. Let $E = D\,J^{\mathsf T}J\,D$ be the
equilibrated normal matrix with eigenpairs $(\lambda_k, v_k)$, and let $G$ be
the discarded set, $|\lambda_k| \le \mathrm{rcond}\cdot\lambda_{\max}$. Column
$i$ is *touched* when the discarded subspace adds at least as much variance as
the solve reports:

<a id="equation-est-discarded-touch"></a>
$$
\frac{1}{\mathrm{rcond}\,\lambda_{\max}}\sum_{k \in G} v_{ik}^2
\;\ge\; \sum_{k \notin G} \frac{v_{ik}^2}{|\lambda_k|}.
$$

*Source:* [`rietx.optimize.statistics.discarded_directions`](https://github.com/yue-here/rietx/blob/v1.10.0/src/rietx/optimize/statistics.py#L221)

A discarded eigenvalue is at most the cut, so the left side is a lower bound on
the variance the solve left out. Where it reaches the right side, the esd is at
least $\sqrt{2}$ short. The sums run over the subspace, never one eigenvector,
because two discarded directions are defined only up to a rotation inside it.
Jacobi scaling makes the test independent of a column’s units.

The magnitude on the right is the one the solve itself inverts. A normal matrix
has no negative eigenvalue, but `eigh` returns an exactly flat direction at a
rounding-level eigenvalue of either sign. Inverted with its sign, a negative one
gives every column on that direction a negative variance, and the esd floors to
0. So the pseudo-inverse inverts each kept eigenvalue by its magnitude. On the
LaB₆ + cBN fit with six Lorentzian widths free, the two flat directions came
back between 1e-16 and 4e-15 of the largest. Before the change Linux raised no
`HIGH_CORRELATION` there. On macOS a 1e-14 change to one start value left six
widths at ± 0.

`COVARIANCE_DIRECTION_DISCARDED` names the touched parameters and leaves their
esds as they are. A touched set that is exactly a pair `FLAT_DIRECTION` already
reported adds no row. On the acceptance fixtures 183 of 1 899 whole-pattern
solves discarded such a direction, most of them within a factor of ten of the
cut, over four to eight columns each. Marking those esds absent would blank the
QPA, SRM 676a and brucite fixtures on a threshold that rounding sets, so the
finding reports and the esd stays.

<a id="staged-strategy-and-series"></a>

## Staged strategy and series

Parameter groups are freed cumulatively in the IUCr-guideline order
[[McCusker *et al.*, 1999](references.html.md#id9)]: scale and background first, then peak positions, then
profile widths. The discrete model state is regenerated between stages and
frozen within them. A series of related patterns (an in-situ ramp, a parametric
sweep) is chained by warm starts, and the result is a parameter trajectory,
path-dependent by construction. The chain can be run in both directions, and
parameters the two runs disagree on are flagged. That is the one check
separating a measured trajectory from an ordering artefact. Parametric
refinement [[Stinton and Evans, 2007](references.html.md#id14)] is the other route: every pattern in one joint
residual, with a parameter held to a declared function of each pattern’s
coordinates ([A parameter as a function of each pattern’s coordinates](parameterisation.html.md#par-parametric-section)).

<a id="solvers"></a>

## Solvers

The default driver is a bounded Levenberg-Marquardt. scipy’s Trust Region
Reflective, the previous default, stays available as `solver="trf"`. The
Levenberg-Marquardt driver implements Coelho’s adaptive Marquardt
constant [[Coelho, 2018](references.html.md#id87)], with the system diagonally pre-conditioned. That
pre-conditioning makes λ dimensionless and lets the published constants
transfer. Each damped step solves

<a id="equation-est-lm-step"></a>
$$
\begin{aligned}
\Delta \;=\; &\mathop{\mathrm{arg\,min}}_{\Delta}\;
  \tfrac12\,\Delta^\top \bigl(A + \lambda D^2\bigr)\,\Delta
  - b^\top\Delta \\
&\text{subject to}\quad \ell \le \theta + \Delta \le u
  \quad\text{and}\quad T(\theta + \Delta) + c \ge 0,
\end{aligned}
$$

*Source:* [`rietx.optimize.lm._solve_constrained_step`](https://github.com/yue-here/rietx/blob/v1.10.0/src/rietx/optimize/lm.py#L779)

with $\ell$ and $u$ the declared limits and the rows $T\theta + c \ge 0$
described below. $D$ is diagonal, and $D_{ii}$ is the largest norm column $i$
of $J$ has had in the stage, the scaling of MINPACK [[Moré, 1978](references.html.md#id93)], or the
square root of a curvature held for the parameter where that is larger (below).
Without
a row binding, the step is solved by bounded-variable least squares
[[Stark and Parker, 1995](references.html.md#id84)], so a parameter the step puts on a bound lands on it. Above
128 parameters, a size only a Pawley block reaches, that solve grows too slow.
There the step is projected conjugate gradients from the Cauchy point
[[Lin and Moré, 1999](references.html.md#id95); [Moré and Toraldo, 1991](references.html.md#id94)], every segment of which lowers the model. The step
there until 1.9 was the bound-constrained conjugate gradient of
[[Coelho, 2005](references.html.md#id86)] (conjugate gradients per [[Hestenes and Stiefel, 1952](references.html.md#id90); [Polak, 1971](references.html.md#id91)]).
On the Pawley fits measured, it raised the model at most damped steps and
stopped the fit on its evaluation cap. Conventions: $A = J^\top J$,
$b = -J^\top r$, and the paper’s objective $S = r^\top r$ is $\chi^2$.

$D$ keeps each column’s largest norm because a column can shrink towards a
point where it vanishes, as a moment angle does at a direction the symmetry
fixes. Scaled to unit norm there, that column took the longest step, and λ
rose to hold it until every parameter crawled. A column that is small from the
start of a stage has no larger norm to keep [[Transtrum and Sethna, 2012](references.html.md#id96)]. Its
Gauss-Newton curvature $A_{ii}$ can lie orders of magnitude under the
curvature of $S$: $2.2\times10^{-5}$ against $3.4\times10^{5}$ for a moment
angle $3\times10^{-8}$ rad from such a direction. So once rejected steps have
raised λ to 100, the driver measures $\tfrac12\,\partial^2 S/\partial\theta_i^2$
by a central second difference along the parameter the last rejected step
moved furthest for its size. Where that curvature is positive and larger than
$A_{ii} + \lambda D_{ii}^2$, it replaces $A_{ii}$ for the rest of the stage, and
λ starts again. Each parameter is measured at most once a stage this way. A
held curvature is measured again before the stage may stop, and dropped if it
has turned negative. Held at a maximum, it would make the maximum read as a
minimum.

It is the default because it settles where TRF crawls. Every parameter refines
in its own units, with its declared limits as box bounds, and TRF’s strictly
feasible iterates approach a bound slowly. On 14 acceptance fits each started
five ways, one fit’s χ² spread by more than 1e-9 under this driver and four
under TRF. From a poor start TRF can also end in a different basin. Issue
#900’s synthetic PbSO₄ frees every parameter in one stage, from a scale 2400
times too high. This driver reaches χ² 7111 with the Pb $B_{iso}$ at 1.48 Å²,
against a true 1.49. TRF stops on its iteration limit at χ² 4.5 × 10⁵, with
the Pb $B_{iso}$ at 133 Å². It reports `STAGE_MAX_ITER` and
`BISO_UNUSUALLY_LARGE`.

It also carries constraints TRF cannot. Linear inequalities on functionals of
θ are rows $T\theta + c \ge 0$, a shape no per-parameter box can express. Two
kinds exist. One is the Stephens positivity cone ([Microstructure](microstructure.html.md#ch-microstructure)).
The other is the floor on the Gaussian width at every reflection,
[(3.3)](profiles.html.md#equation-prof-caglioti-floor). When the bounded step crosses a row, the step is
solved again over the bounds and the rows together. That is a least-squares
problem with linear inequality constraints. It is solved exactly by Lawson and
Hanson’s reduction to a least-distance problem and then to non-negative least
squares [[Lawson and Hanson, 1974](references.html.md#id85)]. A step that crosses no row is the bounded step
unchanged. The rule before truncated the bounded step at the first row it
crossed, and that stalled wherever the optimum lay on a row’s face. An answer that
pressed a row says so. Per-stage counts of the held steps are recorded, and a
`CONSTRAINT_ACTIVE` diagnostic names the block in the answer-producing stage.
It is the one signal that a declared constraint was active.

Speed was measured with `examples/bench_solver.py` on four acceptance
protocols (macOS arm64, three runs). Where both drivers converge, on SRM 660c
and the 11-BM NAC pattern, this driver is 1.19–1.30× faster. It needs 78 and 35
residual evaluations there, against TRF’s 83 and 51. The solve itself is a
minority of the runtime, so a faster solve alone would buy at most about
1.25×. On corundum and the four-phase cpd-2 protocol TRF stops on its
iteration limit after 906 and 3865 evaluations. This driver converges in 116
and 224, 11.6–35× faster. [Reading a paper against its own numbers](method.html.md#ch-method) works through two places where
the Coelho papers disagree with their own text, and how measurement settled
them.

<a id="convergence"></a>

## Convergence

The Rietveld literature judges convergence by the final cycle’s parameter
shifts against their own esds, rather than by a cost decrement
[[McCusker *et al.*, 1999](references.html.md#id9)] (their §7):

<a id="equation-est-convergence"></a>
$$
\max_i \; \frac{|\Delta \theta_i|}{\mathrm{esd}(\theta_i)}
\;\le\; \epsilon_{\mathrm{conv}},
$$

*Source:* [`rietx.optimize.least_squares.run_least_squares`](https://github.com/yue-here/rietx/blob/v1.10.0/src/rietx/optimize/least_squares.py#L2516)

with the paper’s band $\epsilon_{\mathrm{conv}} =$ 0.1
quoted, never tuned. Both sides are measured in external parameter units.
$\Delta\theta$ is decoded exactly through the transform chain of
[Parameterisation and constraints](parameterisation.html.md#ch-parameterisation), and the esd is the chain-ruled physical one of
[(8.5)](#equation-est-cov), because at finite step size an internal-space ratio is a
different number.

The measured value is `Statistics.max_shift_over_esd`, computed by the solver
from the final accepted step and copied onto the answer-producing stage’s
statistics. Nothing else derives it, and it gates nothing. A converged Trust
Region Reflective solve at `ftol` $10^{-9}$ satisfies the criterion with a wide margin
(measured $\sim 3\times 10^{-4}$ on the synthetic LaB₆ round trip), so the
information is on the other branch: a stage stopped on its iteration budget
reports how far it was still moving in esd units, and the same stage starved to
one iteration measures $\approx 14$. The value is withheld wherever it cannot
be measured: with no accepted step, no esds, or an evaluate-only replay.

A joint multi-histogram fit has the same quantity over its joint $\theta$. Each
column has one shift, the joint solve’s last accepted step. It has one esd,
from the covariance of the weighted stacked residual. The histogram weights
enter both, so they change the number and leave the criterion as it is. A
shared column appears in every histogram’s table with the same ratio, because
each table decodes it affinely.

<a id="the-fp64-floor"></a>

## The fp64 floor

The residual used for cost and statistics, and the parameter solve and
covariance, are always fp64 on host. A GPU backend may compute Jacobian columns
in fp32. Three properties set that asymmetry.

- The residual cancels. $\sqrt{w}(y_{\mathrm{obs}} - y_{\mathrm{calc}})$
  subtracts numbers of order 10⁵ counts to leave order 10², and fp32’s ~7
  digits put an absolute error of order 10 counts into everything that reads it,
  which is ~10 % of that difference.
- The solve squares the conditioning [[Higham, 2002](references.html.md#id92)]:

<a id="equation-est-cond"></a>
$$
\operatorname{cond}(J^\top J) \;=\; \operatorname{cond}(J)^2,
$$

*Source:* [`rietx.backend.linalg64`](https://github.com/yue-here/rietx/blob/v1.10.0/src/rietx/backend/linalg64.py#L1)

so a routine Rietveld $\operatorname{cond}(J) \sim 10^4$ leaves
$\operatorname{cond}(J^\top J) \sim 10^8$, which fp32 cannot invert at all.
The bounded LM forms $J^\top J$ explicitly, and is the direct illustration.

- Columns are relative-accuracy tolerant. A column enters through a descent
  direction and a curvature estimate, and the trust region re-measures every
  step against a fresh fp64 cost. Measured on real hardware, an Apple-GPU
  refinement with every column in fp32 lands 3.5×10⁻⁸ Å from the numpy fp64
  cell.

<a id="from-fit-to-report"></a>

## From fit to report

The FitReport reads the converged state in three layers. Layer 0 is model-free:
cumulative-χ² breakpoints localise where misfit lives [[David, 2004](references.html.md#id81)], and
unindexed peaks are flagged against the tick positions of every emission line
(an impurity must clear 8.0σ). Layer 1 attributes per-region
misfit to profile shape derivatives. Each region passes four gates: local-χ²
significance, a local R² showing the shape basis explains the region’s misfit,
resolvability on the scale-normalised Gram, and a validity radius. A global
maturity gate stands above them, so a collinear pair is declared non-separable
rather than resolved into a confident wrong singleton. Layer 2 turns attributions into actions. The
actions are advisory and carry no statistical test of their own. Whether a
freed parameter pays for itself is ΔBIC [[Schwarz, 1978](references.html.md#id76)], predicted before
the fit by `Refinement.suggest` and measured after it by
`report.compare_freed` (next section). The report’s thresholds are versioned
(currently 1.13).

<a id="which-parameter-to-free-next"></a>

## Which parameter to free next

`Refinement.suggest()` ranks every held-but-refinable parameter by the χ²
reduction one Gauss-Newton solve would obtain from freeing it, at the cost
of a single Jacobian evaluation and no solve. With $F$ the currently-free
columns, $P_F$ the orthogonal projector onto their span, $r$ the weighted
residual and $J_j$ a held parameter’s column,

<a id="equation-est-suggest"></a>
$$
\tilde{\jmath} = (I - P_F)\, J_j, \qquad
\tilde{r} = (I - P_F)\, r, \qquad
\Delta\chi^2_j \;=\; \frac{(\tilde{\jmath}^{\top} \tilde{r})^2}
                          {\tilde{\jmath}^{\top} \tilde{\jmath}},
$$

*Source:* [`rietx.optimize.statistics.one_parameter_gains`](https://github.com/yue-here/rietx/blob/v1.10.0/src/rietx/optimize/statistics.py#L798)

which is Rao’s score statistic [[Rao, 1948](references.html.md#id123)] applied to the linearised model,
computed through the Frisch-Waugh-Lovell projection identity
[[Frisch and Waugh, 1933](references.html.md#id124); [Lovell, 1963](references.html.md#id125)]. It is exactly the drop in $\sum w\Delta^2$ that
a least-squares solve of $[F \mid J_j]$ achieves over $F$ alone. It is invariant
under any rescaling of the column, so no per-parameter step heuristics are
needed. At a converged minimum $J^\top r \approx 0$ makes every gain vanish.
GSAS-II answers the same recipe problem [[Toby, 2024](references.html.md#id122)] by ±δ finite
differences with per-type δ heuristics and a sign-consistency test, because its
analytic derivatives are locked inside Hessian assembly. Exact columns at the
current state make all three workarounds unnecessary.

Under the null hypothesis a gain is distributed as
$\chi^2_1 \cdot \chi^2_{\mathrm{red}}$, so a candidate is quotable only above
a noise floor of 9.0 · max(χ²_red, 1). That is the 3σ point
of $\chi^2_1$, with the same floor-at-one convention as the covariance scale.
Two gates keep the ranking honest, the Layer-1 discipline one call over. The
free block absorbs some candidates’ columns, and those are reported
non-separable and carry no score; the same projection caps the $1/(1-R^2)$
inflation of near-collinear gains. Candidates whose projected columns are
pairwise indistinguishable come back as one unresolved group carrying a joint
gain, which is a tie. As with indexing there is no `.best`, and
`best_or_none()` answers `None` whenever the evidence does not choose one
parameter.

A gain that clears the floor is then read as a model-selection answer, ΔBIC
[[Schwarz, 1978](references.html.md#id76)], and that needs a premise the gain does not: Schwarz’s
$N$ counts independent observations. A powder residual is serially
correlated, and at raw channel counts the reward term
$N\ln(\chi^2_r/\chi^2_f)$ outvotes the $\ln N$ penalty for almost any
improvement. The fit already measures the correlation, as the Bérar-Lelann
factor $f$ every esd is inflated by [[Bérar and Lelann, 1991](references.html.md#id72)], and the penalty is
charged at

<a id="equation-est-effective-n"></a>
$$
N_{\mathrm{eff}} = N / f^{2}, \qquad
\Delta\mathrm{BIC} = N_{\mathrm{eff}} \ln\frac{\chi^2_r}{\chi^2_f}
                     - k \ln N_{\mathrm{eff}},
$$

*Source:* [`rietx.optimize.statistics.effective_sample_size`](https://github.com/yue-here/rietx/blob/v1.10.0/src/rietx/optimize/statistics.py#L540)

the count at which, for one added parameter, the reward term is that
parameter’s $t^2$ at its inflated esd, so ΔBIC is positive exactly when
$t^2 > \ln N_{\mathrm{eff}}$. It is a heuristic and not a theorem, and it is
conservative the way the esds are: white residuals still give $f \approx
1.13$. Hamilton’s test [[Hamilton, 1965](references.html.md#id77)] takes the same $N_{\mathrm{eff}}$,
so the two verdicts are read off one count. Measured with the pre-1.6 factor on
four ~49 500-channel synchrotron fits of one occupancy, raw $N$ gave ΔBIC +36 to +211 to a
parameter each fit’s own esd put within 0.76-1.89σ of zero, and
$N_{\mathrm{eff}}$ turned all four negative.

After the fit the same count prices the measured gain. `report.compare_freed`
takes the fit without a parameter and the fit with it, and returns
[(8.18)](#equation-est-effective-n) beside each freed parameter’s $t$. The two agree for one
parameter, so their disagreement is information: a multi-parameter block whose
joint ΔBIC is positive while no member’s $t$ clears 2, or an $f$ that moved
between the fits.