Fitting with the More+2015 HOD Model

Overview

hod_mod/scripts/fitting/run_fit_More15.py (also hod-mod fit) is a command-line tool for fitting the FullHaloModelPrediction forward model to observational data using the More et al. (2015) Halo Occupation Distribution (HOD) model.

The script is a thin wrapper around the validated fitting classes in hod_mod.fitting:

  • WpFitter — projected clustering \(w_p(r_p)\) only

  • DeltaSigmaFitter — excess surface density \(\Delta\Sigma(R)\) only

  • JointFitter — joint \(w_p + \Delta\Sigma + n_\mathrm{gal}\)

The probe mode is detected automatically from the YAML configuration (see Probe Modes).

Quick Start

Step 1 — Write a configuration file (see examples in configs/fitting/):

label: "My galaxy sample"
data:
  file:   path/to/wp_data.csv
  format: bwpd
  rp_min: 0.5
  rp_max: 50.0
cosmology:
  Omega_m: 0.310
  h:       0.700
  sigma8:  0.800
  n_s:     0.965
  Omega_b: 0.045
model:
  hod_model:   MoreHODModel
  hmf_backend: tinker08
  z:           0.50
  pi_max:      100.0
parameters:
  log10mmin:  {free: true,  bounds: [11.0, 15.0], init: 12.5}
  sigma_logm: {free: true,  bounds: [0.05, 1.50], init: 0.4}
  log10m1:    {free: true,  bounds: [12.0, 15.5], init: 13.5}
  alpha:      {free: true,  bounds: [0.50, 2.50], init: 1.0}
  kappa:      {free: true,  bounds: [0.01, 3.00], init: 1.0}
fitting:
  method: map
output:
  dir: results/my_fit/

Step 2 — Run:

python hod_mod/scripts/fitting/run_fit_More15.py  config.yml --map-only

Step 3 — Check outputs in results/my_fit/:

fit_result.json      # best-fit params, χ²/dof
fit_wp.png           # wp data vs best-fit
fit_combined.png     # rp·wp panel
fit_hod.png          # N_c(M), N_s(M) occupation

YAML Configuration Reference

data: block (wp)

Required when fitting \(w_p(r_p)\).

Field

Type

Default

Description

file

string

—

Path to \(w_p\) data file (CSV, HDF5, or FITS).

format

string

—

Data format. bwpd (More+2015-style three-column), hdf5, fits.

rp_min

float

0.0

Minimum \(r_p\) (\(h^{-1}\) Mpc) to include in fit.

rp_max

float

∞

Maximum \(r_p\) (\(h^{-1}\) Mpc) to include in fit.

joint: block (ESD + n_gal)

Include to activate joint \(w_p + \Delta\Sigma + n_\mathrm{gal}\) fitting. Omitting this block makes the fit wp-only.

Field

Type

Default

Description

ds_file

string

—

Path to \(\Delta\Sigma\) data file.

ds_format

string

—

ESD format. bwpd_4col (R, DS, upper, lower), hdf5.

ds_rp_min

float

0.0

Minimum \(R\) (\(h^{-1}\) Mpc) for ESD.

ds_rp_max

float

∞

Maximum \(R\) (\(h^{-1}\) Mpc) for ESD.

ng_obs

float

—

Observed galaxy number density \(h^3\) Mpc:math:^{-3}.

ng_frac_err

float

0.20

Fractional error on \(n_\mathrm{gal}\) (Gaussian term in likelihood).

ds: block (ESD-only mode)

Use ds: instead of joint: to fit \(\Delta\Sigma\) only (without \(w_p\)). The fields are identical to the joint: block above.

cosmology: block

Field

Type

Required

Description

Omega_m

float

yes

Total matter density parameter at z=0.

h

float

yes

Dimensionless Hubble constant (\(H_0 = 100\,h\) km/s/Mpc).

sigma8

float

yes

RMS matter fluctuations at 8 \(h^{-1}\) Mpc.

n_s

float

yes

Primordial spectral index.

Omega_b

float

yes

Baryon density parameter.

model: block

Field

Type

Default

Description

hod_model

string

MoreHODModel

HOD model class. Use MoreHODModel for More+2015.

hmf_backend

string

tinker08

Halo mass function backend. tinker08 recommended.

z

float

—

Galaxy sample effective redshift.

pi_max

float

100.0

Line-of-sight integration limit \(\pi_\mathrm{max}\) (Mpc/\(h\)).

use_bnl

bool

false

Enable beyond-linear halo bias (Mead & Verde 2021). Negligible at \(r_p > 1\,h^{-1}\) Mpc.

parameters: block

Each entry: name: {free: bool, bounds: [lo, hi], init: value}.

If free: false, only init is required (the parameter is held fixed).

fitting: block

Field

Type

Default

Description

method

string

map

map — MAP only; mcmc — MCMC only (starts at init); both — MAP then MCMC.

n_walkers

int

32

Number of emcee walkers.

n_steps

int

2000

Number of MCMC steps per walker.

n_burnin

int

500

Burn-in steps discarded before saving the flat chain.

output: block

Field

Type

Default

Description

dir

string

results/fit/

Directory for all output files.

label: (top-level)

Optional string identifying this fit in figure titles and JSON output. Defaults to the YAML filename stem if absent.

published_params: (top-level, optional)

Dictionary mapping parameter names to [best-fit, error] pairs. When present, a dashed reference curve is overlaid on all data comparison plots, and MCMC corner plots mark the reference values as vertical/horizontal lines.

Example:

published_params:
  log10mmin:  [13.13, 0.13]
  sigma_logm: [0.469, 0.13]
  log10m1:    [14.21, 0.12]
  alpha:      [1.130, 0.09]
  kappa:      [1.250, 0.40]

Probe Modes

The fitter class is selected automatically based on which data sections appear in the YAML:

Probes in config

Fitter class

Likelihood

data: only

WpFitter

\(\chi^2_{w_p}\)

ds: only (no data:)

DeltaSigmaFitter

\(\chi^2_{\Delta\Sigma}\)

data: + joint:

JointFitter

\(\chi^2_{w_p} + \chi^2_{\Delta\Sigma} + \chi^2_{n_\mathrm{gal}}\)

Covariance Options

Diagonal covariance (CSV / bwpd format)

When data.file is a CSV with bwpd or bwpd_4col format, the covariance is diagonal: \(C_{ii} = \sigma_i^2\) from the error column. This is appropriate when data points are independently measured or when off-diagonal correlations are negligible.

Full covariance (HDF5/FITS jackknife)

When the data file is an HDF5/FITS jackknife catalogue, load_config() builds the full jackknife covariance matrix and applies the Hartlap et al. (2007) correction:

\[\hat{C}^{-1} = \frac{n_\mathrm{jk} - n_\mathrm{bins} - 2}{n_\mathrm{jk} - 1}\, C^{-1}_\mathrm{JK}\]

To use jackknife data, provide either:

data:
  file:   path/to/wp_jackknife.h5
  format: hdf5

or the jackknife patch directory:

fits:
  jk_dir:     data/my_survey/jk_patches/
  jk_pattern: NSIDE_04
  h:          0.6736

Command-Line Options

python hod_mod/scripts/fitting/run_fit_More15.py  <config.yml>  [options]

Positional argument:
  config              Path to YAML configuration file.

Options:
  --map-only          Run MAP optimisation only (skip MCMC even if config says
                      method=both or method=mcmc).
  --mcmc              Run MAP then MCMC sampling (overrides method=map in config).
  --plot-only         Skip fitting; reload fit_result.json and regenerate figures.
  --output-dir DIR    Override the output directory from the config.

Output Files

File

Description

fit_result.json

JSON with label, probes, chi2, ndof, chi2/ndof, success flag, all parameter values.

flatchain.npz

MCMC flat chain (after burn-in). Arrays: flatchain (N×n_free), param_names.

fit_wp.png

\(w_p(r_p)\) data vs MAP prediction with residuals; MCMC bands if available.

fit_ds.png

\(\Delta\Sigma(R)\) data vs MAP prediction with residuals (joint/ESD modes).

fit_combined.png

\(r_p\,w_p\) and \(\Delta\Sigma\) side by side with residuals.

fit_hod.png

\(N_c(M)\) and \(N_s(M)\) at MAP; MCMC occupation credible bands if available.

fit_corner.png

MCMC corner plot with 1σ/2σ contours; reference values from published_params: overlaid.

Parameter Reference

More+2015 HOD Parameters

The More et al. (2015) model has five free parameters describing central and satellite occupation:

Name

Units

Physical meaning

log10mmin

\(\log_{10}(M_\odot/h)\)

Halo mass scale at which 50% of halos host a central galaxy. Controls the overall HOD amplitude and galaxy number density.

sigma_logm

dex

Width of the central occupation step function (log-normal scatter in SHMR). Larger values lower the effective halo mass threshold.

log10m1

\(\log_{10}(M_\odot/h)\)

Characteristic satellite mass scale — halos above \(M_1\) host on average one satellite.

alpha

—

Slope of the satellite mean occupation power law \(N_s \propto (M/M_1)^\alpha\).

kappa

—

Threshold multiplicative factor: satellites only occupy halos above \(\kappa\,M_\mathrm{min}\).

Off-centering Parameters

Fixed by default to the More+2015 published MAP values; can be freed.

Name

Units

Physical meaning

p_off

—

Fraction of central galaxies miscentred from the halo centre.

R_off

\(r_s\) units

Off-centering radial scale in units of the NFW scale radius \(r_s(M)\). A value of 2.2 means off-centred galaxies are displaced by \(2.2\,r_s\) on average.

Incompleteness Parameters

Set free: false and init: 1.0 / 13.0 to disable incompleteness correction.

Name

Units

Physical meaning

alpha_inc

—

Incompleteness power-law slope. 1.0 = no incompleteness.

log10m_inc

\(\log_{10}(M_\odot/h)\)

Halo mass scale for incompleteness onset.

Example: wp-only MAP

Configuration file configs/fitting/More15_wp_example.yml:

# Example: More+2015 HOD model fit to wp(rp) only — BOSS CMASS sample A (logM*>11.1)
# Run:
#   python hod_mod/scripts/fitting/run_fit_More15.py \
#       configs/fitting/More15_wp_example.yml --map-only
#
# Data from More+2015 (arXiv:1407.1856) Table 3 sample A, digitised
# Cosmology fixed to published MAP values

label: "More+2015 BOSS CMASS logM11 — wp-only MAP"

data:
  file:   data/more2015_boss_cmass/logM11_12/wp_logM11_12.csv
  format: bwpd
  rp_min: 0.85
  rp_max: 80.0

cosmology:
  Omega_m: 0.310
  h:       0.703
  sigma8:  0.785
  n_s:     0.964
  Omega_b: 0.0451

model:
  hod_model:   MoreHODModel
  hmf_backend: csst
  z:           0.53
  pi_max:      100.0

parameters:
  log10mmin:  {free: true,  bounds: [12.0, 15.0], init: 13.13}
  sigma_logm: {free: true,  bounds: [0.05, 1.50], init: 0.469}
  log10m1:    {free: true,  bounds: [13.0, 15.5], init: 14.21}
  alpha:      {free: true,  bounds: [0.50, 2.50], init: 1.13}
  kappa:      {free: true,  bounds: [0.01, 3.00], init: 1.25}
  alpha_inc:  {free: false, init: 1.0}
  log10m_inc: {free: false, init: 13.0}
  p_off:      {free: false, init: 0.34}
  R_off:      {free: false, init: 2.2}

fitting:
  method:    map
  n_walkers: 32
  n_steps:   2000
  n_burnin:  500

output:
  dir: results/fitting/More15_wp_example/

# Optional: overlay published best-fit values on plots (dashed reference line)
published_params:
  log10mmin:  [13.13, 0.13]
  sigma_logm: [0.469, 0.13]
  log10m1:    [14.21, 0.12]
  alpha:      [1.130, 0.09]
  kappa:      [1.250, 0.40]

Run:

python hod_mod/scripts/fitting/run_fit_More15.py \
    configs/fitting/More15_wp_example.yml --map-only

Expected output:

============================================================
run_fit_More15  [More+2015 BOSS CMASS logM11 — wp-only MAP]
  Config:   configs/fitting/More15_wp_example.yml
  Probes:   wp
  Free params (5): ['log10mmin', 'sigma_logm', 'log10m1', 'alpha', 'kappa']
  Output:   results/fitting/More15_wp_example/

============================================================
Fit: More+2015 BOSS CMASS logM11 — wp-only MAP
============================================================
  chi2 / ndof = 42.892 / 22  →  chi2/dof = 1.950
  Optimizer:  OK

  Best-fit parameters:
    log10mmin            =   13.191  (free)
    sigma_logm           =   0.5110  (free)
    log10m1              =   14.196  (free)
    alpha                =   1.8605  (free)
    kappa                =   2.6753  (free)
    alpha_inc            =      1.0
    log10m_inc           =     13.0
    p_off                =     0.34
    R_off                =      2.2

Result saved → results/fitting/More15_wp_example/fit_result.json

=== Generating figures ===
  Saved: results/fitting/More15_wp_example/fit_wp.png
  ...

Example: joint MAP + MCMC

Configuration file configs/fitting/More15_joint_example.yml:

# Example: More+2015 HOD model joint fit to wp(rp) + ΔΣ(R) + n_gal
# Run MAP only:
#   python hod_mod/scripts/fitting/run_fit_More15.py \
#       configs/fitting/More15_joint_example.yml --map-only
#
# Run MAP + MCMC:
#   python hod_mod/scripts/fitting/run_fit_More15.py \
#       configs/fitting/More15_joint_example.yml --mcmc
#
# Data from More+2015 (arXiv:1407.1856) sample A (logM*>11.1), digitised
# Cosmology fixed to published MAP values

label: "More+2015 BOSS CMASS logM11 — joint wp+ESD MAP+MCMC"

data:
  file:   data/more2015_boss_cmass/logM11_12/wp_logM11_12.csv
  format: bwpd
  rp_min: 0.85
  rp_max: 80.0

joint:
  ds_file:     data/more2015_boss_cmass/logM11_12/ds_logM11_12.csv
  ds_format:   bwpd_4col
  ds_rp_min:   0.1
  ds_rp_max:   20.0
  ng_obs:      3.0e-4
  ng_frac_err: 0.20

cosmology:
  Omega_m: 0.310
  h:       0.703
  sigma8:  0.785
  n_s:     0.964
  Omega_b: 0.0451

model:
  hod_model:   MoreHODModel
  hmf_backend: csst
  z:           0.53
  pi_max:      100.0

parameters:
  log10mmin:  {free: true,  bounds: [12.0, 15.0], init: 13.13}
  sigma_logm: {free: true,  bounds: [0.05, 1.50], init: 0.469}
  log10m1:    {free: true,  bounds: [13.0, 15.5], init: 14.21}
  alpha:      {free: true,  bounds: [0.50, 2.50], init: 1.13}
  kappa:      {free: true,  bounds: [0.01, 3.00], init: 1.25}
  alpha_inc:  {free: false, init: 1.0}
  log10m_inc: {free: false, init: 13.0}
  p_off:      {free: false, init: 0.34}
  R_off:      {free: false, init: 2.2}

fitting:
  method:    both      # MAP first, then MCMC
  n_walkers: 32
  n_steps:   2000
  n_burnin:  500

output:
  dir: results/fitting/More15_joint_example/

# Optional: overlay published best-fit values on plots (dashed reference line)
published_params:
  log10mmin:  [13.13, 0.13]
  sigma_logm: [0.469, 0.13]
  log10m1:    [14.21, 0.12]
  alpha:      [1.130, 0.09]
  kappa:      [1.250, 0.40]

Run MAP then MCMC sampling:

python hod_mod/scripts/fitting/run_fit_More15.py \
    configs/fitting/More15_joint_example.yml --mcmc

Or run MAP only first to check convergence, then add MCMC:

# Quick MAP check
python hod_mod/scripts/fitting/run_fit_More15.py \
    configs/fitting/More15_joint_example.yml --map-only

# Regenerate figures from saved result
python hod_mod/scripts/fitting/run_fit_More15.py \
    configs/fitting/More15_joint_example.yml --plot-only

The joint likelihood is:

\[-2\ln\mathcal{L} = \chi^2_{w_p} + \chi^2_{\Delta\Sigma} + \chi^2_{n_\mathrm{gal}}\]

where

\[\chi^2_{n_\mathrm{gal}} = \left(\frac{n_\mathrm{gal}^\mathrm{pred} - n_\mathrm{gal}^\mathrm{obs}} {f_\mathrm{err}\,n_\mathrm{gal}^\mathrm{obs}}\right)^2\]

References