Monte Carlo Replica Method Fits

In this section, we describe how to run a Monte Carlo (MC) replica method fit in Colibri.

The Monte Carlo replica method introduces a pseudodata distribution, \(d_p \sim N(d_0, \Sigma)\), which is an approximation to the actual distribution from which the measurement \(d_0\) was drawn. It then defines the corresponding ‘best-fit parameter’ values to be those which minimise the \(\chi^2\) evaluated on the pseudodata:

\[\chi^2 = (d_p - t(\theta))^T \, \Sigma^{-1} \, (d_p - t(\theta))\]

where:

\(d_p\) is the pseudodata vector, drawn from the distribution \(d_p \sim N(d_0, \Sigma)\),
\(t(\theta)\) is the theoretical prediction vector, dependent on the model parameters \(\theta\),
\(\Sigma\) is the covariance matrix encoding the uncertainties.

The idea is to repeat this procedure across multiple data replicas, so that the collection of best-fit models for each pseudodata instance approximates a posterior sample.

Runcard

The runcard runcards/lh_monte_carlo.yaml in the examples/les_houches_example directory performs a Monte Carlo replica fit of the Les Houches model. Note that the dependence on the model comes from the model-specific executable.

################################################################################
# Les Houches model: Monte Carlo replica fit (closure test)
#
# Each run fits one replica, selected with -rep:
#   les_houches_exe runcards/lh_monte_carlo.yaml -rep 1
# Loop over -rep (or submit to a batch system) to get more replicas.
#
# Output folder (created in the current directory): lh_monte_carlo/  (replicas in lh_monte_carlo/fit_replicas/)
# Then run mc_postfit on the folder and evolve the fit.
################################################################################

meta: 'Les Houches model: Monte Carlo replica fit, level 0 closure test on a reduced DIS dataset.'

#######################
# Data and theory specs
#######################

dataset_inputs:
# DIS
- {dataset: SLAC_NC_NOTFIXED_P_EM-F2, variant: legacy_dw}
- {dataset: SLAC_NC_NOTFIXED_D_EM-F2, variant: legacy_dw}
- {dataset: BCDMS_NC_NOTFIXED_P_EM-F2, variant: legacy_dw}
- {dataset: BCDMS_NC_NOTFIXED_D_EM-F2, variant: legacy_dw}
# - {dataset: CHORUS_CC_NOTFIXED_PB_NU-SIGMARED, variant: legacy_dw}
# - {dataset: CHORUS_CC_NOTFIXED_PB_NB-SIGMARED, variant: legacy_dw}
# - {dataset: NUTEV_CC_NOTFIXED_FE_NU-SIGMARED, cfac: [MAS], variant: legacy_dw}
# - {dataset: NUTEV_CC_NOTFIXED_FE_NB-SIGMARED, cfac: [MAS], variant: legacy_dw}
# - {dataset: HERA_NC_318GEV_EM-SIGMARED, variant: legacy}
# - {dataset: HERA_NC_225GEV_EP-SIGMARED, variant: legacy}
# - {dataset: HERA_NC_251GEV_EP-SIGMARED, variant: legacy}
# - {dataset: HERA_NC_300GEV_EP-SIGMARED, variant: legacy}
# - {dataset: HERA_NC_318GEV_EP-SIGMARED, variant: legacy}
# - {dataset: HERA_CC_318GEV_EM-SIGMARED, variant: legacy}
# - {dataset: HERA_CC_318GEV_EP-SIGMARED, variant: legacy}
# - {dataset: HERA_NC_318GEV_EAVG_CHARM-SIGMARED, variant: legacy}
# - {dataset: HERA_NC_318GEV_EAVG_BOTTOM-SIGMARED, variant: legacy}
# - {dataset: NMC_NC_NOTFIXED_EM-F2, variant: legacy_dw}
# - {dataset: NMC_NC_NOTFIXED_P_EM-SIGMARED, variant: legacy}

theoryid: 40000000                     # The theory from which the predictions are drawn.
use_cuts: internal                     # The kinematic cuts to be applied to the data.

closure_test_level: 0                  # False: experimental data; 0: pseudodata without noise;
                                       # 1: pseudodata with noise.
closure_test_pdf: LH_PARAM_20250519    # PDF used to generate the pseudodata (ignored if
                                       # closure_test_level is False).

#####################
# Loss function specs
#####################

positivity:                            # Positivity datasets, used in the positivity penalty.
  posdatasets:
  - {dataset: NNPDF_POS_2P24GEV_F2U, variant: None, maxlambda: 1e6}

positivity_penalty_settings:
  positivity_penalty: false
  alpha: 1e-7
  lambda_positivity: 0

integrability_settings:
  integrability: false

use_t0_covmat: true                    # Use the t0 covariance matrix in the chi2 and to
                                       # generate the MC pseudodata.
t0pdfset: NNPDF40_nnlo_as_01180        # The PDF used to build the t0 covariance matrix.

###################
# Methodology specs
###################

# Pseudodata generation
mcseed: 519562661                      # Base seed for the replica pseudodata.
positive_pseudodata: false             # Resample until all pseudodata points are positive.
mc_validation_fraction: 0.2            # Fraction of data held out for validation. With 0.0
                                       # the training loss is used for early stopping.

# Optimiser
optimizer_settings:
  optimizer: adam                      # Any optax optimizer.
  optimizer_hyperparams:               # Keyword arguments of the optax optimizer.
    learning_rate: 0.01
  # clipnorm: 1.0                      # Gradient clipping by global norm.
  # scheduler:                         # Any optax schedule; overrides learning_rate.
  #   name: linear_schedule
  #   params:
  #     init_value: 1e-2
  #     end_value: 1e-5
  #     transition_begin: 1000
  #     transition_steps: 4000

# Training
max_epochs: 5000                       # Maximum number of training epochs.
patience: 500                          # Epochs without validation improvement before stopping.
# min_delta: 1e-5                      # Minimum change counted as an improvement.

# Data batching (optional). Without batch_size the full training set is used each epoch.
# batch_size: 128
# batch_seed: 1
# shuffle_each_epoch: false

param_initialiser_settings:            # Starting point of each replica (types: uniform,
  type: uniform                        # normal, zeros).
  # random_seed: 0                     # Added to the replica index.
  bounds:
    alpha_gluon: [-0.1, 1]
    beta_gluon: [9, 13]
    alpha_up: [0.4, 0.9]
    beta_up: [3, 4.5]
    epsilon_up: [-3, 3]
    gamma_up: [1, 6]
    alpha_down: [1, 2]
    beta_down: [8, 12]
    epsilon_down: [-4.5, -3]
    gamma_down: [3.8, 5.8]
    norm_sigma: [0.1, 0.5]
    alpha_sigma: [-0.2, 0.1]
    beta_sigma: [1.2, 3]

actions_:
- run_monte_carlo_fit

`optimizer_settings`

These settings control the method of gradient descent. You can use any of the Optax optimizers and settings, which you can read more about here. Learning schedulers are also supported, and you can find the available options here.

`param_initialiser_settings`

If it is appropriate for a given model, you may choose to have single, global minimum and maximum values for all parameters, instead of specific bounds for each parameter. In that, case, you may replace

bounds:
    alpha_gluon: [-0.1, 1]
    beta_gluon: [9, 13]
...

with, for example:

min_val: -4.5
max_val: 13

Initialising parameters with a normal distribution

If you want to initialise with parameters that are distributed normally rather than uniformally (as in the runcard above) you may do so with the following settings:

param_initialiser_settings:               # The initialiser for Monte Carlo training.
    type: normal
    means: 2.0
    stds: 0.5

This will set the mean to 2.0 and standard deviation to 0.5 for all parameters. If instead you want to specify different values for each parameter, either for the means or standard deviations or both, you can do so:

param_initialiser_settings:               # The initialiser for Monte Carlo training.
    type: normal
    means:
        alpha_gluon: 0.356
        beta_gluon: 10.9
        [...]
    stds:
        alpha_gluon: 0.06
        beta_gluon: 0.55
        [...]

Note

If you use initialise with a normal distribution and don’t specify means or standard deviations, default values of 0.0 and 1.0 will be used respectively.

Note

param_initialiser_settings covers the PDF parameters only. If the forward map has extra parameters (e.g. SMEFT Wilson coefficients), they are initialised by extra_param_initialiser_settings, with the same options, and to zero by default. See Extra parameters of the forward map.

Using data batching

In Monte Carlo replica fits, it is possible to use data batching during training. This can be done by adding the following settings to the runcard:

batch_size: 128
batch_seed: 3
shuffle_each_epoch: False

batch_size: The size of the data batches. If the dataset is smaller than this size or the batch_size is not specified, the full dataset will be used.
batch_seed: The random seed used to generate the data batches.
shuffle_each_epoch: Whether to reshuffle the data at the start of each epoch. If set to False, the data will be shuffled only once at the start of training and the batches will be fixed.

Running the fit

To perform a Monte Carlo fit for the Les Houches model (presented in this tutorial), you would run:

pip install -e .

from the examples/les_houches_example directory.

Then, from the same directory, you can run the above runcard with the following command:

les_houches_exe runcards/lh_monte_carlo.yaml -rep N

Note that this command will generate one single replica, namely replica number N. For a fit with more than one replica, you should iterate the above or submit the job to a batch system.

You can then run a postfit selection on the fit folders and evolve the fit. Details on how to do this can be found in this section.