HOD Fitting Module

The hod_mod.fitting module provides two-stage HOD fitting of observed galaxy clustering data: a fast MAP estimate via Nelder-Mead optimisation, followed by an MCMC posterior exploration with emcee.

—

Statistical framework

Given observed \(w_p\) values \(\mathbf{d}\) with covariance matrix \(\mathbf{C}\), the log-likelihood under a Gaussian noise model is

\[\ln\mathcal{L}(\boldsymbol{\theta}) = -\frac{1}{2}\left(\mathbf{d} - \mathbf{m}(\boldsymbol{\theta})\right)^T \mathbf{C}^{-1} \left(\mathbf{d} - \mathbf{m}(\boldsymbol{\theta})\right)\]

where \(\mathbf{m}(\boldsymbol{\theta})\) is the model prediction for parameter vector \(\boldsymbol{\theta}\). In terms of a reduced chi-squared:

\[\chi^2_\nu = \frac{1}{\nu} \left(\mathbf{d} - \mathbf{m}\right)^T \mathbf{C}^{-1} \left(\mathbf{d} - \mathbf{m}\right)\]

where \(\nu = N_{\rm data} - N_{\rm free}\) is the number of degrees of freedom.

The log-posterior is

\[\ln P(\boldsymbol{\theta} | \mathbf{d}) = \ln\mathcal{L}(\boldsymbol{\theta}) + \ln\pi(\boldsymbol{\theta})\]

where \(\ln\pi(\boldsymbol{\theta})\) is the log-prior (see Prior distributions below).

Covariance matrix

When data are loaded from a sum_stat HDF5 file, the full covariance matrix from the file is used. A 1% diagonal regularisation is applied before inversion to guard against numerical singularities:

\[\mathbf{C}_{\rm reg} = \mathbf{C} + 0.01 \cdot \mathrm{diag}(\mathbf{C})\]

When data are loaded from a CSV file (legacy), only diagonal errors are available and \(C_{ij} = \sigma_i^2 \delta_{ij}\).

Prior distributions

Uniform prior (default for HOD parameters)

For each free parameter \(\theta_i\) with bounds \([l_i, u_i]\):

\[\begin{split}\ln\pi(\theta_i) = \begin{cases} 0 & l_i \le \theta_i \le u_i \\ -\infty & \text{otherwise}\end{cases}\end{split}\]

Gaussian prior (optional, for cosmological parameters)

For each parameter with Gaussian prior \(\mathcal{N}(\mu_i, \sigma_i^2)\):

\[\begin{split}\ln\pi(\theta_i) = \begin{cases} -\dfrac{(\theta_i - \mu_i)^2}{2\sigma_i^2} & l_i \le \theta_i \le u_i \\ -\infty & \text{otherwise} \end{cases}\end{split}\]

where hard bounds \([l_i, u_i]\) are applied in addition. When param_prior_types[name] = "gaussian" the fitter adds this term on top of the chi-squared.

—

Planck 2018 cosmological prior

(hod_mod.fitting.planck_prior)

The Planck 2018 TT,TE,EE+lowE best-fit values and 1σ uncertainties (Planck Collaboration 2020, Table 2) [PlanckCollaboration2018] are encoded in PLANCK18_MEANS and PLANCK18_SIGMAS:

Parameter	Symbol	Best-fit	1σ	3σ range
Hubble constant	\(h\)	0.6736	0.0054	[0.6574, 0.6898]
Matter density	\(\Omega_m\)	0.3153	0.0073	[0.2934, 0.3372]
Baryon density	\(\Omega_b\)	0.0493	0.0006	[0.0475, 0.0511]
Spectral index	\(n_s\)	0.9649	0.0042	[0.9523, 0.9775]
Log amplitude	\(\ln 10^{10}A_s\)	3.044	0.014	[3.002, 3.086]
\(\sigma_8\)	\(\sigma_8\)	0.8111	0.0060	[0.7931, 0.8291]

The 3σ bounds are used as hard truncation limits for the Gaussian priors. The function planck18_log_prior(theta) returns the sum of all Gaussian log-prior terms; it returns \(-\infty\) if any parameter is outside its 3σ range.

Usage in fitting scripts:

from hod_mod.fitting.planck_prior import PLANCK18_MEANS, PLANCK18_SIGMAS
from hod_mod.fitting import WpFitConfig, WpFitter

cfg = WpFitConfig(
    ...
    param_prior_types  = {"h": "gaussian", "Omega_m": "gaussian"},
    param_prior_means  = {"h": PLANCK18_MEANS["h"],
                          "Omega_m": PLANCK18_MEANS["Omega_m"]},
    param_prior_sigmas = {"h": PLANCK18_SIGMAS["h"],
                          "Omega_m": PLANCK18_SIGMAS["Omega_m"]},
)

—

MAP estimation

WpFitter.map_fit() minimises \(-\ln P\) using the Nelder-Mead simplex algorithm (via scipy.optimize.minimize). Nelder-Mead is gradient-free and robust to the discontinuous derivatives that arise from hard prior bounds. The result is the MAP (maximum a posteriori) point estimate.

A good starting point is critical: param_init in the config should be set to a physically plausible value (e.g. set log10mmin close to the expected characteristic halo mass for the sample).

—

MCMC posterior sampling

WpFitter.mcmc_fit() uses the emcee ensemble sampler (Foreman-Mackey et al. 2013) [Foreman-Mackey2013]. The default configuration uses 64 walkers initialised in a Gaussian ball around the MAP estimate.

Convergence diagnostics:

• Acceptance fraction: should be between 0.2 and 0.5. If too low, the proposal scale moves parameter needs tuning.

• Integrated autocorrelation time \(\hat{\tau}\): the chain is considered converged when the number of steps exceeds \(50\hat{\tau}\). Access via sampler.get_autocorr_time().

• Gelman-Rubin statistic \(\hat{R}\): for multi-chain runs, \(\hat{R} < 1.1\) indicates convergence.

—

Configuration

Fitting is driven by a WpFitConfig dataclass (or equivalently a YAML file parsed by load_config):

data_file:    /path/to/data.h5
data_format:  hdf5          # "csv" or "hdf5"
rp_min:       0.3           # Mpc/h
rp_max:       30.0          # Mpc/h
hod_model:    MoreHODModel
hmf_backend:  csst           # pipeline baseline (default if omitted); use
                              # tinker08 to reproduce literature results
z:            0.15
pi_max:       100.0         # Mpc/h

free_params:    [log10mmin, sigma_logm, log10m1, alpha]
param_bounds:
  log10mmin:  [11.0, 13.5]
  sigma_logm: [0.1,  1.0]
  log10m1:    [12.0, 15.0]
  alpha:      [0.5,  2.0]
param_init:
  log10mmin:  12.5
  sigma_logm: 0.38
  log10m1:    13.5
  alpha:      1.0
  kappa:      1.0
  alpha_inc:  1.0
  log10m_inc: 12.0

# Optional: Gaussian cosmological priors
param_prior_types:  {h: gaussian}
param_prior_means:  {h: 0.6736}
param_prior_sigmas: {h: 0.0054}

output_dir: results/bgs_ls10/mstar10.5/

—

Usage example

from hod_mod.fitting import WpFitter, load_config

cfg    = load_config("configs/hod_fit_more2015_cmass.yml")
fitter = WpFitter(cfg)
result = fitter.map_fit()           # Nelder-Mead MAP estimate
print(result.x, result.fun)

sampler = fitter.mcmc_fit()         # emcee MCMC posterior
flat    = sampler.get_chain(flat=True, discard=100, thin=5)

HOD fitting: MAP + MCMC via emcee, unified FitConfig, Planck 2018 prior.

Planck 2018 cosmological priors for HOD fitting.

Provides best-fit values, 1σ uncertainties, and 3σ flat bounds from the Planck 2018 primary CMB analysis (TT,TE,EE+lowE likelihood, Table 2).

Reference

Planck Collaboration 2020, A&A 641, A6 https://arxiv.org/abs/1807.06209

The primary parameters and their 68% confidence intervals are:

\[\begin{split}h &= 0.6736 \pm 0.0054 \\ \Omega_m &= 0.3153 \pm 0.0073 \\ \Omega_b h^2 &= 0.02237 \pm 0.00015 \\ n_s &= 0.9649 \pm 0.0042 \\ \ln 10^{10} A_s &= 3.044 \pm 0.014\end{split}\]

The 3σ flat bounds are \([\mu - 3\sigma, \mu + 3\sigma]\).

Usage in YAML config

Set prior_type: gaussian for any cosmological parameter to activate the Gaussian prior. The bounds field is still required and acts as hard clipping beyond which the log-prior returns -inf:

parameters:
  h:
    free: true
    init: 0.6736
    bounds: [0.6574, 0.6898]   # 3σ hard bounds
    prior_type: gaussian
    prior_mean: 0.6736
    prior_sigma: 0.0054

hod_mod.fitting.planck_prior.gaussian_log_prior(val: float, mean: float, sigma: float, lo: float = -inf, hi: float = inf) → float

Gaussian log-prior for a single parameter.

\[\ln \pi(\theta) = -\frac{1}{2} \left( \frac{\theta - \mu}{\sigma} \right)^2\]

Returns -inf if val is outside [lo, hi] (hard bounds).

Parameters:

• val (float)

• mean, sigma (float) – Gaussian mean and standard deviation.

• lo, hi (float) – Hard bounds (uniform outside returns -inf).

hod_mod.fitting.planck_prior.planck18_log_prior(theta: dict, params: list | None = None) → float

Sum of Gaussian log-prior terms for Planck 2018 cosmological parameters.

\[\ln \pi(\theta) = -\frac{1}{2} \sum_i \left( \frac{\theta_i - \mu_i}{\sigma_i} \right)^2\]

Parameters:

• theta (dict) – Parameter dict. Only keys present in PLANCK18_MEANS contribute.

• params (list of str, optional) – Restrict to these parameters only. Default: all keys in PLANCK18_MEANS that also appear in theta.

Returns:

float – Log-prior value. Returns -inf if any parameter is outside its 3σ hard bounds.

—

Key references

Lensing and galaxy–matter cross-correlation: [BartelmannSchneider2001], [Mandelbaum2005], [Mandelbaum2006], [Leauthaud2017], [Miyatake2022], [Lange2023], [Heydenreich2025], [Lange2025].

Intrinsic alignments: [Catelan2001], [HirataSeljak2004], [Brown2002], [BridleKing2007], [Blazek2019], [DESI_KP6].

Surveys: [Blanton2003], [BOSS_CMASS], [HSC_Aihara2018], [HSC_Mandelbaum2018], [KiDS_Heymans2021], [DES_Abbott2022], [DESI_EDR], [DESI_BGS_Hahn2023], [Comparat2023], [Lange2024], [Lange2025phz].

Inference: [Phan2019].