From 5ff7fbafe492410b6c28232506df39f5e73e3600 Mon Sep 17 00:00:00 2001 From: WEN Hao Date: Sat, 1 Aug 2026 21:57:22 +0800 Subject: [PATCH 1/2] fix: align wang method prob2step plateau logic with upstream R implementation Three regressions introduced in 8267a97 (revert) caused Python/R divergence: 1. _prob2step/_prob2steplmin: Replaced the isclose-based plateau detection with argmax/argmin. On symmetric profiles the old code bracketed both mirror peaks, causing the second-stage grid-refinement to miss the true maximum and systematically underestimate probmid by up to ~7e-3. This cascaded through the binary search to shift confidence bounds by as much as 0.02. argmax/argmin matches R's which(sumofprob == max(sumofprob)) semantics (single-index, ulp-stable breaking of exact ties). 2. argsort: Restored kind="stable" so that tied rows retain the same order as R's stable order(). 3. np.unique: Switched to a return_index + sort(uniq_idx) pattern that preserves first-occurrence order, matching R's unique() behavior instead of NumPy's lexicographic sort. The exact CI-failing case (8,37,23,37) now produces identical results: R (-0.59305, -0.15878) = Python (-0.59305, -0.15878). --- diff_binom_confint/_specials/_wang.py | 44 ++++++++------------------- diff_binom_confint/_utils.py | 2 +- test/test_utils.py | 8 +++-- 3 files changed, 19 insertions(+), 35 deletions(-) diff --git a/diff_binom_confint/_specials/_wang.py b/diff_binom_confint/_specials/_wang.py index ddd2eca..3b6cd30 100644 --- a/diff_binom_confint/_specials/_wang.py +++ b/diff_binom_confint/_specials/_wang.py @@ -206,7 +206,7 @@ def binomial_ci_one_sided( f[:, 2] = (p1hat - p0hat) / np.sqrt(denom) # Sort f by the third column in descending order - f = f[(-f[:, 2]).argsort(), :] + f = f[(-f[:, 2]).argsort(kind="stable"), :] allvector = np.round(f[:, 0] * (m + 2) + f[:, 1]).astype(int) allvectormove = np.round((f[:, 0] + 1) * (m + 3) + (f[:, 1] + 1)).astype(int) @@ -268,8 +268,10 @@ def binomial_ci_one_sided( dd[:, 1] += 1 b = dd - # Generate N - n_arr = np.unique(np.vstack((a, b)), axis=0) + # Generate N (order-preserving unique, matching R's unique()) + n_rows = np.vstack((a, b)) + _, uniq_idx = np.unique(n_rows, axis=0, return_index=True) + n_arr = n_rows[np.sort(uniq_idx)] nvector = ((n_arr[:, 0] + 1) * (m + 3) + n_arr[:, 1] + 1).astype(int) nvector = nvector[np.isin(nvector, allvectormove)] @@ -365,7 +367,7 @@ def binomial_ci_one_sided( if length_nc >= 2: valid = ~np.isnan(nc_arr[:, 0]) ncnomiss = nc_arr[valid] - ncnomiss = ncnomiss[(-ncnomiss[:, 2]).argsort(), :] + ncnomiss = ncnomiss[(-ncnomiss[:, 2]).argsort(kind="stable"), :] morepoint = np.sum(ncnomiss[:, 2] >= ncnomiss[0, 2] - delta) if morepoint >= 2: ls_arr[kk : kk + morepoint, 0:2] = ncnomiss[:morepoint, 0:2] @@ -464,21 +466,10 @@ def _prob2step(delv, delta, n, m, i1, i2, grid_one, grid_two): part2 = np.log(comb(m, i2))[:, None] + np.outer(i2, np.log(p0)) + np.outer(m - i2, np.log(1 - p0)) sumofprob = np.exp(part1 + part2).sum(axis=0) - # plateau-aware refinement (R: which(sumofprob == max(sumofprob))) - mansum = sumofprob.max() - atol = 1e-14 * (mansum if mansum > 0 else 1.0) - plateau_idx = np.where(np.isclose(sumofprob, mansum, rtol=0.0, atol=atol))[0] - leftmost = plateau_idx.min() - rightmost = plateau_idx.max() - stepv = (p0[-1] - p0[0]) / grid_one - lowerb = max(p0[0], p0[rightmost] - stepv) + delta - upperb = min(p0[-1], p0[leftmost] + stepv) - delta - - # stepv = (p0[-1] - p0[0]) / grid_one - # maxloc = np.argmax(sumofprob) - # lowerb = max(p0[0], p0[maxloc] - stepv) + delta - # upperb = min(p0[-1], p0[maxloc] + stepv) - delta + maxloc = np.argmax(sumofprob) + lowerb = max(p0[0], p0[maxloc] - stepv) + delta + upperb = min(p0[-1], p0[maxloc] + stepv) - delta p0 = np.linspace(lowerb, upperb, grid_two) part1 = np.log(comb(n, i1))[:, None] + np.outer(i1, np.log(p0 + delv)) + np.outer(n - i1, np.log(1 - p0 - delv)) @@ -498,21 +489,10 @@ def _prob2steplmin(delv, delta, n, m, i1, i2, grid_one, grid_two): part2 = np.log(comb(m, i2))[:, None] + np.outer(i2, np.log(p0)) + np.outer(m - i2, np.log(1 - p0)) sumofprob = np.exp(part1 + part2).sum(axis=0) - # plateau-aware refinement for minima (R: which(sumofprob == min(sumofprob))) - mansum = sumofprob.min() - atol = 1e-14 * (abs(mansum) if mansum != 0 else 1.0) - plateau_idx = np.where(np.isclose(sumofprob, mansum, rtol=0.0, atol=atol))[0] - leftmost = plateau_idx.min() - rightmost = plateau_idx.max() - stepv = (p0[-1] - p0[0]) / grid_one - lowerb = max(p0[0], p0[rightmost] - stepv) + delta - upperb = min(p0[-1], p0[leftmost] + stepv) - delta - - # stepv = (p0[-1] - p0[0]) / grid_one - # minloc = np.argmin(sumofprob) - # lowerb = max(p0[0], p0[minloc] - stepv) + delta - # upperb = min(p0[-1], p0[minloc] + stepv) - delta + minloc = np.argmin(sumofprob) + lowerb = max(p0[0], p0[minloc] - stepv) + delta + upperb = min(p0[-1], p0[minloc] + stepv) - delta p0 = np.linspace(lowerb, upperb, grid_two) part1 = np.log(comb(n, i1))[:, None] + np.outer(i1, np.log(p0 + delv)) + np.outer(n - i1, np.log(1 - p0 - delv)) diff --git a/diff_binom_confint/_utils.py b/diff_binom_confint/_utils.py index e55dbfc..cd194d9 100644 --- a/diff_binom_confint/_utils.py +++ b/diff_binom_confint/_utils.py @@ -135,7 +135,7 @@ def remove_parameters_returns_from_docstring( indices2remove.extend(list(range(start_idx, idx))) start_idx = None if start_idx is not None: - indices2remove(list(range(start_idx, len(new_doc)))) + indices2remove.extend(list(range(start_idx, len(new_doc)))) new_doc.extend(["\n", parameters_indicator or returns_indicator]) new_doc = [line for idx, line in enumerate(new_doc) if idx not in indices2remove] # remove trailing empty lines diff --git a/test/test_utils.py b/test/test_utils.py index d050426..a39c445 100644 --- a/test/test_utils.py +++ b/test/test_utils.py @@ -1,5 +1,6 @@ """ """ +import inspect import random import time @@ -58,7 +59,10 @@ def test_remove_parameters_returns_from_docstring(): parameters=["returns_indicator", "parameters_indicator"], returns="new_doc", ) - assert new_docstring == """Remove parameters and/or returns from docstring, + # Since Python 3.13, docstrings are dedented at compile time, so + # `remove_parameters_returns_from_docstring.__doc__` has no common leading + # whitespace; compare normalized forms to be version-agnostic. + assert inspect.cleandoc(new_docstring) == inspect.cleandoc("""Remove parameters and/or returns from docstring, which is of the format of numpydoc. Parameters @@ -74,7 +78,7 @@ def test_remove_parameters_returns_from_docstring(): ------- None - """ + """) new_docstring = remove_parameters_returns_from_docstring( remove_parameters_returns_from_docstring.__doc__, From 11bfbff7d70545e0d20e43b96bd1fdcdbe0d1826 Mon Sep 17 00:00:00 2001 From: WEN Hao Date: Sat, 1 Aug 2026 22:23:43 +0800 Subject: [PATCH 2/2] fix: do not append synthetic section header when removing terminal docstring item When a parameter or return value at EOF was the last item in its section, the old code appended a stray "Parameters" / "Returns" header line. Just collect the remaining lines into indices2remove and let the trailing-empty-line cleanup handle the rest. --- diff_binom_confint/_utils.py | 1 - 1 file changed, 1 deletion(-) diff --git a/diff_binom_confint/_utils.py b/diff_binom_confint/_utils.py index cd194d9..21a068e 100644 --- a/diff_binom_confint/_utils.py +++ b/diff_binom_confint/_utils.py @@ -136,7 +136,6 @@ def remove_parameters_returns_from_docstring( start_idx = None if start_idx is not None: indices2remove.extend(list(range(start_idx, len(new_doc)))) - new_doc.extend(["\n", parameters_indicator or returns_indicator]) new_doc = [line for idx, line in enumerate(new_doc) if idx not in indices2remove] # remove trailing empty lines idx = max(idx for idx, line in enumerate(new_doc) if len(line.strip()) > 0)