Skip to content

Commit c448b70

Browse files
Correct weighted R-factor calculations and document reliability metrics (#224)
* Fix weighted R-factor uncertainty handling * Document reliability factor definitions * Add pending ICSD CIF import regression test * Simplify LaM7O3 tutorial sidebar label * Enable ICSD CIF import regression test
1 parent c4ae7e8 commit c448b70

8 files changed

Lines changed: 345 additions & 17 deletions

File tree

docs/docs/user-guide/analysis-workflow/analysis.md

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -265,6 +265,9 @@ Now, you can inspect the fitted parameters to see how they have changed
265265
during the refinement process, select more parameters to be refined, and
266266
perform additional fits as needed.
267267

268+
The equations and interpretation of the reported fit-quality values are
269+
listed on the [Reliability Factors](reliability-factors.md) page.
270+
268271
To plot the measured and calculated data after the fit, you can use the
269272
`pattern` method of the `display` object:
270273

Lines changed: 156 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,156 @@
1+
---
2+
title: Reliability Factors
3+
icon: material/chart-bell-curve-cumulative
4+
---
5+
6+
# :material-chart-bell-curve-cumulative: Reliability Factors
7+
8+
EasyDiffraction reports several complementary measures of agreement
9+
between observed and calculated diffraction data. They are ratios
10+
internally; values shown with a percent sign are multiplied by 100.
11+
12+
For the equations below, $y_i^{\mathrm{obs}}$ and $y_i^{\mathrm{calc}}$
13+
are the observed and calculated intensities, $\sigma_i$ is the standard
14+
uncertainty of the observed intensity, and
15+
16+
$$
17+
w_i = \frac{1}{\sigma_i^2}
18+
$$
19+
20+
is its inverse-variance weight. The sum runs over the $N$ data points
21+
included in the reported value.
22+
23+
## R-factor (Rf)
24+
25+
The unweighted profile R-factor is the absolute difference between
26+
observed and calculated intensities, normalized by the total absolute
27+
observed intensity:
28+
29+
$$
30+
R_f =
31+
\frac{\sum_i \left|y_i^{\mathrm{obs}}-y_i^{\mathrm{calc}}\right|}
32+
{\sum_i \left|y_i^{\mathrm{obs}}\right|}.
33+
$$
34+
35+
Lower values indicate closer agreement. This factor does not use the
36+
measurement uncertainties, so every data point contributes according to
37+
the magnitude of its absolute residual.
38+
39+
## Squared-residual R-factor (Rf²)
40+
41+
The value labelled `Rf²` in the fit summary is the unweighted
42+
root-squared-residual ratio:
43+
44+
$$
45+
R_{f^2} =
46+
\left[
47+
\frac{\sum_i \left(y_i^{\mathrm{obs}}-y_i^{\mathrm{calc}}\right)^2}
48+
{\sum_i \left(y_i^{\mathrm{obs}}\right)^2}
49+
\right]^{1/2}.
50+
$$
51+
52+
Despite the compact `Rf²` label, this is not the algebraic square of
53+
$R_f$. The label indicates that squared intensities and residuals are
54+
used before taking the square root.
55+
56+
## Weighted R-factor (wR)
57+
58+
The weighted R-factor is the root ratio of weighted squared residuals:
59+
60+
$$
61+
wR =
62+
\left[
63+
\frac{\sum_i w_i
64+
\left(y_i^{\mathrm{obs}}-y_i^{\mathrm{calc}}\right)^2}
65+
{\sum_i w_i \left(y_i^{\mathrm{obs}}\right)^2}
66+
\right]^{1/2},
67+
\qquad
68+
w_i = \frac{1}{\sigma_i^2}.
69+
$$
70+
71+
Consequently, a point with a smaller standard uncertainty has more
72+
influence than a less precise point. EasyDiffraction expects standard
73+
uncertainties as input and converts them to inverse-variance weights; it
74+
does not use $\sigma_i$ itself as the weight.
75+
76+
## Chi-square and reduced chi-square
77+
78+
The uncertainty-weighted sum of squared residuals is
79+
80+
$$
81+
\chi^2 =
82+
\sum_i \left(
83+
\frac{y_i^{\mathrm{obs}}-y_i^{\mathrm{calc}}}{\sigma_i}
84+
\right)^2
85+
= \sum_i w_i
86+
\left(y_i^{\mathrm{obs}}-y_i^{\mathrm{calc}}\right)^2.
87+
$$
88+
89+
If $p$ free parameters were fitted, the number of degrees of freedom is
90+
$\nu=N-p$, and the reported goodness-of-fit is the reduced chi-square:
91+
92+
$$
93+
\chi_\nu^2 = \frac{\chi^2}{\nu}.
94+
$$
95+
96+
A value near 1 means that the size of the residuals is consistent with
97+
the stated standard uncertainties. A much larger value can indicate a
98+
poor model or underestimated uncertainties; a much smaller value can
99+
indicate overestimated uncertainties or an over-flexible model.
100+
101+
For a joint fit, EasyDiffraction also multiplies each experiment's
102+
squared normalized residuals by its normalized joint-fit weight. Those
103+
experiment weights are normalized so that their sum equals the number of
104+
experiments.
105+
106+
## Expected weighted profile R-factor
107+
108+
For powder fits, the expected weighted profile R-factor is
109+
110+
$$
111+
wR_{\mathrm{expected}} =
112+
\left[
113+
\frac{\nu}
114+
{\sum_i w_i \left(y_i^{\mathrm{obs}}\right)^2}
115+
\right]^{1/2}.
116+
$$
117+
118+
It is the weighted profile R-factor expected when $\chi_\nu^2=1$.
119+
Therefore, $wR / wR_{\mathrm{expected}} = \sqrt{\chi_\nu^2}$ when the
120+
same data points and weights are used for both values.
121+
122+
## Bragg R-factor (BR)
123+
124+
When observed and calculated structure-factor magnitudes are available,
125+
EasyDiffraction can report the Bragg R-factor:
126+
127+
$$
128+
BR =
129+
\frac{\sum_h \left|F_h^{\mathrm{obs}}-F_h^{\mathrm{calc}}\right|}
130+
{\sum_h F_h^{\mathrm{obs}}}.
131+
$$
132+
133+
Here $h$ indexes reflections and $F_h$ is a structure-factor magnitude.
134+
Lower values indicate closer agreement between observed and calculated
135+
reflection amplitudes.
136+
137+
## Names and data subsets
138+
139+
The fit summary uses the short labels `Rf`, `Rf²`, `wR`, and `BR`. Saved
140+
deterministic fit results also expose IUCr-style names:
141+
142+
| Saved result | Definition and scope |
143+
| -------------------- | ---------------------------------------------------------------- |
144+
| `R_factor_all` | $R_f$ for all included observations |
145+
| `wR_factor_all` | $wR$ for all included observations |
146+
| `R_factor_gt` | $R_f$ for observations satisfying $y_i^{\mathrm{obs}}>3\sigma_i$ |
147+
| `wR_factor_gt` | $wR$ for observations satisfying $y_i^{\mathrm{obs}}>3\sigma_i$ |
148+
| `prof_R_factor` | $R_f$ for all included powder-profile points |
149+
| `prof_wR_factor` | $wR$ for all included powder-profile points |
150+
| `prof_wR_expected` | $wR_{\mathrm{expected}}$ for all included powder-profile points |
151+
| `reduced_chi_square` | $\chi_\nu^2$ for the fitted residual vector |
152+
153+
Only finite observations with finite calculated values and positive,
154+
finite standard uncertainties are included in saved deterministic
155+
statistics. A metric is unavailable when its denominator is zero or when
156+
it does not apply to the fitted data.

docs/mkdocs.yml

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -222,6 +222,7 @@ nav:
222222
- Structure: user-guide/analysis-workflow/model.md
223223
- Experiment: user-guide/analysis-workflow/experiment.md
224224
- Analysis: user-guide/analysis-workflow/analysis.md
225+
- Reliability Factors: user-guide/analysis-workflow/reliability-factors.md
225226
- Report: user-guide/analysis-workflow/report.md
226227
- Tutorials:
227228
- Tutorials: tutorials/index.md
@@ -235,7 +236,7 @@ nav:
235236
- HS pd-neut-cwl: tutorials/refine-hs-hrpt.ipynb
236237
- Si pd-neut-tof: tutorials/refine-si-sepd.ipynb
237238
- PbSO4 pd-xray-cwl: tutorials/refine-pbso4-xray.ipynb
238-
- LaM(7)O3 P02.1 pd-xray-cwl: tutorials/refine-lam7o3-p021.ipynb
239+
- LaM(7)O3 pd-xray-cwl: tutorials/refine-lam7o3-p021.ipynb
239240
- LMO pd-neut-cwl: tutorials/refine-lmo-echidna.ipynb
240241
- Without Measured Data:
241242
- LBCO pd-neut-cwl: tutorials/simulate-lbco-cwl.ipynb

src/easydiffraction/analysis/analysis.py

Lines changed: 2 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -44,6 +44,7 @@
4444
from easydiffraction.analysis.fit_helpers.bayesian import PosteriorSamples
4545
from easydiffraction.analysis.fit_helpers.bayesian import posterior_predictive_cache_key
4646
from easydiffraction.analysis.fit_helpers.metrics import calculate_r_factor
47+
from easydiffraction.analysis.fit_helpers.metrics import calculate_weighted_r_factor
4748
from easydiffraction.analysis.fit_helpers.reporting import FitResults
4849
from easydiffraction.analysis.fitting import Fitter
4950
from easydiffraction.analysis.fitting import FitterFitOptions
@@ -2035,12 +2036,7 @@ def _weighted_r_factor_or_none(
20352036
"""Return a weighted R factor when inputs are available."""
20362037
if observed.size == 0:
20372038
return None
2038-
weights = 1.0 / uncertainties**2
2039-
denominator = float(np.sum(weights * observed**2))
2040-
if denominator <= 0.0:
2041-
return None
2042-
numerator = float(np.sum(weights * (observed - calculated) ** 2))
2043-
value = np.sqrt(numerator / denominator)
2039+
value = calculate_weighted_r_factor(observed, calculated, uncertainties)
20442040
return float(value) if np.isfinite(value) else None
20452041

20462042
@staticmethod

src/easydiffraction/analysis/fit_helpers/metrics.py

Lines changed: 6 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -44,7 +44,7 @@ def calculate_r_factor(
4444
def calculate_weighted_r_factor(
4545
y_obs: np.ndarray,
4646
y_calc: np.ndarray,
47-
weights: np.ndarray,
47+
standard_uncertainties: np.ndarray,
4848
) -> float:
4949
"""
5050
Calculate weighted R-factor between observed and calculated data.
@@ -55,8 +55,9 @@ def calculate_weighted_r_factor(
5555
Observed data points.
5656
y_calc : np.ndarray
5757
Calculated data points.
58-
weights : np.ndarray
59-
Weights for each data point.
58+
standard_uncertainties : np.ndarray
59+
Standard uncertainties for the observed data points. The
60+
inverse-variance weights are calculated as ``1 / sigma**2``.
6061
6162
Returns
6263
-------
@@ -65,7 +66,8 @@ def calculate_weighted_r_factor(
6566
"""
6667
y_obs = np.asarray(y_obs)
6768
y_calc = np.asarray(y_calc)
68-
weights = np.asarray(weights)
69+
standard_uncertainties = np.asarray(standard_uncertainties)
70+
weights = 1.0 / standard_uncertainties**2
6971
numerator = np.sum(weights * (y_obs - y_calc) ** 2)
7072
denominator = np.sum(weights * y_obs**2)
7173
return np.sqrt(numerator / denominator) if denominator != 0 else np.nan
Lines changed: 151 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,151 @@
1+
# SPDX-FileCopyrightText: 2026 EasyScience contributors <https://github.com/easyscience>
2+
# SPDX-License-Identifier: BSD-3-Clause
3+
"""Regression coverage for importing an underscore-style ICSD CIF."""
4+
5+
import numpy as np
6+
import pytest
7+
8+
from easydiffraction.analysis.calculators.cryspy import CryspyCalculator
9+
from easydiffraction.datablocks.structure.item.factory import StructureFactory
10+
11+
ZRW2O8_CIF = """\
12+
data_83267-ICSD
13+
_database_code_ICSD 83267
14+
_audit_creation_date 1998-06-26
15+
_chemical_name_systematic
16+
'Zirconium bis(tungstate)'
17+
_chemical_formula_structural
18+
'Zr (W O4)2'
19+
_chemical_formula_sum
20+
'O8 W2 Zr1'
21+
_publ_section_title
22+
'O8 W2 Zr1'
23+
loop_
24+
_citation_id
25+
_citation_journal_abbrev
26+
_citation_year
27+
_citation_journal_volume
28+
_citation_page_first
29+
_citation_page_last
30+
_citation_journal_id_ASTM
31+
primary 'Science' 1996 272 90 92 SCIEAS
32+
loop_
33+
_publ_author_name
34+
Mary, T.A.;Evans, J.S.O.;Vogt, T.;Sleight, A.W.
35+
_cell_length_a 9.15993(5)
36+
_cell_length_b 9.15993(5)
37+
_cell_length_c 9.15993(5)
38+
_cell_angle_alpha 90.
39+
_cell_angle_beta 90.
40+
_cell_angle_gamma 90.
41+
_cell_volume 768.56
42+
_cell_formula_units_Z 4
43+
_symmetry_space_group_name_H-M 'P 21 3'
44+
_symmetry_Int_Tables_number 198
45+
_refine_ls_R_factor_all 0.024000
46+
loop_
47+
_symmetry_equiv_pos_site_id
48+
_symmetry_equiv_pos_as_xyz
49+
1 '-z+1/2, -x, y+1/2'
50+
2 '-y+1/2, -z, x+1/2'
51+
3 '-x+1/2, -y, z+1/2'
52+
4 '-z, x+1/2, -y+1/2'
53+
5 '-y, z+1/2, -x+1/2'
54+
6 '-x, y+1/2, -z+1/2'
55+
7 'z+1/2, -x+1/2, -y'
56+
8 'y+1/2, -z+1/2, -x'
57+
9 'x+1/2, -y+1/2, -z'
58+
10 'z, x, y'
59+
11 'y, z, x'
60+
12 'x, y, z'
61+
loop_
62+
_atom_type_symbol
63+
_atom_type_oxidation_number
64+
Zr4+ 4
65+
W6+ 6
66+
O2- -2
67+
loop_
68+
_atom_site_label
69+
_atom_site_type_symbol
70+
_atom_site_symmetry_multiplicity
71+
_atom_site_Wyckoff_symbol
72+
_atom_site_fract_x
73+
_atom_site_fract_y
74+
_atom_site_fract_z
75+
_atom_site_B_iso_or_equiv
76+
_atom_site_occupancy
77+
Zr1 Zr4+ 4 a 0.0003(4) 0.0003(4) 0.0003(4) 0.010(1) 1.
78+
W1 W6+ 4 a 0.3412(3) 0.3412(3) 0.3412(3) 0.012(1) 1.
79+
W2 W6+ 4 a 0.6008(3) 0.6008(3) 0.6008(3) 0.010(1) 1.
80+
O1 O2- 12 b 0.2071(3) 0.4378(4) 0.4470(3) 0.022(1) 1.
81+
O2 O2- 12 b 0.7876(3) 0.5694(4) 0.5565(3) 0.020(1) 1.
82+
O3 O2- 4 a 0.4916(5) 0.4916(5) 0.4916(5) 0.023(1) 1.
83+
O4 O2- 4 a 0.2336(3) 0.2336(3) 0.2336(3) 0.037(1) 1.
84+
"""
85+
86+
87+
def test_icsd_cif_import_preserves_structure_and_ionic_symbols_for_cryspy():
88+
"""Import the ICSD structure and preserve ionic atom types through Cryspy."""
89+
from cryspy.H_functions_global.function_1_cryspy_objects import str_to_globaln
90+
91+
structure = StructureFactory.from_cif_str(ZRW2O8_CIF)
92+
93+
assert structure.name == '83267-icsd'
94+
assert structure.as_cif.startswith('data_83267-icsd\n')
95+
assert structure.space_group.name_h_m.value == 'P 21 3'
96+
np.testing.assert_allclose(
97+
[
98+
structure.cell.length_a.value,
99+
structure.cell.length_b.value,
100+
structure.cell.length_c.value,
101+
structure.cell.angle_alpha.value,
102+
structure.cell.angle_beta.value,
103+
structure.cell.angle_gamma.value,
104+
],
105+
[9.15993, 9.15993, 9.15993, 90.0, 90.0, 90.0],
106+
)
107+
np.testing.assert_allclose(
108+
[
109+
structure.cell.length_a.uncertainty,
110+
structure.cell.length_b.uncertainty,
111+
structure.cell.length_c.uncertainty,
112+
],
113+
[0.00005, 0.00005, 0.00005],
114+
)
115+
116+
expected_sites = {
117+
'Zr1': ('Zr4+', 4, 'a', (0.0003, 0.0003, 0.0003), (0.0004,) * 3, 0.010),
118+
'W1': ('W6+', 4, 'a', (0.3412, 0.3412, 0.3412), (0.0003,) * 3, 0.012),
119+
'W2': ('W6+', 4, 'a', (0.6008, 0.6008, 0.6008), (0.0003,) * 3, 0.010),
120+
'O1': ('O2-', 12, 'b', (0.2071, 0.4378, 0.4470), (0.0003, 0.0004, 0.0003), 0.022),
121+
'O2': ('O2-', 12, 'b', (0.7876, 0.5694, 0.5565), (0.0003, 0.0004, 0.0003), 0.020),
122+
'O3': ('O2-', 4, 'a', (0.4916, 0.4916, 0.4916), (0.0005,) * 3, 0.023),
123+
'O4': ('O2-', 4, 'a', (0.2336, 0.2336, 0.2336), (0.0003,) * 3, 0.037),
124+
}
125+
126+
assert structure.atom_sites.names == list(expected_sites)
127+
for label, expected in expected_sites.items():
128+
type_symbol, multiplicity, wyckoff, coordinates, coordinate_sus, adp_iso = expected
129+
site = structure.atom_sites[label]
130+
assert site.type_symbol.value == type_symbol
131+
assert site.multiplicity.value == multiplicity
132+
assert site.wyckoff_letter.value == wyckoff
133+
assert site.occupancy.value == 1.0
134+
np.testing.assert_allclose(
135+
[site.fract_x.value, site.fract_y.value, site.fract_z.value],
136+
coordinates,
137+
)
138+
np.testing.assert_allclose(
139+
[site.fract_x.uncertainty, site.fract_y.uncertainty, site.fract_z.uncertainty],
140+
coordinate_sus,
141+
)
142+
assert site.adp_iso.value == pytest.approx(adp_iso)
143+
assert site.adp_iso.uncertainty == pytest.approx(0.001)
144+
145+
cryspy_cif = CryspyCalculator()._convert_structure_to_cryspy_cif(structure)
146+
cryspy_structure = str_to_globaln(cryspy_cif).items[0]
147+
148+
assert cryspy_structure.data_name == '83267-icsd'
149+
assert [site.type_symbol for site in cryspy_structure.atom_site.items] == [
150+
expected[0] for expected in expected_sites.values()
151+
]

0 commit comments

Comments
 (0)