Research: the hardest problems, with proofs

morie.research carries the research programme on the hardest open problems in criminology and sociolegal studies. Every function rests on a theorem checked in Lean 4 with Mathlib (research/lean in the rmorie repository: 0 sorry, axioms propext, Classical.choice and Quot.sound only), names the theorems it uses in its theorems field, and says where the empirical assumption enters. Lean certifies the implication, never the antecedent: no proof shows that two lists are independent, that a noise box is right, or that spillover stops at the ring the analyst named.

The R package (rmorie) carries the same functions with the morie_ prefix. The two arms agree to rounding on the same inputs, and the stochastic ones draw from the shared Philox stream, so a seed gives the same path in both.

from morie import research as R

R.dark_figure_two_source(400, 250, 80, kappa=2)   # Lincoln-Petersen with a dependence box
R.feedback_loop_limit(0.3, 0.2, 10, 10)            # the proved limit of the patrol share
R.meta_random_effects([0.20, 0.35, 0.10], [0.010, 0.020, 0.015])  # DerSimonian-Laird, with the truncation flag

Problems and functions

Problem	Theorems (Research.*)	Functions
P1 dark figure	P1.TwoSource.petersen_bounds, P1.true_rate_bounds, P1.dark_figure_bounds, P1.conclusion_holds_below_breakdown, P1.three_list_saturated_fits, P1.offence_count_bounds	dark_figure_two_source(), dark_figure_bounds(), dark_figure_breakdown(), dark_figure_three_list(), dark_figure_hierarchy()
P1 Le Cam two-point bound	P1LeCam.sum_min, P1LeCam.two_point, P1LeCam.minimax	two_point_bound()
P2 selection in police records	P2.rate_bounds, P2.disparity_bounds, P2.benchmark_product, P2.offset_shift, P2.rr_between, P2.dyad_ratio, P2.collider_or_eq_background	disparity_exposure_bounds(), disparity_benchmark(), relative_risk_from_or(), interracial_rates(), collider_arrest()
P3 interference	P3.Model.exposure_adjustment, P3.Design.ht_unbiased, P3.Design.ht_variance, P3.SYG.syg_eq_ht, P3.cheeger_easy	spillover_exposure(), spillover_effects(), spillover_ht(), spillover_ht_variance(), spillover_exposure_probs(), cheeger_bound()
P4 predictive-policing feedback	P4.naive_step_drift, P4.naiveShare_tendsto_one, P4.corrected_share_tendsto, P4.naiveShare_rate_bound, P4.rho_cap, P4.Urn.polya_uniform	feedback_loop_meanfield(), feedback_loop_limit(), feedback_loop_urn_law(), feedback_loop_sim(), feedback_loop_bound()
P5 risk scores under label bias	P5.Table.chouldechova, P5.impossibility, P5.true_base_rate_bounds, P5.compare_decided, P5.reduced_coefficient, P5.rank_stable_of_gap, P5.hr2_gt_one_of_depletion, P5.no_mle	fairness_rates(), fairness_implied_fpr(), fairness_base_rate_bounds(), fairness_true_rate(), fairness_compare_groups(), logit_rescale(), ranking_resolution(), hazard_selection(), logit_separation()
P6 age-crime curve	P6.aggregate_not_identifying	age_crime_aggregate()
P7 crime concentration	P7.gini_zero_decomposition, P7.poisson_zero_prob, P7.Mixture.mixture_var_ge_mean, P7.expectedDistinct_bounds	concentration_gini(), concentration_decompose(), concentration_dispersion(), concentration_distinct_growth()
P8 deterrence	P8.constant_dimension_not_identified, P8.certainty_monotone, P8.necessity_bounds	deterrence_design_check(), deterrence_response(), probability_of_necessity()
P9 recording as a linear map	P9.total_invariant_of_colStochastic, P9.detection_rate_rises	recording_map(), detection_rate_shift()
P10 near-repeat contagion	P10.cluster_size_of_lt_one, P10.endogeneity_share, P10.extinction_le_fixed, P10.supercritical_extinction_lt_one	contagion_branching(), contagion_extinction()
P11 sentencing effects as intervals	P11.Pop.outcome_bounds, P11.clean_bounds, P11.Pop.mtr_upper, P11.Pop.mts_ate_le_naive, P11.im_cutoff_antitone, P11.two_sided_overcovers	sentence_effect_bounds(), contaminated_bounds(), sentence_effect_mtr(), sentence_effect_mts(), bounds_confidence()
P12 ecological inference	P12.cov_decomp, P12.ecological_ge, P12.dd_bounds	ecological_decompose(), ecological_bounds()
P13 pooling evaluations	P13.truncation_bias, P13.dl_biased_under_homogeneity, P13.re_var_ge	meta_random_effects(), meta_dl_bias()
P13 HKSJ interval	P13HKSJ.hksj_wider_iff, P13HKSJ.Q_eq_zero_iff, P13HKSJ.hksj_equal_weights	meta_hksj()
P14 judge leniency	P14.itt_decomposition, P14.first_stage_decomposition, P14.late_identification, P14.wald_with_defiers, P14.defiers_can_flip	judge_iv_population(), judge_iv()
P14 many-judge slope test	P14Slope.propensity_mono, P14Slope.outcome_diff, P14Slope.slope_bound, P14Slope.violation_refutes_monotonicity	judge_slope_test()
P15 disparity decomposition	P15.twofold_B, P15.threefold, P15.reference_dependence, P15.attribution_shift	disparity_decomposition()
P15 DFL reweighting	P15Reweight.reweighting_matches, P15Reweight.reweighted_mass, P15Reweight.counterfactual_outcome	dfl_reweight()
P16 court backlog	P16.occupancy_integral, P16.little, P16.little_backlog, P16.little_target	court_backlog()
P16 disposed-cases mean as a bound	P16Censoring.true_mean_ge, P16Censoring.bias_lower, P16Censoring.disposed_understates, P16Censoring.no_upper_bound	backlog_censoring()
P17 incapacitation	P17.steady_state_rate, P17.prevented_share_lt_one, P17.marginal_prevention_eq, P17.high_rate_more_prevented	incapacitation()
P17 desistance and replacement	P17Replacement.prevented_le_const, P17Replacement.prevented_ge_const, P17Replacement.later_sentence_prevents_less, P17Replacement.prevented_net_le	incapacitation_career()
P18 selective labels	P18.nested_rate_identified, P18.unobserved_bounds, P18.unobserved_width	selective_labels()
P19 regression to the mean	P19.exchange_cross, P19.indicator_bound, P19.selected_change_nonpos	regression_to_mean()
P19 empirical-Bayes shrinkage	P19Shrinkage.loss_eq, P19Shrinkage.loss_min, P19Shrinkage.loss_bstar_le_raw, P19Shrinkage.predicted_fall	hotspot_shrinkage(), shrinkage_loss()

Reference

The research programme: methods for the hardest problems in criminology, each resting on a Lean 4 theorem.

Every function names the theorems in research/lean (Lean 4 + Mathlib, 0 sorry, standard axioms only) that cover its arithmetic, and says where the empirical assumption enters. Lean certifies the implication, never the antecedent. The R package carries the same functions with the morie_ prefix; both arms agree to rounding on the same inputs, and the stochastic ones share the Philox stream.

morie.research.age_crime_aggregate(shares, curves, ages=None)

Aggregate age-crime curve of a mixture of latent types (Research.P6.aggregate_not_identifying).

Examples

>>> a = age_crime_aggregate([0.5, 0.5], [[1, 3], [2, 2], [3, 1]], ages=[10, 11, 12])
>>> list(a["aggregate"])
[2.0, 2.0, 2.0]

morie.research.backlog_censoring(disposed, pending_ages)

The disposed-cases mean as a bound: what the pending cases imply.

Research.P16Censoring.true_mean_ge / lower_bound_sub / bias_lower / disposed_understates / no_upper_bound.

Examples

>>> b = backlog_censoring([30, 45, 60, 90, 120], [100, 150, 200])
>>> (b["disposed_mean"], b["lower_bound"], b["bias_lower"], b["understates"], b["upper_bound"])
(69.0, 99.375, 30.375, True, inf)

Return type:

dict

morie.research.bounds_confidence(lower, upper, se_lower, se_upper, level=0.95)

Imbens-Manski confidence interval for a partially identified parameter.

The cutoff c solves Phi(c + Delta) - Phi(-c) = 1 - alpha with Delta = (upper - lower)/max(se) it is non-increasing in Delta (Research.P11.im_cutoff_antitone), lies between the one- and two-sided quantiles (im_cutoff_between), and the two-sided quantile over-covers any region of positive width (two_sided_overcovers).

Examples

>>> r = bounds_confidence(lower=0.10, upper=0.35, se_lower=0.03, se_upper=0.04)
>>> round(r["cutoff"], 6)
1.644854

Return type:

dict

morie.research.cheeger_bound(adjacency, S, exhaustive=False)

Conductance of a set and the Cheeger bound lambda_2 <= E(f)/D(f) <= 2 h(S).

Examples

>>> A = [[0, 1, 1, 0, 0, 0], [1, 0, 1, 0, 0, 0], [1, 1, 0, 1, 0, 0],
...      [0, 0, 1, 0, 1, 1], [0, 0, 0, 1, 0, 1], [0, 0, 0, 1, 1, 0]]
>>> r = cheeger_bound(A, S=[0, 1, 2], exhaustive=True)
>>> (round(r["conductance"], 12), round(r["rayleigh_test"], 12), r["bound_holds"], r["argmin_set"])
(0.142857142857, 0.285714285714, True, [0, 1, 2])

Return type:

dict

morie.research.collider_arrest(a, e, pi_bg)

Collider bias from conditioning on arrest (collider_or_eq_background, collider_or_lt_one).

Examples

>>> r = collider_arrest(a=0.3, e=0.2, pi_bg=0.1)
>>> (round(r["arrestee_or"], 12), r["population_or"])
(0.1, 1)

Return type:

dict

morie.research.concentration_decompose(x)

Decompose crime concentration into crime-free places and the rest (gini_zero_decomposition).

Examples

>>> r = concentration_decompose([0, 0, 0, 1, 9])
>>> (r["zero_share"], round(r["gini_positive"], 12), round(r["identity_check"], 15))
(0.6, 0.4, 0.0)

Return type:

dict

morie.research.concentration_dispersion(x)

Dispersion of place counts against the Poisson null and its mixtures.

Examples

>>> r = concentration_dispersion([0, 0, 1, 3, 0, 2])
>>> (r["mean_count"], round(r["dispersion_index"], 12))
(1.0, 1.6)

Return type:

dict

morie.research.concentration_distinct_growth(place)

Growth of the number of distinct places under Polya allocation (expectedDistinct_bounds).

Examples

>>> r = concentration_distinct_growth([1, 2, 1, 3, 2, 1, 4, 1])
>>> (r["n"], r["distinct"][-1], round(r["M_hat"], 6))
(8, 4, 2.500624)

Return type:

dict

morie.research.concentration_gini(x)

Gini coefficient of a count vector (mean-absolute-difference form).

Examples

>>> round(concentration_gini([0, 0, 0, 1, 9]), 12)
0.76

Return type:

float

morie.research.contagion_branching(n, mu=1.0, generations=10)

Branching-ratio arithmetic of a self-exciting (Hawkes) crime process.

Examples

>>> r = contagion_branching(n=0.4, mu=2)
>>> (round(r["expected_cluster_size"], 12), round(r["stationary_rate"], 12), r["endogeneity_share"])
(1.666666666667, 3.333333333333, 0.4)
>>> contagion_branching(n=1.1, mu=2)["expected_cluster_size"]
inf

Return type:

dict

morie.research.contagion_extinction(p, tol=1e-14, max_iter=100000)

Extinction probability of a near-repeat chain: the smallest fixed point of the generating function.

Iterates s_0 = 0, s_{n+1} = f(s_n) (Research.P10.iter_tendsto, extinction_fixed, extinction_le_fixed); mean offspring below one gives certain extinction (subcritical_extinction_one), above one a survival probability strictly positive (supercritical_extinction_lt_one).

Examples

>>> r = contagion_extinction([0.3, 0.3, 0.4])
>>> (round(r["extinction"], 9), r["regime"])
(0.75, 'supercritical')
>>> round(contagion_extinction([0.5, 0.3, 0.2])["extinction"], 9)
1.0

Return type:

dict

morie.research.contaminated_bounds(q, p)

Contaminated-sample bounds for a recorded proportion (Research.P11.clean_bounds, clean_width, clean_informative).

Examples

>>> b = contaminated_bounds([0.05, 0.5, 0.97], 0.1)
>>> [round(v, 12) for v in b["lower"]]
[0.0, 0.444444444444, 0.966666666667]

morie.research.court_backlog(arrivals, dispositions, horizon=None, target_backlog=None)

Little’s law on a docket: time-average pending cases = filing rate x mean disposition time.

Examples

>>> r = court_backlog([0, 1, 2, 4, 5, 7], [3, 2.5, 6, 5, 9, 10], horizon=10, target_backlog=1)
>>> (round(r["average_backlog"], 12), round(r["required_mean_time"], 12))
(1.65, 1.666666666667)

Return type:

dict

morie.research.dark_figure_bounds(v_obs, r, alpha_max, beta_max)

Sharp bounds on a true victimisation rate and its dark figure.

max(r, (v_obs - a)/(1 - a)) <= v <= v_obs/(1 - b) with the boxes alpha in [0, a], beta in [0, b], a + b < 1 (Research.P1.true_rate_bounds, dark_figure_bounds). Returns a frame with v_obs, r, v_lower, v_upper, dark_lower, dark_upper, ratio_upper.

Examples

>>> b = dark_figure_bounds(0.06, 0.02, alpha_max=0.01, beta_max=0.30)
>>> round(float(b["v_upper"][0]), 10)
0.0857142857

morie.research.dark_figure_breakdown(v_obs, threshold)

Breakdown analysis for a dark-figure conclusion.

“The true rate is below threshold” survives every admissible noise pair while the under-reporting box stays below 1 - v_obs/threshold (Research.P1.conclusion_holds_below_breakdown) and fails at or above it (conclusion_fails_above_breakdown).

Examples

>>> round(dark_figure_breakdown(0.06, threshold=0.10)["breakdown_beta_max"], 12)
0.4

Return type:

dict

morie.research.dark_figure_hierarchy(offences_per_incident)

Incident versus offence counting: the hierarchy-rule arithmetic.

With N incidents carrying k_i >= 1 offences, bounded by K, the offence count lies in [N, KN] and equals N(1 + mean extra) (Research.P1.offence_count_bounds, offence_count_eq); it is not a function of N (category_not_identified).

Examples

>>> r = dark_figure_hierarchy([1, 1, 2, 1, 3, 1, 1, 2])
>>> (r["incidents"], r["offences"], r["mean_extra"])
(8, 12, 0.5)

Return type:

dict

morie.research.dark_figure_three_list(counts, candidate_missing=None)

Three-list capture-recapture: what is and is not identified.

The saturated model reproduces any positive eight-cell table (Research.P1.three_list_saturated_fits), so the missing cell is free (missing_cell_unconstrained); setting the three-way interaction to zero gives m000 = m111 m100 m010 m001 / (m110 m101 m011). Pairwise Petersen and Chapman estimates respect the floor (petersen_ge_floor, chapman_ge_floor).

Examples

>>> r = dark_figure_three_list({"100": 120, "010": 90, "001": 70, "110": 40, "101": 30, "011": 25, "111": 15})
>>> (r["observed"], round(r["missing_no_three_way"], 6))
(390.0, 378.0)

Return type:

dict

morie.research.dark_figure_two_source(n1, n2, m, kappa=1.0, chapman=False)

Two-source (capture-recapture) count with a dependence box.

Under independence the true count is the Lincoln-Petersen ratio n1 n2 / m. With the dependence factor only known to lie in [1/kappa, kappa] the true count lies in [n1 n2 / (kappa m), kappa n1 n2 / m], both ends attainable (Research.P1.TwoSource.petersen_bounds).

Examples

>>> dark_figure_two_source(400, 250, 80)["point"]
1250.0
>>> r = dark_figure_two_source(400, 250, 80, kappa=2)
>>> (r["lower"], r["upper"])
(625.0, 2500.0)

Return type:

dict

morie.research.detection_rate_shift(detected, total, n, d, q)

Detection-rate arithmetic of downgrade-and-caution (Research.P9.detection_rate_rises).

Examples

>>> r = detection_rate_shift(detected=2000, total=10000, n=1500, d=0.2, q=0.5)
>>> (r["rate_before"], r["rate_after"], r["rise"])
(0.2, 0.26, 0.06)

Return type:

dict

morie.research.deterrence_design_check(p, s, c)

Which deterrence partial effects a design can identify (Research.P8.constant_dimension_not_identified).

Examples

>>> r = deterrence_design_check(p=[0.1, 0.3, 0.5], s=[2, 2, 2], c=[30, 30, 30])
>>> (r["rank"], r["identified"])
(2, {'certainty': True, 'severity': False, 'celerity': False})

Return type:

dict

morie.research.deterrence_response(x, benefit, sanction, p)

Direction of deterrence without convexity (certainty_monotone, severity_monotone, aggregate_monotone).

Examples

>>> r = deterrence_response([0, 1, 2, 3], [0, 2, 3, 3.5], [0, 1, 3, 6], p=[0.1, 0.5, 1])
>>> list(r["x_opt"])
[3.0, 2.0, 1.0]

morie.research.dfl_reweight(group, x, y, weights=None)

DiNardo-Fortin-Lemieux reweighting: composition and structure without a linear model.

Research.P15Reweight.reweighting_matches / reweighted_mass / counterfactual_outcome / decomposition.

Examples

>>> g = [True] * 6 + [False] * 6
>>> x = ["a", "a", "a", "a", "b", "b",  "a", "a", "b", "b", "b", "b"]
>>> y = [10, 12, 11, 13, 20, 22,  8, 9, 15, 16, 14, 17]
>>> r = dfl_reweight(g, x, y)
>>> tuple(round(r[k], 12) for k in ("mean_1", "mean_0", "counterfactual", "structure", "composition"))
(14.666666666667, 13.166666666667, 10.833333333333, 3.833333333333, -2.333333333333)
>>> (r["psi"]["x"], r["psi"]["psi"], r["max_composition_gap"])
(['a', 'b'], [2.0, 0.5], 0.0)

Return type:

dict

morie.research.disparity_benchmark(pop, contact, force, reference, exposure_error_factor=None)

Disparity benchmarks: the product identity and the exposure-offset shift (Research.P2.benchmark_product, benchmark_not_additive, offset_shift, disparity_ratio_shift).

Examples

>>> d = disparity_benchmark({"A": 100, "B": 100}, {"A": 30, "B": 10}, {"A": 12, "B": 2}, reference="B")
>>> ([round(v, 12) for v in d["resident_disparity"]], [round(v, 12) for v in d["product_check"]])
([6.0, 1.0], [0.0, 0.0])

morie.research.disparity_decomposition(y, X, group, reference, names=None, shift=1.0)

Oaxaca-Blinder decomposition of a gap in means, with both references and the interaction.

X holds the covariates (one column each, no constant: it is added), group the two-level label per row, reference the level whose coefficients price the explained part (group B).

Examples

>>> y = [10, 12, 13, 15, 9, 8, 11, 7]; X = [[1, 0], [2, 1], [3, 1], [4, 0], [1, 1], [2, 0], [2, 0], [3, 1]]
>>> r = disparity_decomposition(y, X, ["A"] * 4 + ["B"] * 4, reference="B", names=["x1", "x2"])
>>> (round(r["gap"], 12), [round(abs(v), 12) for v in r["identity_checks"].values()])
(3.75, [0.0, 0.0, 0.0, 0.0])
>>> ({k: round(v, 12) for k, v in r["twofold_B"].items()}, round(r["threefold"]["interaction"], 12))
({'explained': -0.5, 'unexplained': 4.25}, 1.3)

Return type:

dict

morie.research.disparity_exposure_bounds(y, m, gamma=1.0, reference=None)

Bounds on a police-outcome disparity when exposure is measured by a proxy (Research.P2.rate_bounds, disparity_bounds, disparity_sign_identified).

y and m are dicts keyed by group.

Examples

>>> d = disparity_exposure_bounds({"A": 300, "B": 100}, {"A": 1000, "B": 1000}, gamma=1.5, reference="B")
>>> ([round(v, 12) for v in d["ratio_proxy"]], [round(v, 12) for v in d["ratio_lower"]], list(d["direction_identified"]))
([3.0, 1.0], [1.333333333333, 0.444444444444], [True, False])

morie.research.ecological_bounds(p, q, weights=None)

Duncan-Davis bounds: what neighbourhood marginals say about an individual rate.

max(0, (p+q-1)/p) <= P(y | x) <= min(1, q/p), both ends attained (Research.P12.dd_bounds, Cells.ends_attained); the complement rate follows from q = p r + (1-p) r' (dd_complement); the aggregate rate inherits the m_g p_g-weighted mean of the intervals (dd_aggregate_bounds).

Examples

>>> r = ecological_bounds([0.2, 0.5, 0.8], [0.1, 0.3, 0.6], weights=[1000, 2000, 500])
>>> ([round(float(v), 12) for v in r["neighbourhoods"]["upper"]], round(r["aggregate"]["upper"], 12))
([0.5, 0.6, 0.75], 0.625)

Return type:

dict

morie.research.ecological_decompose(x, y, group)

Between/within decomposition of a correlation across groups.

Population (divide by n) conventions throughout, matching the Lean definitions.

Examples

>>> r = ecological_decompose([0, 2, 1, 3], [1, 3, 0, 2], ["a", "a", "b", "b"])
>>> (round(r["corr_individual"], 12), r["corr_ecological"], r["sign_reversed"])
(0.6, -1.0, True)

Return type:

dict

morie.research.fairness_base_rate_bounds(p_obs, alpha_max, beta_max)

Sharp bounds on a true base rate from a noisy recorded one (Research.P5.true_base_rate_bounds).

Examples

>>> b = fairness_base_rate_bounds([0.35, 0.55], alpha_max=0.10, beta_max=0.20)
>>> [round(v, 12) for v in b["upper"]]
[0.4375, 0.6875]

morie.research.fairness_compare_groups(p_obs_a, p_obs_b, alpha_max, beta_max)

Can two groups’ true base rates be ordered under label noise? (compare_decided, compare_undecided)

Examples

>>> fairness_compare_groups(0.35, 0.55, alpha_max=0.05, beta_max=0.20)["order"]
'a < b'
>>> fairness_compare_groups(0.35, 0.55, alpha_max=0.05, beta_max=0.40)["order"]
'undecided'

Return type:

dict

morie.research.fairness_implied_fpr(p, ppv, fnr)

Chouldechova’s identity fpr = p/(1-p) (1-ppv)/ppv (1-fnr) (Research.P5.Table.chouldechova).

Examples

>>> r = fairness_rates(120, 60, 40, 280)
>>> round(fairness_implied_fpr(r["p"], r["ppv"], r["fnr"]) - r["fpr"], 15)
0.0

Return type:

float

morie.research.fairness_rates(tp, fp, fn, tn)

Group-wise rates from a confusion table.

Examples

>>> r = fairness_rates(tp=120, fp=60, fn=40, tn=280)
>>> (round(r["p"], 12), round(r["ppv"], 12), round(r["fpr"], 12), r["fnr"])
(0.32, 0.666666666667, 0.176470588235, 0.25)

Return type:

dict

morie.research.fairness_true_rate(p_obs, alpha, beta)

Recover a true base rate when the noise rates are known (Research.P5.trueRate_observedRate).

Examples

>>> round(float(fairness_true_rate(0.4 * 0.8 + 0.6 * 0.1, alpha=0.1, beta=0.2)[0]), 12)
0.4

morie.research.feedback_loop_bound(lam_a, lam_b, c_a0, c_b0, n_steps=1000)

Proved upper bound on how far the naive loop can move in N steps (naiveShare_rate_bound).

Examples

>>> round(feedback_loop_bound(0.3, 0.2, 10, 10, n_steps=4)["bound"], 12)
0.004926706191

Return type:

dict

morie.research.feedback_loop_limit(lam_a, lam_b, c_a0, c_b0, update='naive', rho=0.0)

Proved limit of the two-region feedback loop, with the theorem that proves it.

Examples

>>> feedback_loop_limit(0.3, 0.2, 10, 10)["share_a"]
1
>>> feedback_loop_limit(0.3, 0.2, 10, 10, "corrected")["share_a"]
0.6
>>> round(feedback_loop_limit(0.3, 0.2, 10, 10, rho=0.5)["cap"], 12)
0.75

Return type:

dict

morie.research.feedback_loop_meanfield(lam_a, lam_b, c_a0, c_b0, n_steps=100, update='naive', rho=0.0)

Mean-field feedback-loop recursion for two regions.

Returns a frame with step, share_a, c_a, c_b; frame.attrs["limit"] holds the proved limit (see feedback_loop_limit()).

Examples

>>> mf = feedback_loop_meanfield(0.3, 0.2, 7, 3, n_steps=1)
>>> x = 0.7; drift = x * (1 - x) * 0.1 / (10 + 0.3 * x + 0.2 * (1 - x))
>>> round(float(mf["share_a"][1] - mf["share_a"][0]) - drift, 15)
0.0

morie.research.feedback_loop_sim(lam_a, lam_b, c_a0, c_b0, n_steps=1000, update='naive', rho=0.0, n_sims=100, seed=0)

Stochastic urn simulation of the two-region feedback loop on the shared Philox stream.

Examples

>>> s = feedback_loop_sim(0.3, 0.2, 10, 10, n_steps=50, n_sims=3, seed=1)
>>> (s["share_a"].shape, s["limit"]["theorem"])
((3, 51), 'Research.P4.naiveShare_tendsto_one')

Return type:

dict

morie.research.feedback_loop_urn_law(n_steps)

Exact law of the stochastic two-region urn with equal rates (polya_uniform).

Examples

>>> law = feedback_loop_urn_law(10)
>>> (round(float(law["prob"][0]), 12), law.attrs["prob_middle"] >= 0.25)
(0.090909090909, True)

morie.research.hazard_selection(s, h, l, survive_high, survive_low)

Built-in selection in period-by-period hazard ratios (survivor_hazard_mono, hr2_gt_one_of_depletion).

Examples

>>> r = hazard_selection(s=0.5, h=0.5, l=0.1, survive_high=(0.5, 0.8), survive_low=(0.9, 0.9))
>>> round(r["period2_hazard_ratio"], 12)
1.186851211073

Return type:

dict

morie.research.hotspot_shrinkage(y, noise_variance, weights=None, B=None)

Empirical-Bayes shrinkage of hot-spot counts: the expected size of the fall.

Research.P19Shrinkage.loss_min / bstar_mem / predicted_fall.

Examples

>>> s = hotspot_shrinkage([40, 12, 9, 25, 7, 31, 5, 18], noise_variance=18.375)
>>> (round(s["B"], 12), round(s["predicted_fall"][0], 12), round(s["shrunk"][0], 12), s["total_variance"])
(0.132686449284, 2.869344465757, 37.130655534243, 138.484375)

Return type:

dict

morie.research.incapacitation(lam, q, S, shares=None)

Crime rate, prevented share and the marginal sentence year under the incapacitation model.

Examples

>>> r = incapacitation(lam=[2, 10], q=0.1, S=1, shares=[0.8, 0.2])
>>> ([round(float(v), 12) for v in r["rate"]], {k: round(v, 12) for k, v in r.attrs["aggregate"].items()})
([1.666666666667, 5.0], {'free_rate': 3.6, 'incapacitated_rate': 2.333333333333, 'prevented_share': 0.351851851852})

morie.research.incapacitation_career(lambda_path, t0, S, replacement=0.0)

Incapacitation under desistance (a non-increasing rate path) and replacement.

Research.P17Replacement.prevented_le_const / prevented_ge_const / later_sentence_prevents_less / replaced_antitone / prevented_net_le.

Examples

>>> r = incapacitation_career([12, 10, 8, 6, 5, 4, 3, 2, 2, 1], t0=2, S=3, replacement=0.25)
>>> (r["prevented"], r["prevented_net"], r["upper"], r["lower"], r["later"])
(19.0, 14.25, 24.0, 12.0, 15.0)

Return type:

dict

morie.research.interracial_rates(offences, population)

Interracial offending rates against the random-mixing null (rate_per_offender_group, dyad_ratio).

Examples

>>> d = interracial_rates({"A_on_B": 120, "B_on_A": 200, "A_on_A": 900, "B_on_B": 300}, {"A": 80000, "B": 20000})
>>> [round(v, 12) for v in d["rate_per_pair_exposure"]]
[0.0075, 0.0125, 0.0140625, 0.075]

morie.research.judge_iv(z, d, y, weights=None, defier_share=(0, 0.05, 0.1), defier_effect=(0,))

Observational Wald ratio, the type shares under monotonicity, and the defier sensitivity.

Examples

>>> z = [0, 0, 0, 0, 1, 1, 1, 1]; d = [0, 0, 1, 0, 1, 1, 1, 0]; y = [1, 2, 5, 1, 6, 4, 5, 2]
>>> r = judge_iv(z, d, y)
>>> (round(r["first_stage"], 12), round(r["itt"], 12), round(r["wald"], 12))
(0.5, 2.0, 4.0)

Return type:

dict

morie.research.judge_iv_population(d0, d1, y0, y1, weights=None)

The exact decomposition of a binary-instrument design on a known population.

Examples

>>> r = judge_iv_population([0, 0, 1, 0, 1], [1, 1, 1, 0, 0], [0, 0, 1, 0, 0], [1, 0, 1, 1, 1])
>>> ({k: round(v, 12) for k, v in r["shares"].items()}, round(r["wald"], 12))
({'complier': 0.4, 'defier': 0.2, 'always': 0.2, 'never': 0.2}, 0.0)

Return type:

dict

morie.research.judge_slope_test(judge, d, y, weights=None, lo=None, hi=None)

The many-judge slope test of monotonicity (Frandsen, Lefgren & Leslie 2023).

Research.P14Slope.propensity_mono / outcome_diff / slope_bound / violation_refutes_monotonicity.

Examples

>>> judge = ["A"] * 6 + ["B"] * 6 + ["C"] * 6
>>> d = [0, 0, 0, 1, 1, 0,  0, 1, 1, 1, 0, 1,  1, 1, 1, 1, 1, 0]
>>> y = [1, 0, 0, 1, 0, 0,  0, 1, 1, 0, 0, 1,  1, 1, 1, 1, 0, 1]
>>> s = judge_slope_test(judge, d, y)
>>> (s["judges"]["judge"], [round(p, 12) for p in s["judges"]["propensity"]], s["pairs"]["violation"], s["violations"])
(['A', 'B', 'C'], [0.333333333333, 0.666666666667, 0.833333333333], [False, False, True], 1)
>>> [round(v, 12) for v in s["pairs"]["late"]]
[0.5, 1.0, 2.0]

Return type:

dict

morie.research.logit_rescale(beta, omitted_var, error_var=3.289868133696453)

Rescaling of logit coefficients across nested models (Research.P5.reduced_coefficient and kin).

Examples

>>> r = logit_rescale(beta=0.8, omitted_var=1)
>>> round(r["rescale"], 12)
0.875724044222

Return type:

dict

morie.research.logit_separation(y, x, intercept=True, method='auto')

Detect separation in a logistic regression and show why the fit cannot converge.

The exact check is the linear programme: maximise sum_i s_i x_i . d subject to 0 <= s_i x_i . d <= 1; a positive optimum is a separating direction (complete when every margin is positive, quasi-complete when some are zero). method="glm" inspects an iteratively reweighted fit for fitted probabilities at 0 or 1 instead, a heuristic.

Examples

>>> x = [[1, 0.2], [1, -1], [1, 0.5], [0, 0.1], [0, -0.4], [0, 1.2], [0, 0.3]]
>>> r = logit_separation([1, 1, 1, 0, 0, 0, 0], x)
>>> (r["separation"], all(v < 0 for v in r["loglik_along"].values()), r["loglik_along"]["t=100"] > -1e-6)
('complete', True, True)

Return type:

dict

morie.research.meta_dl_bias(variances, n_draws=2000, seed=0)

Size of the DerSimonian-Laird truncation bias under homogeneity (shared Philox stream).

Examples

>>> b = meta_dl_bias([0.010, 0.020, 0.015, 0.030, 0.012], n_draws=500, seed=1)
>>> (b["mean_tau2"] > 0, 0 < b["share_positive"] < 1)
(True, True)

Return type:

dict

morie.research.meta_hksj(estimates, variances, tau2='DL', level=0.95)

Hartung-Knapp-Sidik-Jonkman interval with DerSimonian-Laird or REML heterogeneity.

Research.P13HKSJ.hksj_wider_iff (wider than Wald iff q >= 1), Q_eq_zero_iff (degenerate iff every site equals the pooled value), hksj_equal_weights (equal weights: the one-sample t variance).

Examples

>>> h = meta_hksj([-0.25, -0.10, -0.40, 0.05, -0.30], [0.010, 0.020, 0.015, 0.030, 0.012])
>>> (round(h["tau2"], 12), round(h["q"], 12), round(h["hksj"]["se"], 12), round(h["wald"]["se"], 12), h["wider_than_wald"])
(0.006967741935, 1.086177559643, 0.070056783817, 0.067220197281, True)
>>> r = meta_hksj([-0.25, -0.10, -0.40, 0.05, -0.30], [0.010, 0.020, 0.015, 0.030, 0.012], tau2="REML")
>>> (round(r["tau2"], 12), round(r["q"], 12))
(0.002908035914, 1.271198964824)

Return type:

dict

morie.research.meta_random_effects(estimates, variances, level=0.95)

Fixed-effect and DerSimonian-Laird random-effects pooling, with what the truncation does.

Examples

>>> m = meta_random_effects([-0.25, -0.10, -0.40, 0.05, -0.30], [0.010, 0.020, 0.015, 0.030, 0.012])
>>> (round(m["tau2"], 12), m["variance_ratio"] >= 1, m["truncated"])
(0.006967741935, True, False)

Return type:

dict

morie.research.probability_of_necessity(p_treated, p_control)

Probability of necessity bounds (Research.P8.necessity_bounds, necessity_of_monotone).

Examples

>>> r = probability_of_necessity(p_treated=0.6, p_control=0.4)
>>> ({k: round(v, 12) for k, v in r["necessary_share_bounds"].items()}, round(r["pn_monotone"], 12))
({'lower': 0.2, 'upper': 0.6}, 0.333333333333)

Return type:

dict

morie.research.ranking_resolution(estimate, half_width)

Resolution of a ranking built from noisy risk scores (rank_reversal_exists, rank_stable_of_gap).

Examples

>>> r = ranking_resolution([0.2, 0.35, 0.8], half_width=[0.1, 0.1, 0.05])
>>> (r["n_pairs"], round(r["share_unidentified"], 12), r["resolution"])
(3, 0.333333333333, 0.2)

Return type:

dict

morie.research.recording_map(M, counts, names=None)

Recorded counts from true counts through a recording matrix.

M[i, j] is the share of true category-j offences recorded as category i.

Examples

>>> r = recording_map([[0.7, 0], [0.3, 1]], [100, 400], names=["robbery", "theft"])
>>> (r["recorded"], r["regime"], r["recorded_total"])
({'robbery': 70.0, 'theft': 430.0}, 'reclassification', 500.0)

Return type:

dict

morie.research.regression_to_mean(x1, x2, threshold, weights=None)

Observed, mirror and symmetrised change on the places selected for a high first-period count.

The symmetrised change is the same statistic on the data plus its period-swapped copy, an exchangeable population, so it is non-positive by Research.P19.selected_change_nonpos.

Examples

>>> r = regression_to_mean([5, 1, 7, 2, 9, 3], [3, 2, 6, 2, 4, 5], threshold=4)
>>> (r["n_selected"], round(r["selected_change"], 12), r["symmetrised_change"] <= 0)
(3, -2.666666666667, True)

Return type:

dict

morie.research.relative_risk_from_or(odds_ratio, base_rate=None, exposed_share=None)

Relative risk from an odds ratio: the interval arrest-only data allow (rr_between, or_overstates).

Examples

>>> relative_risk_from_or(3)["rr_bounds"]
{'lower': 1, 'upper': 3.0}
>>> round(relative_risk_from_or(3, base_rate=0.2, exposed_share=0.3)["risks"]["relative_risk"], 9)
2.333333333

Return type:

dict

morie.research.selective_labels(y, released, rule_released, weights=None)

What a release rule’s failure rate can be known from the cases a judge released.

Examples

>>> y = [0, 1, 0, 1, 1, 0]; released = [True, True, True, True, False, False]
>>> r = selective_labels(y, released, [True, True, True, False, False, False])
>>> (r["identified"], round(r["rule_rate"], 12))
(True, 0.333333333333)
>>> r2 = selective_labels(y, released, [True, True, True, True, True, False])
>>> ({k: round(v, 12) for k, v in r2["bounds"].items()}, round(r2["width"], 12))
({'lower': 0.4, 'upper': 0.6}, 0.2)

Return type:

dict

morie.research.sentence_effect_bounds(y, z, weights=None, contrast=None)

Worst-case identification bounds for a binary outcome under two sentences.

P(y=1, z=t) <= P[y(t)=1] <= P(y=1, z=t) + P(z != t), both ends attained (Research.P11.Pop.outcome_bounds) the contrast interval has width one and contains zero (ate_width_one, ate_contains_zero).

Examples

>>> b = sentence_effect_bounds([1, 0, 1, 0, 1, 1], ["a", "a", "a", "b", "b", "b"])
>>> (round(b["ate_width"], 12), round(b["ate_bounds"]["lower"], 12), round(b["ate_bounds"]["upper"], 12))
(1.0, -0.5, 0.5)

Return type:

dict

morie.research.sentence_effect_mtr(y, z, weights=None, contrast=None, direction='non-decreasing')

Monotone-treatment-response bounds (Research.P11.Pop.mtr_lower, mtr_upper, mtr_upper_attained).

Examples

>>> r = sentence_effect_mtr([1, 0, 1, 0, 1, 1], ["a", "a", "a", "b", "b", "b"])
>>> (r["bounds"]["lower"], round(r["bounds"]["upper"], 12))
(0, 0.5)

Return type:

dict

morie.research.sentence_effect_mts(y, z, weights=None, contrast=None)

Monotone-treatment-selection bounds for a sentencing contrast (Manski & Pepper 2000).

The observed mean of the harsher group bounds E[y(b)] from above and the lighter group’s bounds E[y(a)] from below (Research.P11.Pop.mts_mean_b_le, mts_mean_a_ge); the naive difference overstates the effect (mts_ate_le_naive) with MTR as well the contrast lies in [0, naive] (mtr_mts_bounds).

Examples

>>> r = sentence_effect_mts([1, 0, 1, 0, 1, 1], ["a", "a", "a", "b", "b", "b"])
>>> ({k: round(v, 12) for k, v in r["ate_bounds_mts"].items()}, r["ate_bounds_mtr_mts"])
({'lower': -0.5, 'upper': 0.0}, {'lower': 0, 'upper': 0.0})

Return type:

dict

morie.research.shrinkage_loss(theta, noise, weights=None, B=None)

Loss of the shrinkage estimator against a known truth, next to the closed form of loss_eq.

Research.P19Shrinkage.loss_eq / loss_min / loss_bstar_eq / loss_bstar_le_raw.

Examples

>>> l = shrinkage_loss([10, 20, 30, 40], [3, -3, -3, 3], B=0.5)
>>> (l["loss"], l["closed_form"], round(l["B_star"], 12), round(l["loss_star"], 12), l["loss_raw"], l["noise_law"])
(134.0, 134.0, 0.067164179104, 33.582089552239, 36.0, True)

Return type:

dict

morie.research.spillover_effects(y, exposure, stratum, weights=None)

Direct, spillover and net effects under a stated exposure mapping (stratified exposure adjustment).

Examples

>>> y = [5, 4.6, 4, 6, 5.6, 5]; e = [0, 1, 2, 0, 1, 2]; s = ["a", "a", "a", "b", "b", "b"]
>>> r = spillover_effects(y, e, s)
>>> (round(r["spillover"], 12), round(r["direct"], 12), r["positivity"])
(-0.4, -0.6, True)

Return type:

dict

morie.research.spillover_exposure(treated, edges)

Three-level exposure from a treatment on a place network (edges are 1-based pairs, as in R).

Examples

>>> e = spillover_exposure([1, 0, 0, 0, 1], [[1, 2], [2, 3], [3, 4], [4, 5]])
>>> (list(e["exposure"]), list(e["treated_neighbours"]))
([2, 1, 0, 1, 2], [0, 1, 0, 1, 0])

morie.research.spillover_exposure_probs(n, edges, n_treated, n_draws=5000, exact_max=20000, joint=False, seed=0)

Exposure probabilities of a completely randomised deployment (exact enumeration, else Philox Monte Carlo).

Examples

>>> pr = spillover_exposure_probs(4, [[1, 2], [2, 3], [3, 4]], 1)
>>> [round(float(v), 12) for v in pr[:, 2]]
[0.25, 0.25, 0.25, 0.25]

morie.research.spillover_ht(y, exposure, probs, joint=None)

Horvitz-Thompson totals and contrasts under a randomised deployment (ht_unbiased, ht_variance).

Examples

>>> pr = spillover_exposure_probs(6, [[1, 2], [2, 3], [3, 4], [4, 5], [5, 6]], 2, joint=True)
>>> e = spillover_exposure([0, 1, 0, 0, 0, 1], [[1, 2], [2, 3], [3, 4], [4, 5], [5, 6]])["exposure"]
>>> r = spillover_ht([5, 3.5, 4.5, 5, 4.5, 3.5], list(e), pr["marginal"], pr["joint"])
>>> (round(r["total_effect"], 12), sorted(r.keys())[:3])
(-0.666666666667, ['direct', 'means', 'se'])

Return type:

dict

morie.research.spillover_ht_variance(y_pot, pi, pij, form='ht')

Design variance of the Horvitz-Thompson total (ht_variance; SYG form, syg_eq_ht).

Examples

>>> pr = spillover_exposure_probs(6, [[1, 2], [2, 3], [3, 4], [4, 5], [5, 6]], 2, joint=True)
>>> v = spillover_ht_variance([1] * 6, pr["marginal"][:, 2], pr["joint"][:, :, 2], form="both")
>>> (v["fixed_size"], round(v["level_count"], 12))
(True, 2.0)

morie.research.two_point_bound(p, q, theta_p, theta_q, estimator=None)

Le Cam’s two-point lower bound on a finite sample space.

Examples

>>> p = [0.0625, 0.25, 0.375, 0.25, 0.0625]
>>> q = [0.0256, 0.1536, 0.3456, 0.3456, 0.1296]
>>> b = two_point_bound(p, q, theta_p=2, theta_q=3, estimator=[0, 1.25, 2.5, 3.75, 5])
>>> (round(b["tv"], 12), round(b["bound"], 12), round(b["minimax_risk"], 12), b["satisfied"], round(b["sum_min"], 12))
(0.1627, 0.41865, 1.125, True, 0.8373)

Return type:

dict