diff --git a/openavmkit/utilities/stats.py b/openavmkit/utilities/stats.py index dbebed80..11ed866b 100644 --- a/openavmkit/utilities/stats.py +++ b/openavmkit/utilities/stats.py @@ -583,12 +583,12 @@ def calc_prb( model = sm.OLS(left, right).fit() # Guard against degenerate fit (rare but better to be explicit) - if model.df_resid <= 0 or not np.isfinite(model.params[0]): + if model.df_resid <= 0 or not np.isfinite(model.params[1]): return np.nan, np.nan, np.nan - prb = float(model.params[0]) + prb = float(model.params[1]) prb_lower, prb_upper = ( - model.conf_int(alpha=1.0 - confidence_interval)[0].tolist() + model.conf_int(alpha=1.0 - confidence_interval)[1].tolist() ) return prb, prb_lower, prb_upper diff --git a/tests/test_stats.py b/tests/test_stats.py index 0373d64d..e7dca590 100644 --- a/tests/test_stats.py +++ b/tests/test_stats.py @@ -37,4 +37,31 @@ def test_cod_bootstrap(): print(expected) print("***") - assert objects_are_equal(results, expected) \ No newline at end of file + assert objects_are_equal(results, expected) + +def test_calc_prb_returns_slope(): + import numpy as np + import statsmodels.api as sm + from openavmkit.utilities.stats import calc_prb + + n = 500 + truth = np.linspace(100_000, 1_000_000, n) + rank = np.linspace(0.0, 1.0, n) + preds = truth * (0.9 + 0.2 * rank) # ratios rise with parcel value + + prb, lo, hi = calc_prb(preds, truth) + + # fit the identically transformed regression by hand + ratios = preds / truth + med = np.median(ratios) + left = (ratios - med) / med + right = sm.add_constant(np.log2(preds / med + truth), has_constant="add") + slope = sm.OLS(left, right).fit().params[1] + + assert abs(prb - slope) < 1e-9 + assert prb > 0 + assert lo <= prb <= hi + + # unbiased control: a flat multiplier of ground truth gives PRB near zero + prb0, _, _ = calc_prb(truth * 0.95, truth) + assert abs(prb0) < 1e-6