diff --git a/docs/CPU_DATA_ANALYSIS.md b/docs/CPU_DATA_ANALYSIS.md
index 31c3b969..129bce8b 100644
--- a/docs/CPU_DATA_ANALYSIS.md
+++ b/docs/CPU_DATA_ANALYSIS.md
@@ -604,6 +604,6 @@ A linear regression of $\log\langle I\rangle$ vs $1/d^2$ provides an estimate of
 - **Rotation vs still modes** differ substantially in prediction and scaling: partiality is angle-driven in rotation data, while stills are predicted (within an excitation-error window) and scaled with unit partiality.
 - **Space-group determination.** When no space group is supplied, a POINTLESS-like search scores Laue-group symmetry (CC of $I(h)$ vs $I(Rh)$ plus merge self-consistency) and detects screw/centering absences from the $P1$-merged intensities.
 - **Twinning check.** A Padilla–Yeates $L$-test ($\langle|L|\rangle$, $\langle L^2\rangle$) and the second moment $\langle I^2\rangle/\langle I\rangle^2$ are computed and written to the merged mmCIF as a twinning diagnostic.
-- **Outlier rejection.** Merging applies an optional per-observation median-based $N\sigma$ cut (default 6σ for `rot3d`) and an optional per-crystal $\Delta\mathrm{CC}_{1/2}$ image rejection (`--reject-delta-cchalf`, CrystFEL-style, off by default).
+- **Outlier rejection.** Merging applies an optional per-observation median-based $N\sigma$ cut (default 6σ for `rot3d`) and an optional per-crystal $\Delta\mathrm{CC}_{1/2}$ image rejection (`--reject-delta-cchalf`, CrystFEL-style, off by default). The same $N\sigma$ cut is fed back into the error model: after an initial $a,b$ fit the parameters are re-fit once on the reflections that survive rejection (dropping any whose squared deviation exceeds $N\sigma^2\,[a\,\sigma^2 + (b\,\langle I\rangle)^2]$), so the calibrated errors describe the reflections that actually enter the merge rather than the pre-rejection pool.
 - **Automatic resolution cutoff.** By default the reported/written high-resolution limit is trimmed where $\mathrm{CC}_{1/2}$ falls off (logistic, target 0.30); `--scaling-high-resolution` overrides it and `--resolution-cutoff off` disables it.
 - **Intensities only.** The merged output carries intensities (mmCIF `intensity_meas`, MTZ `IMEAN`/`SIGIMEAN`); it does not convert to amplitudes $|F|$ (no French–Wilson / truncate step) — do that downstream.
diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp
index 52aa9f66..6c6aa543 100644
--- a/image_analysis/scale_merge/RotationScaleMerge.cpp
+++ b/image_analysis/scale_merge/RotationScaleMerge.cpp
@@ -913,17 +913,20 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool
                 }
         }
         constexpr int n_bins = 16;
-        if (samples.size() >= static_cast<size_t>(8 * n_bins)) {
-            std::sort(samples.begin(), samples.end(),
-                      [](const Sample &a, const Sample &b) { return a.I2 < b.I2; });
+        // Fit (a, b) from the intensity-binned median deviations. Factored into a lambda so it can be
+        // re-run on a misfit-free pool below; takes the samples by value (it sorts them in place).
+        auto fit_ab = [&](std::vector<Sample> smp) {
+            if (smp.size() < static_cast<size_t>(8 * n_bins))
+                return;
+            std::sort(smp.begin(), smp.end(), [](const Sample &a, const Sample &b) { return a.I2 < b.I2; });
             std::vector<double> bs2, bI2, bd2;
-            const size_t per = samples.size() / n_bins;
+            const size_t per = smp.size() / n_bins;
             for (int bin = 0; bin < n_bins; ++bin) {
                 const size_t lo = bin * per;
-                const size_t hi = (bin == n_bins - 1) ? samples.size() : lo + per;
+                const size_t hi = (bin == n_bins - 1) ? smp.size() : lo + per;
                 std::vector<double> vs2, vI2, vd2;
                 for (size_t i = lo; i < hi; ++i) {
-                    vs2.push_back(samples[i].s2); vI2.push_back(samples[i].I2); vd2.push_back(samples[i].dev2);
+                    vs2.push_back(smp[i].s2); vI2.push_back(smp[i].I2); vd2.push_back(smp[i].dev2);
                 }
                 bs2.push_back(median_of(vs2));
                 bI2.push_back(median_of(vI2));
@@ -946,13 +949,29 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool
                 error_model_b = std::sqrt(b2);
                 error_model_active = true;
                 std::vector<double> chi2;
-                chi2.reserve(samples.size());
-                for (const auto &s : samples) {
+                chi2.reserve(smp.size());
+                for (const auto &s : smp) {
                     const double v = error_model_a * s.s2 + b2 * s.I2;
                     if (v > 0.0) chi2.push_back(s.dev2 / v);
                 }
                 error_model_chi2 = chi2.empty() ? 0.0 : median_of(chi2) / CHI2_1_MEDIAN;
             }
+        };
+        fit_ab(samples);
+        // Refit on a misfit-free pool: the merge drops symmetry outliers (|I - median| > reject_nsigma *
+        // sigma) from the merged intensity, so drop the equivalent samples (dev2 > reject_nsigma^2 * model
+        // variance) from the error-model fit too, keeping the fitted sigmas consistent with the reflections
+        // that actually survive. Operates on the shared `samples`, so CPU and GPU stay bit-identical.
+        if (reject_outliers && error_model_active) {
+            const double ns2 = reject_nsigma * reject_nsigma, b2 = error_model_b * error_model_b;
+            std::vector<Sample> kept;
+            kept.reserve(samples.size());
+            for (const auto &s : samples) {
+                const double v = error_model_a * s.s2 + b2 * s.I2;
+                if (v > 0.0 && s.dev2 <= ns2 * v) kept.push_back(s);
+            }
+            if (kept.size() >= static_cast<size_t>(8 * n_bins) && kept.size() < samples.size())
+                fit_ab(std::move(kept));
         }
     }
     if (error_model_active)