Bayesian Inference

Bayesian Statistics

Hubble constant posterior evaluation.

BayesianStatistics loads saved Cramér-Rao bounds and orchestrates the full Hubble-constant posterior evaluation using the real GLADE galaxy catalog, simulation-based SimulationDetectionProbability, full Fisher-matrix covariance, and multiprocessing.

Invoked via main.py:evaluate() / --evaluate CLI flag. Output is written to simulations/posteriors/ as JSON.

class master_thesis_code.bayesian_inference.bayesian_statistics.BayesianStatistics[source]

Bases: object

Hubble constant posterior evaluation.

Loads saved Cramér-Rao bounds from CSV, constructs a simulation-based SimulationDetectionProbability, builds multivariate-normal GW likelihoods from the full Fisher-matrix covariance, and evaluates per-detection posteriors over an H₀ grid using a multiprocessing pool.

Invoked via main.py:evaluate() (--evaluate CLI flag). Output is written to simulations/posteriors/ as JSON.

Omega_DE: float
Omega_m: float
additional_galaxies_without_bh_mass: dict[str, dict[str, list[float]]]
cosmological_model: LamCDMScenario
cramer_rao_bounds: DataFrame
detection: Detection
evaluate(galaxy_catalog, cosmological_model, h_value, num_workers=None, catalog_only=False, pdet_dl_bins=60, pdet_mass_bins=40, pdet_estimator='local_linear', fisher_cond_threshold=1e16, normalization_mode='generator_marginal', base_seed=0, allow_low_pdet_coverage=False, h_values=None, smear_global_selection=False, dgen_catalog_selection='4d_exact', pdet_z_resolved=True, pdet_wbh_z_resolved=False, host_z_kernel='auto', host_mass_kernel='auto')[source]
Parameters:

  • galaxy_catalog (GalaxyCatalogueHandler)

  • cosmological_model (Model1CrossCheck)

  • h_value (float)

  • num_workers (int | None)

  • catalog_only (bool)

  • pdet_dl_bins (int)

  • pdet_mass_bins (int)

  • pdet_estimator (str)

  • fisher_cond_threshold (float)

  • normalization_mode (str)

  • base_seed (int)

  • allow_low_pdet_coverage (bool)

  • h_values (Sequence[float] | None)

  • smear_global_selection (bool)

  • dgen_catalog_selection (str)

  • pdet_z_resolved (bool)

  • pdet_wbh_z_resolved (bool)

  • host_z_kernel (str)

  • host_mass_kernel (str)

Return type:

None

galaxy_weights: dict[str, dict[str, list[float]]]
h: float
h_values: list[float]
h_values_with_bh_mass: list[float]
p_D(galaxy_catalog, redshift_upper_limit, pool, completeness, detection_probability_obj)[source]
Parameters:

  • galaxy_catalog (GalaxyCatalogueHandler)

  • redshift_upper_limit (float)

  • pool (Pool)

  • completeness (CompletenessModel)

  • detection_probability_obj (SimulationDetectionProbability)

Return type:

None

p_Di(possible_host_galaxies, possible_host_galaxies_with_bh_mass, detection_index, pool, completeness, detection_probability_obj, redshift_upper_limit=HOST_DRAW_Z_MAX)[source]
Parameters:

  • possible_host_galaxies (list[HostGalaxy])

  • possible_host_galaxies_with_bh_mass (list[HostGalaxy])

  • detection_index (int)

  • pool (Pool)

  • completeness (CompletenessModel)

  • detection_probability_obj (SimulationDetectionProbability)

  • redshift_upper_limit (float)

Return type:

tuple[float, float]

posterior_data: dict[int, list[float]]
posterior_data_with_bh_mass: dict[int | str, Any]
w_0: float
w_a: float
master_thesis_code.bayesian_inference.bayesian_statistics.child_process_init(redshift_lower_limit, redshift_upper_limit, bh_mass_lower_limit, bh_mass_upper_limit, current_detection_probability, current_means_3d, current_cov_inv_3d, current_log_norm_3d, current_means_4d, current_cov_inv_4d, current_log_norm_4d, current_det_index_to_slot, current_sigma2_cond_arr, current_proj_arr, current_det_d_L_arr, current_det_d_L_unc_arr, current_det_M_arr, current_det_phi_arr, current_det_theta_arr, current_D_h_table=None)[source]
Parameters:

  • redshift_lower_limit (float)

  • redshift_upper_limit (float)

  • bh_mass_lower_limit (float)

  • bh_mass_upper_limit (float)

  • current_detection_probability (SimulationDetectionProbability)

  • current_means_3d (ndarray[tuple[Any, ...], dtype[float64]])

  • current_cov_inv_3d (ndarray[tuple[Any, ...], dtype[float64]])

  • current_log_norm_3d (ndarray[tuple[Any, ...], dtype[float64]])

  • current_means_4d (ndarray[tuple[Any, ...], dtype[float64]])

  • current_cov_inv_4d (ndarray[tuple[Any, ...], dtype[float64]])

  • current_log_norm_4d (ndarray[tuple[Any, ...], dtype[float64]])

  • current_det_index_to_slot (dict[int, int])

  • current_sigma2_cond_arr (ndarray[tuple[Any, ...], dtype[float64]])

  • current_proj_arr (ndarray[tuple[Any, ...], dtype[float64]])

  • current_det_d_L_arr (ndarray[tuple[Any, ...], dtype[float64]])

  • current_det_d_L_unc_arr (ndarray[tuple[Any, ...], dtype[float64]])

  • current_det_M_arr (ndarray[tuple[Any, ...], dtype[float64]])

  • current_det_phi_arr (ndarray[tuple[Any, ...], dtype[float64]])

  • current_det_theta_arr (ndarray[tuple[Any, ...], dtype[float64]])

  • current_D_h_table (dict[float, float] | None)

Return type:

None

master_thesis_code.bayesian_inference.bayesian_statistics.compute_catalog_draw_weight_total(galaxy_catalog, z_max=HOST_DRAW_Z_MAX)[source]

Total draw-eligible catalogue rate weight W_cat (h-independent scalar).

\[W_\mathrm{cat} = \sum_{g:\, z_g < z_\mathrm{max}} w_g, \qquad w_g = \frac{R_\mathrm{eff}(M_g)}{1 + z_g},\]

over the SAME pruned catalogue rows and the SAME eligibility mask (z_g < HOST_DRAW_Z_MAX, no other cut) that the generator’s in-catalogue host draw normalizes over — this is exactly total_weight in draw_rate_weighted_hosts(). It is the draw-side companion normalizer of the completeness-weighted population volume precompute_completeness_population_volume(); their ratio n_hat_w = W_cat / V_f(h) is the generator-consistent rate-weight density that replaces the Option-A calibration n_bar_w = Sigma_glob/beta_G in the generator_marginal normalization mode.

W_cat carries NO P_det and NO h dependence: it normalizes the draw, not the detection (domain note in the derivation packet §3.2). Any analysis-depth cap (issue-#30 z_max_cap) must be applied to z_max HERE and in V_f together with the candidate window (f29a5e7 principle: numerator and denominator move together).

Parameters:

  • galaxy_catalog (GalaxyCatalogueHandler) – Loaded catalogue handler (its reduced_galaxy_catalog is summed over; same rows the rate-weighted draw uses).

  • z_max (float) – Exclusive upper redshift bound of the draw eligibility. Defaults to HOST_DRAW_Z_MAX.

Returns:

W_cat in yr^-1 (the emri_rate.C_NORM scale cancels in every ratio it enters).

Return type:

float

References

  • results/lcat_h_dependence_20260725/DERIVATION_GENERATOR_CONSISTENT_NORM.md §2.3 Eq. (4) (spec; W_cat anchor: 6.3477e8 over 9,060,017 pruned rows).

  • master_thesis_code/galaxy_catalogue/handler.py, draw_rate_weighted_hosts (the generator draw this normalizer replicates).

  • Babak et al. (2017), arXiv:1703.09722 — per-MBH rate R_eff.

master_thesis_code.bayesian_inference.bayesian_statistics.compute_sigma_deviation(sigma, sigma_error, h_mean, h_mean_error)[source]
Parameters:

  • sigma (float)

  • sigma_error (float)

  • h_mean (float)

  • h_mean_error (float)

Return type:

tuple[float, float]

master_thesis_code.bayesian_inference.bayesian_statistics.eddington_shifted_host_mass(host_M, host_M_error)[source]

Effective host mass under the rate-weighted (Eddington-in-M) prior.

The per-galaxy mass prior N(M; M_g, sigma_M^2) * R_eff(M) / Z_M is, under a locally log-linear R_eff (exponential-tilt identity), EXACTLY the shifted Gaussian N(M; M_g (1 + alpha sigma_rel^2), sigma_M^2) with alpha = dln R_eff / dln M |_{M_g} and sigma_rel = sigma_M / M_g. Classic Eddington (1913) correction; derivation and curvature-residual control in docs/derivations/G2d_host_mass_rate_prior.md (G7 row 9).

Parameters:

  • host_M (float) – Catalogue (source-frame) host BH mass estimate [M_sun].

  • host_M_error (float) – 1-sigma mass uncertainty [M_sun].

Returns:

The shifted effective mass M_g^eff [M_sun]; equals host_M when the uncertainty is zero/invalid (bare-Gaussian limit).

Return type:

float

master_thesis_code.bayesian_inference.bayesian_statistics.precompute_completeness_population_volume(h_values, completeness, *, z_min=1e-6, z_max=HOST_DRAW_Z_MAX, n_grid=4096)[source]

Completeness-weighted population volume V_f(h) (per steradian).

\[V_f(h) = \int_{z_\mathrm{min}}^{z_\mathrm{max}} \bar f(z, h)\, \frac{1}{1+z}\,\frac{dV_c}{dz\,d\Omega}\, dz ,\]

the SAME integral the generator’s in-catalogue mixture fraction F uses (F = V_f / V_tot, master_thesis_code.dark_siren_injection.compute_global_catalog_fraction()), with the SAME sky-averaged completeness f_bar and the SAME trapezoid quadrature convention (linspace grid, _DEFAULT_Z_GRID_POINTS = 4096 nodes there; the completeness curve is piecewise linear, so the trapezoid rule is exact-in-shape and more robust than Gauss-Legendre on the kinks). Under the frozen pixel m_th map f_bar is h-invariant, so V_f(h) = V_f(0.73) (0.73/h)^3 exactly (dV_c/dz carries the h^-3); the table is still evaluated per-h so a future h-dependent completeness model flows through unchanged (derivation packet §2.2 note).

Role: n_hat_w(h) = W_cat / V_f(h) is the generator draw-side rate-weight density of the generator_marginal mode — the calibration contains NO P_det, so the Option-A identity Sigma_glob = n_hat_w * beta_G is never invoked (the whole point of the mode).

Parameters:

  • h_values (list[float]) – Hubble parameter values to evaluate.

  • completeness (CompletenessModel) – Catalogue completeness model (f_bar accessor), the SAME frozen-cache object the generator uses (C1 consistency).

  • z_min (float) – Lower integration bound (matches the generator’s _DEFAULT_Z_MIN = 1e-6).

  • z_max (float) – Upper integration bound — the DRAW depth HOST_DRAW_Z_MAX, NOT the detection horizon z_max(h) (domain note, derivation §3.2). An issue-#30 depth cap must move this together with W_cat.

  • n_grid (int) – Trapezoid grid nodes (default 4096, the generator convention).

Returns:

Dict mapping h -> V_f(h) in Mpc^3 sr^-1 (same measure as D(h)/beta_Gbar(h)). Validation anchor: V_f(0.73) = 2.3237e8 (generator_norm_Vf_tables.json).

Return type:

dict[float, float]

References

  • results/lcat_h_dependence_20260725/DERIVATION_GENERATOR_CONSISTENT_NORM.md §2.2-§2.3, Eq. (4) (spec and numeric anchors).

  • master_thesis_code/dark_siren_injection.py, compute_global_catalog_fraction — the generator integral replicated.

  • Gray et al. (2020), arXiv:1908.06050, Eq. (9) — completeness f(z).

master_thesis_code.bayesian_inference.bayesian_statistics.precompute_completion_denominator(h_values, detection_probability_obj, Omega_m, Omega_DE, *, completeness=None, quad_n=_DH_QUAD_ORDER, z_max_cap=None)[source]

Precompute the completion-term denominator D(h) for each h value.

Gray et al. (2020), arXiv:1908.06050, Eqs. 33 / A.19: the out-of-catalogue selection denominator integrates the detection probability against the EMRI population prior over the detectable volume.

\[D(h) = \int_{z_{\min}}^{z_{\max}(h)} P_{\det}(d_L(z,h)) \,\frac{1}{1+z}\,\frac{dV_c}{dz\,d\Omega}\, dz\]

where z_max(h) is the redshift corresponding to the P_det grid’s maximum d_L at the given h, and 1/(1+z) is the source-to-detector time dilation (matching comp_num and the event sampler master_thesis_code.emri_rate.p_pop_unnormalized()).

Role in the partition-norm likelihood:

D(h) is the FULL-volume selection normalisation D(h) = beta_G(h) + beta_Gbar(h) – the denominator of the single per-event ratio p_i = (beta_G L_cat + B_num) / D(h) (p_Di()). It carries no (1-f) factor: the incompleteness lives in its missing-volume partner precompute_missing_completion_denominator() (beta_Gbar = INTEGRAL (1-f) P_det dVc/(1+z)), and the in-catalogue share is recovered by beta_G = D(h) - beta_Gbar. The selection-weighted catalog membership weight w_G = beta_G/D(h) = beta_G/(beta_G+beta_Gbar) (Gray Eq. 29) is now computed EXACTLY – it replaced the earlier scalar narrow-window approximation completeness(z_det).

Modeling assumption (still in force): constant comoving number density

for the missing galaxies – the galaxy number density n_gal(z) and the mass-integrated rate INTEGRAL dM R_EMRI(z,M) are taken z-independent (the latter exact under the p0=1 surrogate), so they are overall constants that cancel between the discrete catalogue sums and the continuous integrals (Option A; see precompute_global_catalog_selection()). Departures (clustering, rate/MF evolution) are second order.

Parameters:

  • h_values (list[float]) – List of Hubble parameter values to evaluate.

  • detection_probability_obj (SimulationDetectionProbability) – SimulationDetectionProbability instance (must have get_dl_max and detection_probability_without_bh_mass_interpolated_zero_fill).

  • Omega_m (float) – Matter density parameter.

  • Omega_DE (float) – Dark energy density parameter.

  • quad_n (int) – Gauss-Legendre quadrature order (default 100).

  • completeness (CompletenessModel | None)

  • z_max_cap (float | None)

Returns:

Dict mapping h -> D(h) in units of Mpc^3/sr.

Return type:

dict[float, float]

master_thesis_code.bayesian_inference.bayesian_statistics.precompute_global_catalog_selection(h_values, galaxy_catalog, detection_probability_obj, *, with_bh_mass, z_max_cap=None, smear_sigma_z=False)[source]

Precompute the GLOBAL in-catalogue selection denominator (Option A).

The partition-norm restructure forms the in-catalogue likelihood as L_cat = (sum_local w_g N_g) / (sum_global w_g D_g) where the SELECTION denominator runs over the FULL catalogue out to the detection horizon z_max(h), NOT the per-event candidate ball. Globalising the denominator makes L_cat scale-free, so the per-galaxy <-> per-volume number-density factor n_gal cancels against the continuous beta_G(h) = D(h) - beta_Gbar(h) and no calibration constant is needed (Gray et al. 2020, arXiv:1908.06050, Eq. 29: the discrete catalogue sum is the Monte-Carlo realisation of beta_G = INTEGRAL f P_det dVc/(1+z)).

\[\Sigma_{\mathrm{global}}(h) = \sum_{g:\, z_g < z_{\max}(h)} w_g\, P_{\det}\bigl(d_L(z_g, h)\bigr), \qquad w_g = \frac{R_\mathrm{eff}(M_g)}{1 + z_g}.\]

The weight w_g is IDENTICAL to the rate-weighted host draw (draw_rate_weighted_hosts()) and the in-catalogue likelihood weight (_rate_weight()). P_det is evaluated SKY-MARGINALISED (phi = theta = 0), on the same footing as the completion D(h) / beta_Gbar (the per-galaxy sky dependence is deferred to the pixelated-completeness change). D_g ~= P_det(z_g) uses the narrow galaxy-redshift-PDF limit. The sum is event-INDEPENDENT, so it is precomputed once per h like D(h).

Parameters:

  • h_values (list[float]) – Hubble parameter values to evaluate.

  • galaxy_catalog (GalaxyCatalogueHandler) – Loaded catalogue handler (its reduced_galaxy_catalog is summed over; same rows the rate-weighted draw uses).

  • detection_probability_obj (SimulationDetectionProbability) – Detection probability (provides get_dl_max and the 3D / 4D P_det accessors).

  • with_bh_mass (bool) – False uses the 3D (sky+distance) P_det (the without-BH-mass channel); True uses the 4D (sky+distance+observer-frame mass M_z = M_g(1+z_g)) P_det, the global companion of the with-BH-mass catalogue sum.

  • z_max_cap (float | None)

  • smear_sigma_z (bool)

Returns:

Dict mapping h -> sum_global w_g D_g(h) (dimensionless rate-weighted detection count).

Return type:

dict[float, float]

References

Gray et al. (2020), arXiv:1908.06050, Eq. (29) — beta_G selection

integral (here its discrete catalogue realisation).

Babak et al. (2017), arXiv:1703.09722 — per-MBH rate R_eff

(master_thesis_code.emri_rate.R_eff_per_mbh()).

master_thesis_code.bayesian_inference.bayesian_statistics.precompute_missing_completion_denominator(h_values, detection_probability_obj, completeness, *, quad_n=_DH_QUAD_ORDER, z_max_cap=None)[source]

Precompute the missing-volume selection integral beta_Gbar(h).

The (1-f(z)) companion of precompute_completion_denominator() (which returns the unchanged full-volume D(h) = beta_G + beta_Gbar). Gray et al. (2020), arXiv:1908.06050, Eq. (33): the out-of-catalogue selection integral weights the full detection denominator by the incompleteness 1 - f(z), i.e. it integrates only over the galaxies the catalogue is missing:

\[\beta_{\bar G}(h) = \int_{z_{\min}}^{z_{\max}(h)} \bigl(1 - f(z)\bigr)\, P_{\det}(d_L(z,h))\,\frac{1}{1+z}\,\frac{dV_c}{dz}\, dz .\]

The in-catalogue selection normalisation is then beta_G(h) = D(h) - beta_Gbar(h) = INTEGRAL f(z) P_det (1/(1+z)) dVc. f(z) = completeness.get_completeness_at_redshift(z, h) is the SAME completeness call the generator uses (master_thesis_code.dark_siren_injection.compute_global_catalog_fraction() and _draw_dark_redshifts), so the inference completion population and the injected dark population are bit-for-bit identical.

Parameters:

  • h_values (list[float]) – Hubble parameter values to evaluate.

  • detection_probability_obj (SimulationDetectionProbability) – Same object passed to precompute_completion_denominator() (provides get_dl_max and detection_probability_without_bh_mass_interpolated_zero_fill).

  • completeness (CompletenessModel) – Catalogue completeness f(z) (Gray Eq. 9). Evaluated sky-marginalised, identically to the generator.

  • quad_n (int) – Gauss-Legendre quadrature order (default _DH_QUAD_ORDER), matching D(h).

  • z_max_cap (float | None)

Returns:

Dict mapping h -> beta_Gbar(h) in units of Mpc^3/sr (same as D(h)).

Return type:

dict[float, float]

References

Gray et al. (2020), arXiv:1908.06050, Eq. (33) — out-of-catalogue

selection denominator (here the missing (1-f) fraction).

master_thesis_code.bayesian_inference.bayesian_statistics.resolve_host_mass_kernel(host_mass_kernel, normalization_mode, host_z_kernel)[source]

Resolve the 2D host-mass kernel selection to ‘gaussian’ or ‘trunc_lognormal’.

Decomposition flag for the #40 remainder (RATIFY-M3/M4, docs/derivations/mass_marginal_2d_kernel.md): makes the truncated lognormal x R_eff mass kernel separately selectable from the normalization leg. "auto" reproduces the historical bundling — the truncated kernel if and only if normalization_mode == "mass_trunc" — so the production default path is unchanged.

Guard (derivation §3.3): the delta-kernel (point) host-z numerator path always evaluates the analytic Gaussian mass product at the catalogue host_M (issue #24 point-M treatment), while the trunc_lognormal denominator carries the LN x R_eff prior — N_g and D_g would silently use DIFFERENT mass priors, violating the counted-once-in-M invariant. That combination raises instead of running silently.

Parameters:

  • host_mass_kernel (str) – One of HOST_MASS_KERNEL_CHOICES.

  • normalization_mode (str) – The in-catalogue normalization mode (see p_Di).

  • host_z_kernel (str) – The (unresolved) numerator host-z kernel selection; resolved internally via resolve_host_z_kernel() for the prior-consistency guard.

Returns:

"gaussian" (analytic Gaussian mass product + G2d moment-matched shift in the calibrated kernels) or "trunc_lognormal" (the ratified truncated lognormal x R_eff kernel, GH numerator with small-sigma crossover + GL-in-lnM denominator).

Raises:

ValueError – Unknown choice, or the prior-inconsistent combination of a point host-z numerator with the trunc_lognormal mass kernel.

Return type:

str

master_thesis_code.bayesian_inference.bayesian_statistics.resolve_host_z_kernel(host_z_kernel, normalization_mode)[source]

Resolve the numerator host-z kernel selection to ‘point’ or ‘volume_deconv’.

Decomposition flag for issue #40(a): makes the delta-kernel (point/point) in-catalogue numerator separately selectable from the normalization leg. "auto" reproduces the historical bundling — the delta kernel if and only if normalization_mode == "generator_marginal" — so the production default path is unchanged. Explicit "point" / "volume_deconv" override the numerator kernel only; the selection-normalization machinery (n_hat_w/D_gen vs n_bar_w/D) remains governed by normalization_mode.

Parameters:

  • host_z_kernel (str) – One of HOST_Z_KERNEL_CHOICES.

  • normalization_mode (str) – The in-catalogue normalization mode (see p_Di).

Returns:

"point" (delta kernel at the catalogue z_g) or "volume_deconv" (the mode’s own quadrature kernel — volume-deconvolved in the *_marginal modes, bare Gaussian in “global”/”local_ratio”).

Return type:

str

master_thesis_code.bayesian_inference.bayesian_statistics.single_host_likelihood(host_phiS, host_qS, host_z, host_z_error, host_M, host_M_error, detection_index, h, evaluate_with_bh_mass, normalization_mode='generator_marginal', base_seed=0, host_z_kernel='auto', host_mass_kernel='auto')[source]
Parameters:

  • host_phiS (float)

  • host_qS (float)

  • host_z (float)

  • host_z_error (float)

  • host_M (float)

  • host_M_error (float)

  • detection_index (int)

  • h (float)

  • evaluate_with_bh_mass (bool)

  • normalization_mode (str)

  • base_seed (int)

  • host_z_kernel (str)

  • host_mass_kernel (str)

Return type:

list[float]

master_thesis_code.bayesian_inference.bayesian_statistics.single_host_likelihood_batch(host_phiS, host_qS, host_z, host_z_error, host_M, host_M_error, detection_index, h, evaluate_with_bh_mass, normalization_mode='generator_marginal', host_z_kernel='auto', host_mass_kernel='auto')[source]

Host-batched twin of single_host_likelihood().

Computes the per-host likelihood integrals for n candidate hosts of one detection in a single vectorized pass. Row i of the result equals single_host_likelihood(...) called with host i’s scalars — the same physics, the same quadrature (fixed_quad’s exact affine node map and reduction, see _batched_gl_nodes()/_batched_gl_reduce()), the same Gaussian pdf operation order (_gaussian_pdf()) — with the host loop moved from Python/starmap into the array axis. Eliminated per-host costs: scipy.stats.norm frozen-distribution construction, the event-level dist_to_redshift window calls (now once per batch), and per-host p_det interpolator calls (now one call over all hosts’ nodes).

Reads the child_process_init worker globals (the subset the scalar kernel actually uses). base_seed is intentionally absent: it was a dead parameter of the scalar signature (vestigial from the removed MC denominator).

Parameters:

  • host_phiS (ndarray[tuple[Any, ...], dtype[float64]]) – Host ecliptic azimuths, shape (n,).

  • host_qS (ndarray[tuple[Any, ...], dtype[float64]]) – Host ecliptic polar angles, shape (n,).

  • host_z (ndarray[tuple[Any, ...], dtype[float64]]) – Host redshifts, shape (n,).

  • host_z_error (ndarray[tuple[Any, ...], dtype[float64]]) – Host redshift 1-sigma errors, shape (n,).

  • host_M (ndarray[tuple[Any, ...], dtype[float64]]) – Host BH masses [M_sun], shape (n,).

  • host_M_error (ndarray[tuple[Any, ...], dtype[float64]]) – Host BH mass 1-sigma errors, shape (n,).

  • detection_index (int) – CRB row index of the detection.

  • h (float) – Dimensionless Hubble parameter.

  • evaluate_with_bh_mass (bool) – Include the with-BH-mass channel.

  • normalization_mode (str) – In-catalogue normalization mode (see p_Di).

  • host_z_kernel (str) – Numerator host-z kernel selection (issue #40a); "auto" reproduces the historical mode bundling.

  • host_mass_kernel (str) – 2D host-mass kernel selection (#40 remainder); "auto" reproduces the historical mass_trunc bundling.

Returns:

Array of shape (n, 6) when evaluate_with_bh_mass else (n, 4); columns match the scalar kernel’s return list.

Return type:

ndarray[tuple[Any, …], dtype[float64]]

master_thesis_code.bayesian_inference.bayesian_statistics.single_host_likelihood_integration_testing(possible_host, detection, detection_index, h, evaluate_with_bh_mass)[source]
Parameters:

  • possible_host (HostGalaxy)

  • detection (Detection)

  • detection_index (int)

  • h (float)

  • evaluate_with_bh_mass (bool)

Return type:

list[float]

master_thesis_code.bayesian_inference.bayesian_statistics.use_detection(detection)[source]
Parameters:

detection (Detection)

Return type:

bool

master_thesis_code.bayesian_inference.bayesian_statistics.weighted_ratio_of_sums(numerators, denominators, weights)[source]

Weighted in-catalog ratio-of-sums likelihood (Σ w·N) / (Σ w·D).

Generalizes the equal-weight Gray et al. (2020) in-catalog term L_cat = (Σ_g N_g) / (Σ_g D_g) (Eq. A.9/A.10) by weighting each candidate host galaxy g by an astrophysical rate prior w(g):

\[L_\mathrm{cat} = \frac{\sum_g w(g)\,N_g}{\sum_g w(g)\,D_g}.\]

The weight enters numerator and denominator identically, so

  • any overall rescaling of w cancels (SCALING INVARIANCE), and

  • constant weights reproduce the plain ratio of sums exactly (the equal-weight Change-2 limit).

This is the inference-side counterpart of the rate-weighted host draw draw_rate_weighted_hosts().

Parameters:

  • numerators (Sequence[float]) – Per-host likelihood numerators N_g (host-aligned).

  • denominators (Sequence[float]) – Per-host selection denominators D_g (host-aligned, same order as numerators).

  • weights (Sequence[float]) – Per-host rate weights w(g) (host-aligned, same order as numerators / denominators).

Returns:

The weighted ratio of sums, or 0.0 when the weighted denominator Σ w·D is non-positive (matching the unweighted guard).

Return type:

float

References

Gray et al. (2020), arXiv:1908.06050, Eqs. (A.9)/(A.10) — in-catalog

ratio-of-sums likelihood, here weighted by a galaxy rate prior.

master_thesis_code.bayesian_inference.bayesian_statistics.weighted_sum(values, weights)[source]

Rate-weighted sum Σ_g w_g · v_g (the in-catalogue numerator building block).

The partition-norm in-catalogue likelihood is L_cat = (Σ_local w_g N_g) / (Σ_global w_g D_g) (Gray et al. 2020, arXiv:1908.06050, Eqs. A.10 / 29) where the GW-likelihood numerator sum runs over the local candidate ball but the SELECTION denominator runs over the full catalogue (precompute_global_catalog_selection()). This helper returns the weighted sum of either; an empty input yields 0.0.

Parameters:

  • values (Sequence[float]) – Per-host values v_g (host-aligned).

  • weights (Sequence[float]) – Per-host rate weights w_g (same order as values).

Returns:

Σ_g w_g · v_g (0.0 for empty inputs).

Return type:

float

Simulation Detection Probability

Simulation-based detection probability from injection campaign data.

Replaces DetectionProbability (KDE-based) with a detection-horizon survival-function approach that loads raw injection CSVs and applies an SNR threshold at evaluation time.

Detection is deterministic optimal SNR (SNR = sqrt(<h|h>), no detector noise; the threshold is passed as snr_threshold). Because the GW strain amplitude scales as 1/d_L, each injection has an h-invariant detection horizon

d_hor_k = SNR_k * d_L_k / snr_threshold [Gpc]

(the d_L at which event k would sit exactly at threshold), and the detection probability is exactly the survival function of the horizon distribution:

p_det(d_L) = P(d_hor >= d_L) = fraction of injections detectable at d_L.

The horizon set is independent of the trial Hubble parameter h (the 1/d_L amplitude scaling and the d_L rescaling cancel), so the survival grid is built once and reused for every h.

All injection data is pooled regardless of the Hubble parameter value used during the injection campaign. The legacy SNR-rescaling helper _rescale_snr() is retained (it is exact physics and is unit-tested directly) but the survival estimator does not require it.

References

  • Finn & Chernoff (1993), arXiv:gr-qc/9301003 — detection horizon / SNR threshold for inspirals.

  • Finn (1996), arXiv:gr-qc/9601048 — p_det = P(Theta > Theta_thr) as a survival function of the orientation/distance factor.

  • Gray et al. (2020), arXiv:1908.06050, Section III.B-C — selection function structure for the numerator/denominator.

  • Mandel, Farr & Gair (2019), arXiv:1809.02063 — fixed-injection selection function evaluated at hypothesis observer-frame parameters.

  • SNR ~ 1/d_L: Hogg (1999), arXiv:astro-ph/9905116, Eq. (16).

class master_thesis_code.bayesian_inference.simulation_detection_probability.SimulationDetectionProbability(injection_data_dir, snr_threshold, h_grid=None, *, dl_bins=_DEFAULT_DL_BINS, mass_bins=_DEFAULT_M_BINS, h_prior_range=(_DEFAULT_H_PRIOR_MIN, _DEFAULT_H_PRIOR_MAX), bandwidth_scale=1.0, estimator=_DEFAULT_ESTIMATOR, n_sky_bands=_DEFAULT_N_SKY_BANDS, _force_unit_weights=False, expected_z_max=None, allow_shallow_pool=False, pdet_z_resolved=False, pdet_wbh_z_resolved=False)[source]

Bases: object

Simulation-based detection probability from injection campaign data.

Loads raw injection CSVs (z, M, phiS, qS, SNR, h_inj, luminosity_distance), pools ALL events regardless of h_inj, and builds the detection-horizon survival grids once (they are h-invariant).

For a source with measured optimal SNR_raw at luminosity distance d_L_k, the detection horizon is the distance at which it would sit exactly at threshold:

d_hor_k = SNR_raw_k * d_L_k / snr_threshold [Gpc]

Detection is deterministic (SNR >= threshold ⇔ d_L <= d_hor), so the detection probability is the survival function of the horizon distribution:

p_det(d_L) = P(d_hor >= d_L).

The horizon set is independent of h, so the survival grids are built once and the same cached interpolators are returned for every queried h.

Parameters:

  • injection_data_dir (str) – Directory containing injection CSV files matching injection_h_*_task_*.csv or injection_h_*.csv.

  • snr_threshold (float) – SNR threshold for detection. Events with SNR >= snr_threshold are considered detected.

  • h_grid (list[float] | None) – Deprecated. Previously used to specify h grid points for pre-computed grids. Now ignored (a single h-invariant survival grid serves all h). Passing this parameter emits a deprecation warning.

  • dl_bins (int) – Number of d_L grid centers for the survival grids.

  • mass_bins (int) – Number of observer-frame M_z grid centers (2D grid).

  • h_prior_range (tuple[float, float]) – Deprecated for the d_L support. Accepted for API compatibility; the d_L support is now derived from the (compact, h-invariant) detection-horizon distribution.

  • bandwidth_scale (float) – Multiplier on Scott’s-rule bandwidth for the observer-frame log-mass kernel of the 2D survival estimator.

  • estimator (Literal['local_linear', 'nadaraya_watson']) – Irrelevant to the d_L treatment (the survival function is exact in d_L). Accepted for API compatibility / the NW regression escape hatch.

  • _force_unit_weights (bool) – Internal flag for testing. When True, passes explicit weights=np.ones(N) to _build_grid_2d to verify IS estimator backward compatibility.

  • pdet_z_resolved (bool) – FIX-2 (default False = pooled, byte-identical to pre-FIX-2 behaviour). When True, every 3D (without-BH-mass) survival query returns the z-CONDITIONAL survival S(d_L | z) = P(d_hor >= d_L | z) (Gaussian kernel in u = ln(1+z), Scott d=1 bandwidth, Abramson-adaptive; exact suffix-survival in d_L per node) and the 3D accessors REQUIRE the z keyword. The 2D (M_z-conditioned) grid keeps its current form. DERIVATION_ZRESOLVED_SURVIVAL.md.

  • pdet_wbh_z_resolved (bool) – FIX-3 §7.1 (default False = pooled-in-z 2D grid, byte-identical to pre-FIX-3 behaviour). When True, the with-BH (2D) survival query returns the joint conditional S(d_L | z, M_z) (product Gaussian kernel in u = ln(1+z) and m = log10 M_z, Scott d=2 bandwidths, Abramson-adaptive on u only; exact suffix-survival in d_L; ESS-weighted (K5) shrinkage toward S(d_L | M_z)) and the 2D accessor REQUIRES the z keyword. Requires pdet_z_resolved=True (RATIFY-Z7 guard). docs/derivations/fix3_zmz_catalog_selection.md.

  • n_sky_bands (int)

  • expected_z_max (float | None)

  • allow_shallow_pool (bool)

References

Finn & Chernoff (1993), arXiv:gr-qc/9301003. Finn (1996), arXiv:gr-qc/9601048. Gray et al. (2020), arXiv:1908.06050, Section III.B-C. Mandel, Farr & Gair (2019), arXiv:1809.02063.

band_centers_sin_beta()[source]

Band centres in |sin beta| (length n_sky_bands).

Return type:

ndarray[tuple[Any, …], dtype[float64]]

band_edges_sin_beta()[source]

Equal-|sin beta| band edges (length n_sky_bands + 1) in [0, 1].

The SAME edges the inference must use to bin pixels into bands so the sky marginal is invariant (PHYSICS-CHANGE-PROTOCOL test T3).

Return type:

ndarray[tuple[Any, …], dtype[float64]]

detection_probability_with_bh_mass_interpolated(d_L, M_z, phi, theta, *, h, z=None)[source]

Detection probability including BH mass dependence (survival form).

Interpolates the 2D detection-horizon survival grid p_det(d_L, M_z) = K_M-weighted P(d_hor >= d_L) with a linear RegularGridInterpolator (bounds_error=False, fill_value=None).

Boundary handling (no extrapolation machinery — the survival grid is naturally boundary-correct):

  • d_L below the first center → clamp to the first center (survival ≈ 1 there, since the grid starts near d_L = 0).

  • d_L above the last center → 0 (no injection’s horizon reaches there).

  • M_z outside the grid range → clamped to the nearest grid edge. (fill_value=None alone would LINEARLY extrapolate — made-up but plausible-looking values; the explicit clip enforces true nearest.)

The result is monotone non-increasing in d_L and bounded in [0, 1].

Sky angles (phi, theta) are accepted for API compatibility but are marginalized over internally (D-02).

Parameters:

  • d_L (float | ndarray[tuple[Any, ...], dtype[float64]]) – Luminosity distance in Gpc.

  • M_z (float | ndarray[tuple[Any, ...], dtype[float64]]) – Observer-frame (redshifted) BH mass in solar masses.

  • phi (float | ndarray[tuple[Any, ...], dtype[float64]]) – Sky angle phi (unused, marginalized over).

  • theta (float | ndarray[tuple[Any, ...], dtype[float64]]) – Sky angle theta (unused, marginalized over).

  • h (float) – Dimensionless Hubble parameter (accepted; horizon is h-invariant).

  • z (float | ndarray[tuple[Any, ...], dtype[float64]] | None) – Conditioning redshift per query point. REQUIRED when the FIX-3 §7.1 joint estimator is active (wbh_z_resolved, atomic-switch rule); IGNORED otherwise (flag-off behaviour is byte-identical with or without z).

Returns:

Detection probability in [0, 1].

Return type:

float | ndarray[tuple[Any, …], dtype[float64]]

References

Finn & Chernoff (1993), arXiv:gr-qc/9301003; Finn (1996), arXiv:gr-qc/9601048. Gray et al. (2020), arXiv:1908.06050, Eq. (8). Mandel, Farr & Gair (2019), arXiv:1809.02063. docs/derivations/fix3_zmz_catalog_selection.md (flag-on path).

detection_probability_without_bh_mass_interpolated(d_L, phi, theta, *, h, z=None)[source]

Detection probability marginalized over BH mass (exact survival).

Drop-in replacement for DetectionProbability.detection_probability_without_bh_mass_interpolated with an additional h keyword. Returns the EXACT detection-horizon survival p_det(d_L) = P(d_hor >= d_L) (identical to the _zero_fill accessor — the survival is naturally boundary-correct, so the two accessors coincide).

# p_det = survival function of the detection horizon, P(d_hor >= d_L), # with d_hor = SNR·d_L/threshold. # Finn & Chernoff (1993), arXiv:gr-qc/9301003; Finn (1996), # arXiv:gr-qc/9601048 (p_det = P(Theta > Theta_thr)).

Sky angles (phi, theta) are accepted for API compatibility but are marginalized over internally (D-02).

Parameters:

  • d_L (float | ndarray[tuple[Any, ...], dtype[float64]]) – Luminosity distance in Gpc.

  • phi (float | ndarray[tuple[Any, ...], dtype[float64]]) – Sky angle phi (unused, marginalized over).

  • theta (float | ndarray[tuple[Any, ...], dtype[float64]]) – Sky angle theta (unused, marginalized over).

  • h (float) – Dimensionless Hubble parameter.

  • z (float | ndarray[tuple[Any, ...], dtype[float64]] | None)

Returns:

Detection probability in [0, 1].

Return type:

float | ndarray[tuple[Any, …], dtype[float64]]

References

Finn & Chernoff (1993), arXiv:gr-qc/9301003; Finn (1996), arXiv:gr-qc/9601048. Gray et al. (2020), arXiv:1908.06050, Eq. (8).

detection_probability_without_bh_mass_interpolated_zero_fill(d_L, phi, theta, *, h, z=None)[source]

Detection probability marginalized over BH mass (exact survival).

Returns the EXACT detection-horizon survival p_det(d_L) = P(d_hor >= d_L) via np.searchsorted on the stored sorted horizon with (uniform) suffix weights. This guarantees by construction:

  • p(0) = 1 (every injection’s horizon is >= 0),

  • p(d_L > max d_hor) = 0,

  • monotone non-increasing in d_L.

The function name retains the legacy _zero_fill suffix for backward compatibility with the >=6 call sites in bayesian_statistics. No bridge / slope-matched-clamp extrapolation is needed: the survival is naturally boundary-correct.

# p_det = survival function of the detection horizon, P(d_hor >= d_L), # with d_hor = SNR·d_L/threshold. # Finn & Chernoff (1993), arXiv:gr-qc/9301003; Finn (1996), # arXiv:gr-qc/9601048 (p_det = P(Theta > Theta_thr)).

Sky angles (phi, theta) are accepted for API compatibility but are marginalized over internally (D-02).

Parameters:

  • d_L (float | ndarray[tuple[Any, ...], dtype[float64]]) – Luminosity distance in Gpc.

  • phi (float | ndarray[tuple[Any, ...], dtype[float64]]) – Sky angle phi (unused, marginalized over).

  • theta (float | ndarray[tuple[Any, ...], dtype[float64]]) – Sky angle theta (unused, marginalized over).

  • h (float) – Dimensionless Hubble parameter (accepted; horizon is h-invariant).

  • z (float | ndarray[tuple[Any, ...], dtype[float64]] | None)

Returns:

Detection probability in [0, 1].

Return type:

float | ndarray[tuple[Any, …], dtype[float64]]

References

Finn & Chernoff (1993), arXiv:gr-qc/9301003; Finn (1996), arXiv:gr-qc/9601048. Gray et al. (2020), arXiv:1908.06050, Eq. (A.19).

detection_probability_without_bh_mass_sky(d_L, phi, theta, *, h, z=None)[source]

Sky-resolved detection probability p_det(d_L | Omega) (Change 1).

Maps the ecliptic sky direction (phi, theta) to the ecliptic latitude band via |sin beta| = |cos theta| (beta = pi/2 - theta) and returns that band’s detection-horizon survival, interpolated linearly in |sin beta| across band centres. phi is accepted but unused (azimuthal symmetry of the orbit-averaged response, Cutler 1998, arXiv:gr-qc/9703068).

# p_det = P(d_hor >= d_L | ecliptic-latitude band); empirical per-band # survival, azimuthally symmetric (Cutler 1998, arXiv:gr-qc/9703068).

Reduces to the pooled isotropic survival when n_sky_bands == 1 (the regression fallback; PHYSICS-CHANGE-PROTOCOL test T1).

Parameters

d_Lfloat or ndarray

Luminosity distance [Gpc].

phifloat or ndarray

Ecliptic azimuth [rad] (unused; azimuthal symmetry).

thetafloat or ndarray

Ecliptic colatitude [rad]; beta = pi/2 - theta.

hfloat

Dimensionless Hubble parameter (horizon is h-invariant).

Returns

float or ndarray

Detection probability in [0, 1].

Parameters:

  • d_L (float | ndarray[tuple[Any, ...], dtype[float64]])

  • phi (float | ndarray[tuple[Any, ...], dtype[float64]])

  • theta (float | ndarray[tuple[Any, ...], dtype[float64]])

  • h (float)

  • z (float | ndarray[tuple[Any, ...], dtype[float64]] | None)

Return type:

float | ndarray[tuple[Any, …], dtype[float64]]

get_dl_max(h)[source]

Return the maximum d_L of the 1D P_det grid for the given h value.

This is the upper edge of the d_L support, i.e. the maximum detection horizon padded by _DL_PADDING_FACTOR. Needed to compute z_max(h) for the full-volume denominator integral.

Parameters:

h (float) – Dimensionless Hubble parameter (horizon is h-invariant).

Returns:

Maximum d_L in Gpc.

Return type:

float

quality_flags(h)[source]

Return per-bin quality metadata for the given h value.

Quality flags are diagnostic metadata. They do not affect the P_det survival result. If the grid for this h has not been built yet, it will be built (the single h-invariant survival grid).

The returned dict contains:

  • n_total: float array (dl_bins, M_bins) – injection count per cell.

  • n_detected: float array (dl_bins, M_bins) – detected count per cell (d_hor >= cell lower dl edge).

  • reliable: bool array (dl_bins, M_bins) – True where n_total >= 10.

  • dl_edges: float array (dl_bins+1,) – d_L bin edges in Gpc.

  • M_edges: float array (M_bins+1,) – M_z bin edges.

  • n_eff: float array (dl_bins, M_bins) – effective sample size (= n_total under uniform weights).

Parameters:

h (float) – Hubble parameter value.

Returns:

Dict of quality flag arrays.

Raises:

ValueError – If no quality flags are available after construction.

Return type:

dict[str, ndarray[tuple[Any, …], dtype[float64]] | ndarray[tuple[Any, …], dtype[bool]]]

survival_per_band(d_L, z=None)[source]

Per-band detection-horizon survival S_b(d_L); shape (n_sky_bands, Nq).

The building block of the sky-resolved selection integrals: the inference forms the sky sum (1/Npix) sum_k p_det(d_L, Omega_k) as sum_b (n_pix_b/Npix) S_b(d_L) (each pixel takes its band’s flat survival), and the missing-completion integral weights S_b by the per-band incompleteness (1/Npix) sum_{k in b}(1 - f_k(z)).

When the FIX-2 z-resolved estimator is active (pdet_z_resolved), z is REQUIRED and each band returns the z-conditional S(d_L | z, band b) (with the per-cell ESS-floor fallback to the z-only conditional, packet §7).

Parameters

d_Lfloat or ndarray

Luminosity distance query points [Gpc].

zfloat or ndarray, optional

Conditioning redshift per query point (FIX-2 only; required when z_resolved is True, ignored otherwise).

Returns

ndarray, shape (n_sky_bands, Nq)

Band-resolved survival in [0, 1].

Parameters:

  • d_L (float | ndarray[tuple[Any, ...], dtype[float64]])

  • z (float | ndarray[tuple[Any, ...], dtype[float64]] | None)

Return type:

ndarray[tuple[Any, …], dtype[float64]]

validate_coverage(h, crb_df)[source]

Compute fraction of events whose 4-sigma d_L bounds fall within the P_det grid.

For each event, compute d_L +/- 4*sigma_dL from the Cramer-Rao bounds. Check if both bounds fall within the grid’s d_L range.

Parameters:

  • h (float) – Hubble parameter value (to build/retrieve grid).

  • crb_df (DataFrame) – DataFrame with columns luminosity_distance and delta_luminosity_distance_delta_luminosity_distance (variance).

Returns:

Coverage fraction in [0, 1].

Return type:

float

wbh_joint_knot_values(d_L, z)[source]

M_z knots and shrunk-survival values for the erf-sum path (§3.3-C).

Returns the m-node knots lifted to observer-frame mass M_z,j = 10^{m_j} together with S̃ evaluated at (d_L_i, u(z_i), m_j) for every knot j — two-u-node linear blend plus linear interpolation along DLQ, NO m-interpolation (the values ARE the knots). The erf-sum consumer treats the interpolant as PIECEWISE-LINEAR IN M_z between these lifted knots ([RATIFY-Z3] §3.3-C convention 2 choice (a)), keeping its closed form exact (fix3_zmz_catalog_selection.md §3.5 erf-sum correction).

Parameters:

  • d_L (ndarray[tuple[Any, ...], dtype[float64]]) – query luminosity distances [Gpc], shape (n,).

  • z (ndarray[tuple[Any, ...], dtype[float64]]) – conditioning redshifts, shape (n,) (broadcastable to d_L).

Returns:

(M_z_knots, S_values) with shapes (n_m,) and (n, n_m).

Raises:

ValueError – if the joint estimator is not active.

Return type:

tuple[ndarray[tuple[Any, …], dtype[float64]], ndarray[tuple[Any, …], dtype[float64]]]

property wbh_z_resolved: bool

True iff the joint z x M_z-resolved with-BH (FIX-3 §7.1) estimator is active.

property z_resolved: bool

True iff the z-resolved (FIX-2) survival estimator is active.

Posterior Combination

Posterior combination module for merging per-event likelihoods.

Combines per-h-value posterior JSON files from the evaluation campaign into a joint posterior over the Hubble constant, with multiple zero-handling strategies and diagnostic reporting.

Four zero-handling strategies: - naive: Replace 0.0 with np.finfo(float).tiny before log. - exclude: Remove events that have any zero likelihood. - per-event-floor: Replace 0.0 with min(nonzero) / 100 per event. - physics-floor: Per-event minimum nonzero likelihood as floor value; all-zero events excluded.

class master_thesis_code.bayesian_inference.posterior_combination.CombinationStrategy(*values)[source]

Bases: StrEnum

Zero-handling strategy for posterior combination.

EXCLUDE = 'exclude'
NAIVE = 'naive'
PER_EVENT_FLOOR = 'per-event-floor'
PHYSICS_FLOOR = 'physics-floor'
master_thesis_code.bayesian_inference.posterior_combination.apply_strategy(likelihoods, strategy)[source]

Apply a zero-handling strategy to the likelihood array.

Parameters

likelihoodsndarray of shape (n_events, n_h_values)

Raw likelihood array (may contain 0.0 entries).

strategyCombinationStrategy

Which strategy to apply.

Returns

processedndarray

Likelihood array with zeros handled.

excluded_countint

Number of events removed (0 for strategies that keep all events).

Parameters:
Return type:

tuple[ndarray[tuple[Any, …], dtype[float64]], int]

master_thesis_code.bayesian_inference.posterior_combination.build_likelihood_array(h_values, event_likelihoods)[source]

Build a 2-D likelihood array from the nested dict structure.

Parameters

h_valueslist[float]

Sorted h-values (columns).

event_likelihoodsdict[int, dict[float, float]]

Nested dict from load_posterior_jsons().

Returns

arrayndarray of shape (n_events, n_h_values)

Likelihood values. NaN where an event is missing a particular h-value; 0.0 where the likelihood was explicitly zero.

detection_indiceslist[int]

Sorted detection indices (row labels).

Parameters:

  • h_values (list[float])

  • event_likelihoods (dict[int, dict[float, float]])

Return type:

tuple[ndarray[tuple[Any, …], dtype[float64]], list[int]]

master_thesis_code.bayesian_inference.posterior_combination.combine_log_space(likelihoods, log_D_h=None, n_events_used=0)[source]

Combine per-event likelihoods into a joint posterior using log-space.

The per-event likelihoods \(L_i(h)\) are already conditional on detection: L_comp = num/D(h) (Gray et al. 2020 Eq. 31, where D(h) plays the role of the prior normalization for p_galaxy ∝ p_det · dV_c) and L_cat = (1/N_g) Σ_g (N_g_local/D_g_local) (per-galaxy detection-conditional expectation). The joint posterior is therefore just Π_i L_i(h) with no additional β(h)^N selection correction (Loredo 2004; Mandel et al. 2019, arXiv:1809.02063 §3): conditioning on the observed N_detected makes the β^N factor in the unconditional likelihood cancel against the Poisson rate prior.

The log_D_h and n_events_used arguments are retained for backward compatibility with older callers, but are now ignored — applying −N · log D(h) here would double-count the selection correction already inside each per-event L_comp.

Parameters

likelihoodsndarray of shape (n_events, n_h_values)

Likelihood array with zeros already handled (all values > 0).

log_D_hndarray of shape (n_h_values,) or None

Ignored. Kept for backward compatibility; pass None for new callers. Phase 43-H1 commit 2853c32 introduced an −N · log D(h) subtraction here that we have since shown (Tier 3 audit, 2026-05-04) over-applies the selection correction.

n_events_usedint

Ignored. Kept for backward compatibility.

Returns

posteriorndarray of shape (n_h_values,)

Normalized joint posterior.

Parameters:

  • likelihoods (ndarray[tuple[Any, ...], dtype[float64]])

  • log_D_h (ndarray[tuple[Any, ...], dtype[float64]] | None)

  • n_events_used (int)

Return type:

ndarray[tuple[Any, …], dtype[float64]]

master_thesis_code.bayesian_inference.posterior_combination.combine_posteriors(posteriors_dir, strategy, output_dir, d_h_table=None, allow_shallow_pool=False)[source]

Combine per-event posteriors into a joint posterior.

This is the main entry point called from CLI.

Parameters

posteriors_dirstr

Path to directory containing h_*.json files.

strategystr

Zero-handling strategy name (one of the CombinationStrategy values).

output_dirstr

Path to write output files.

d_h_tabledict[float, float] or None

Pre-computed {h: D(h)} mapping (Gray et al. 2020, Eq. A.19). When None, D(h) is computed automatically using precompute_completion_denominator().

Returns

resultdict

Combined posterior result with keys h_values, posterior, strategy, n_events_total, n_events_used, n_events_excluded, n_events_empty, map_h, map_posterior, variant, D_h_per_h.

Parameters:

  • posteriors_dir (str)

  • strategy (str)

  • output_dir (str)

  • d_h_table (dict[float, float] | None)

  • allow_shallow_pool (bool)

Return type:

dict[str, object]

master_thesis_code.bayesian_inference.posterior_combination.compute_canonical_combined_posterior(posteriors_dir)[source]

Compute the canonical joint H0 posterior from per-h JSONs.

Aggregates per-event log-likelihoods via raw Σ log L_i (no physics-floor, no outer D(h) correction — the Tier 3 fix removed the latter from combine_log_space). This is the combination quoted in Phase 48 verdict JSON and throughout docs/H0_BIAS_RESOLUTION.md.

Parameters

posteriors_dir:

Directory of h_*.json files (1D or 2D variant).

Returns

dict with keys

  • h_values: list[float] — sorted h-grid

  • posterior: list[float] — peak-normalised posterior

  • log_posterior: list[float] — un-normalised Σ log L_i

  • n_events_used: int — events surviving the full-mask filter

  • discrete_map: float — argmax on the discrete grid

  • continuous_map: float — parabolic-refined sub-grid MAP

  • strategy: literal string "raw-sum-log"

Parameters:

posteriors_dir (Path)

Return type:

dict[str, object]

master_thesis_code.bayesian_inference.posterior_combination.generate_comparison_table(h_values, likelihoods, detection_indices, variant, log_D_h=None)[source]

Generate a markdown comparison table for all strategies.

Parameters

h_valuesndarray

Sorted h-values.

likelihoodsndarray of shape (n_events, n_h_values)

Raw likelihood array.

detection_indiceslist[int]

Detection indices.

variantstr

Name of the posterior variant (e.g. "posteriors").

log_D_hndarray of shape (n_h_values,) or None

Ignored (kept for backward compatibility). See combine_log_space docstring — D(h) enters via L_comp = num/D inside each per-event likelihood (Gray Eq. 31), not as an outer correction.

Returns

tablestr

Markdown-formatted comparison table.

Parameters:

  • h_values (ndarray[tuple[Any, ...], dtype[float64]])

  • likelihoods (ndarray[tuple[Any, ...], dtype[float64]])

  • detection_indices (list[int])

  • variant (str)

  • log_D_h (ndarray[tuple[Any, ...], dtype[float64]] | None)

Return type:

str

master_thesis_code.bayesian_inference.posterior_combination.generate_diagnostic_report(h_values, likelihoods, detection_indices)[source]

Generate a markdown diagnostic report about zero-likelihood events.

Parameters

h_valueslist[float]

Sorted h-values.

likelihoodsndarray of shape (n_events, n_h_values)

Raw likelihood array (before strategy application).

detection_indiceslist[int]

Detection indices corresponding to rows.

Returns

reportstr

Markdown-formatted diagnostic report.

Parameters:

  • h_values (list[float])

  • likelihoods (ndarray[tuple[Any, ...], dtype[float64]])

  • detection_indices (list[int])

Return type:

str

master_thesis_code.bayesian_inference.posterior_combination.load_per_h_likelihoods(directory)[source]

Load per-event log-likelihoods across all h-value JSONs in directory.

Reads every h_*.json file, extracts the scalar per-event likelihood for each integer event key, drops events that have missing data at any h-value (full_mask filter), and returns (h_values, log_L) where log_L has shape (n_events_full, n_h_values).

This is the reference implementation used by the bias-investigation suite and is the canonical loader for all paper-grade H0 posterior figures.

Parameters

directory:

Path to a directory containing h_*.json files (typically posteriors/ or posteriors_with_bh_mass/).

Returns

h_values:

Sorted list of h-values (one per JSON file).

log_L:

Float64 array of shape (n_events, n_h_values); entries are log(max(L_event_at_h, 1e-300)). Events with NaN at any h are excluded.

Parameters:

directory (Path)

Return type:

tuple[list[float], ndarray[tuple[Any, …], dtype[float64]]]

master_thesis_code.bayesian_inference.posterior_combination.load_posterior_jsons(posteriors_dir)[source]

Load per-h-value posterior JSON files from a directory.

Parameters

posteriors_dirPath

Directory containing h_*.json files.

Returns

h_valueslist[float]

Sorted list of h values found across all files.

event_likelihoodsdict[int, dict[float, float]]

Nested dict mapping detection_index -> {h_value: likelihood}. Events with empty lists (missing data) are excluded.

Parameters:

posteriors_dir (Path)

Return type:

tuple[list[float], dict[int, dict[float, float]]]

master_thesis_code.bayesian_inference.posterior_combination.parabolic_refine_map(h_grid, log_posterior)[source]

Refine the discrete MAP via 3-point parabolic interpolation.

Given a log-posterior on a discrete h-grid, fits a parabola through the discrete argmax and its two neighbours and returns the analytic vertex of the parabola. Falls back to the discrete argmax when the peak is on the grid boundary or the parabola degenerates.

Parameters

h_grid:

Monotonically increasing h-values.

log_posterior:

Log-posterior values on h_grid.

Returns

float

Sub-grid continuous MAP estimate.

Parameters:

  • h_grid (ndarray[tuple[Any, ...], dtype[float64]])

  • log_posterior (ndarray[tuple[Any, ...], dtype[float64]])

Return type:

float

Evaluation Report

Evaluation report module for baseline extraction and comparison.

Provides tools to extract a baseline H0 posterior snapshot from an h-sweep, compute credible intervals, and generate before/after comparison reports.

Used by Phases 31-34 to measure the effect of each fix on the H0 posterior.

class master_thesis_code.bayesian_inference.evaluation_report.BaselineSnapshot(map_h, ci_lower, ci_upper, ci_width, bias_percent, n_events, h_values=<factory>, log_posteriors=<factory>, per_event_summaries=<factory>, n_excluded_fisher=0, median_cond_3d=0.0, median_cond_4d=0.0, created_at=<factory>, git_commit=<factory>)[source]

Bases: object

Snapshot of H0 posterior metrics extracted from a posteriors directory.

Parameters:

  • map_h (float) – MAP (maximum a posteriori) Hubble constant value.

  • ci_lower (float) – Lower bound of the 68% credible interval.

  • ci_upper (float) – Upper bound of the 68% credible interval.

  • ci_width (float) – Width of the 68% credible interval (ci_upper - ci_lower).

  • bias_percent (float) – Relative bias as (MAP - true_h) / true_h * 100.

  • n_events (int) – Number of detection events contributing to the posterior.

  • h_values (list[float]) – Sorted list of h values in the sweep.

  • log_posteriors (list[float]) – Log-posterior values at each h.

  • per_event_summaries (list[dict[str, float]]) – Per-event diagnostic data (d_L, SNR, etc.).

  • created_at (str) – ISO 8601 timestamp of when this snapshot was created.

  • git_commit (str) – Git commit hash at time of creation.

  • n_excluded_fisher (int)

  • median_cond_3d (float)

  • median_cond_4d (float)

bias_percent: float
ci_lower: float
ci_upper: float
ci_width: float
created_at: str
classmethod from_json(data)[source]

Deserialize from a JSON-compatible dictionary.

Parameters:

data (dict[str, object]) – Dictionary as produced by to_json().

Returns:

A new BaselineSnapshot instance.

Return type:

BaselineSnapshot

git_commit: str
h_values: list[float]
log_posteriors: list[float]
map_h: float
median_cond_3d: float = 0.0
median_cond_4d: float = 0.0
n_events: int
n_excluded_fisher: int = 0
per_event_summaries: list[dict[str, float]]
to_json()[source]

Serialize to a JSON-compatible dictionary.

Returns:

Dictionary representation suitable for json.dumps.

Return type:

dict[str, object]

master_thesis_code.bayesian_inference.evaluation_report.compute_credible_interval(h_values, log_posteriors, level=0.68)[source]

Compute a credible interval from log-posterior values.

Converts log-posteriors to a normalized probability distribution, then computes the CDF and finds where it crosses (1-level)/2 and (1+level)/2.

Parameters:

  • h_values (list[float]) – Sorted list of h values.

  • log_posteriors (list[float]) – Log-posterior value at each h.

  • level (float) – Credible interval level (default 0.68 for 68% CI).

Returns:

Tuple (lower, upper) bounding the credible interval.

Return type:

tuple[float, float]

master_thesis_code.bayesian_inference.evaluation_report.extract_baseline(posteriors_dir, crb_csv_path=None, true_h=0.73)[source]

Extract baseline H0 posterior metrics from a posteriors directory.

Warning

This function does NOT apply the -N log D(h) selection-effect normalization from Gray et al. (2020) arXiv:1908.06050 Eq. A.19. The MAP returned here is systematically biased toward h_max in standard ΛCDM cosmology (where D(h) grows with h), because the log-posterior Σ_i log L_i(h) is a monotone-increasing function of h without the D(h) correction.

Use --evaluate (BayesianStatistics.evaluate) for production MAP reporting. extract_baseline is diagnostic-only — suitable for comparing relative shifts between h-sweeps run under identical conditions, but not for absolute MAP recovery.

Parameters:

  • posteriors_dir (Path) – Directory containing h_*.json files from an h-sweep.

  • crb_csv_path (Path | None) – Optional path to CRB CSV for per-event summaries.

  • true_h (float) – True Hubble constant value for bias computation. Default 0.73.

Returns:

BaselineSnapshot with MAP h, 68% CI, bias %, and event count.

Raises:

ValueError – If fewer than 3 h-value files are found (insufficient for credible interval computation per D-02).

Return type:

BaselineSnapshot

master_thesis_code.bayesian_inference.evaluation_report.generate_comparison_report(baseline, current, output_dir, label='current')[source]

Generate a comparison report between a baseline and current posterior snapshot.

Writes both a human-readable Markdown report and a machine-readable JSON sidecar.

Parameters:

  • baseline (BaselineSnapshot) – The reference (pre-change) BaselineSnapshot.

  • current (BaselineSnapshot) – The current (post-change) BaselineSnapshot.

  • output_dir (Path) – Directory to write output files into (created if needed).

  • label (str) – Label suffix for output file names (default “current”).

Returns:

Path to the generated Markdown report.

Return type:

Path

master_thesis_code.bayesian_inference.evaluation_report.generate_diagnostic_summary(diagnostic_csv_path, output_dir, label='diagnostic')[source]

Analyze per-event diagnostic CSV and generate explanatory summary.

Reads the diagnostic CSV produced by BayesianStatistics.evaluate() and computes statistics explaining WHY the posterior bias changes: - Mean/median f_i across events (catalog completeness fraction) - L_comp contribution statistics - Fraction of events where L_comp pulls toward lower h

The “L_comp pulls toward lower h” metric compares L_comp at h=0.66 vs h=0.73 per event. If L_comp(h=0.66) > L_comp(h=0.73), the completion term biases that event toward lower h.

Parameters:

  • diagnostic_csv_path (Path) – Path to event_likelihoods.csv from evaluate().

  • output_dir (Path) – Directory to write the summary report.

  • label (str) – Label suffix for output file name.

Returns:

Dict with summary statistics (also written to JSON).

Return type:

dict[str, object]

master_thesis_code.bayesian_inference.evaluation_report.load_posteriors(posteriors_dir)[source]

Load h_*.json posterior files from a directory.

Each file should contain a dict with key “h” (float) and integer-string keys mapping to [likelihood_value] arrays for each detection event.

Parameters:

posteriors_dir (Path) – Path to directory containing h_*.json files.

Returns:

“h”, “log_posterior”, “n_detections”. Sorted by h value in ascending order.

Return type:

List of dicts with keys