pyfc package

Submodules

pyfc.binned module

Binned Likelihood and Feldman-Cousins Toy Generation Module

This file contains the core mathematical and statistical routines required for performing binned frequentist analyses, tailored for the Feldman-Cousins construction. It implements expected rate calculations, standard Poisson and finite-Monte Carlo (Poisson-Gamma mixture) likelihood evaluations, grid-based profile likelihood optimizations, and parametric bootstrap pseudo-experiment (toy) generation.

Created: v0.1.0 (July 24, 2026) Author: Mauricio Bustamante (mbustamante@gmail.com)

This file was released as part of the PyFC code, stored at https://github.com/mbustama/FeldmanCousins, which exists under a GNU GPL v3 License.

pyfc.binned.calc_nll(params, N_obs, S_sumw2, B_sumw2, use_finite_mc, compute_rates_func)

[BINNED] Computes the negative log-likelihood for arbitrary N parameters.

Statistical Theory & Assumptions: This function computes the negative log-likelihood (NLL) for comparing observed data counts against expected model counts. It handles two distinct statistical regimes:

1. Standard Poisson Likelihood (use_finite_mc = False): Assumes template predictions have zero statistical uncertainty (infinite MC). It uses the Baker-Cousins chi-square equivalent form (based on the likelihood ratio with a saturated model). Expression: NLL = 2 * SUM [ mu_i - n_i + n_i * ln(n_i / mu_i) ]

2. Finite Monte Carlo Likelihood (use_finite_mc = True): Accounts for statistical uncertainty in the simulation templates. It marginalizes over the unknown true Poisson rate using a Gamma conjugate prior, resulting in a Negative Binomial (Poisson-Gamma mixture) likelihood. Expression: alpha = (mu^2 / sigma^2) + 1 beta = mu / sigma^2 ln(L) = alpha*ln(beta) + ln(Gamma(n+alpha)) - (n+alpha)*ln(1+beta) - ln(Gamma(alpha)) NLL = -2 * ln(L)

   This is exactly L_Eff, Eq. (3.16) of Arguelles, Schneider & Yuan, “A binned likelihood for stochastic models”, JHEP 06 (2019) 030 [arXiv:1901.04645] (their main, recommended result; see also docs/source/refs.bib’s Arguelles2019izp entry and docs/source/methodology.rst). alpha/beta here are their Eq. (3.12) with the uniform-prior choice a=1, b=0 in their Eq. (3.17) – the “+1” in alpha is a deliberate, load-bearing part of that choice, not a typo: it is what distinguishes L_Eff from the paper’s alternative L_Mean parameterization (a=b=0, alpha = mu^2/sigma^2 with no “+1”), which matches mean-and-variance moment-matching but has worse coverage properties in the paper’s own toy-experiment tests (their Sec. 4.2, Fig. 5). The paper’s own ln(k!) term is omitted here, as it is a data-only constant that cancels in any NLL difference (the same convention documented in unbinned.calc_nll_unbinned’s own dropped ln(N_obs!) term).

   Numerically stable evaluation (not the naive formula above): for well-simulated templates (small sigma^2 relative to mu^2), alpha grows large, and evaluating ln(Gamma(n+alpha)) - ln(Gamma(alpha)) as a literal difference of two lgamma calls loses catastrophic precision – both terms individually scale like alpha*ln(alpha) while their difference is only O(n*ln(alpha)), so double precision’s ~16 significant digits get eaten by the leading digits both terms share. Verified directly against a 50-digit mpmath reference: at mu=1000, sigma^2=1e-9 (a realistic well-simulated-template regime), the naive difference is off by 7.5 in absolute NLL – large enough to distort t = NLL_cond - NLL_uncond and bias interval boundaries. Fixed via two exact reformulations, not an approximation: (a) ln(Gamma(n+alpha)) - ln(Gamma(alpha)) is computed as the exact integer-recurrence sum SUM_{k=0}^{n-1} ln(alpha + k) (valid since n_obs is always a non-negative integer count), which adds n well-conditioned terms instead of subtracting two enormous ones; (b) alpha*ln(beta) - (n+alpha)*ln(1+beta) is rewritten as -alpha*log1p(1/beta) - n*log1p(beta), avoiding the same kind of cancellation between ln(beta) and ln(1+beta) when beta is large. Both reformulations are algebraically exact (zero added approximation error), not truncated series – see tests/test_finite_mc_likelihood.py’s high-alpha regression tests, which fail against the naive formula and pass against this one.

   Note on the sigma2 <= 1e-10 Poisson fallback below: this is not a numerical workaround for the instability above – it is the mathematically correct limit. As sigma^2 -> 0, alpha -> infinity and this Negative Binomial collapses onto Poisson(mu) (zero simulation variance is exactly the “infinite MC” assumption the plain-Poisson branch above already encodes), so evaluating the finite-MC formula at sigma^2 = 0 exactly would only divide by zero, not disagree with the limit. The two branches use different additive normalizations (this branch drops only ln(n!); the use_finite_mc=False branch is the fully saturated Baker-Cousins form), so their raw NLL values differ by a data-only constant even in this limit – only t = NLL_cond - NLL_uncond is meaningful to compare, and that constant cancels there regardless of which branch is taken, as long as both the conditional and unconditional fit for a given call use the same use_finite_mc setting (they always do).

Parameters:

paramsarray_like, 1D

The physical model parameters to evaluate.

N_obsarray_like, any shape

The observed data counts (or toy data) in each bin. May be an arbitrary N-dimensional histogram (e.g. shape (E_bins, cos_theta_bins)); it is flattened internally, so the NLL is a sum over all bins regardless of how they are laid out spatially. Non-contiguous views (e.g. a transposed histogram) are handled by copying to a contiguous layout before flattening, since numba’s nopython .reshape(-1) only accepts C-contiguous input.

S_sumw2array_like, same shape as N_obs

Sum of squared MC weights (sum(w_i^2)) per bin for the signal template – the standard per-bin variance estimator for a weighted MC sample (unweighted MC is the w_i=1 special case).

B_sumw2array_like, same shape as N_obs

Sum of squared MC weights per bin for the background template.

use_finite_mcbool

If True, applies the Poisson-Gamma mixture likelihood. If False, uses standard Poisson.

compute_rates_funccallable

User-provided mapping function (params, S_sumw2, B_sumw2) -> (mu, sigma2) that calculates bin expectations (mu) and variances (sigma2), of the same shape as N_obs. Any fixed template arrays it needs (signal/background shapes) should be referenced via closure or module-global rather than passed in as arguments.

Returns:

nllfloat

The calculated negative log-likelihood value. For bins where the model yields an unphysical or near-zero expected count mu_i (<= 1e-12) while data was observed there (n_obs > 0), a smooth quadratic barrier is added to that bin’s contribution instead of returning a flat constant (see note below).

Note on unphysical (mu_i <= mu_floor = 1e-12) bins:

Earlier versions of this function returned a flat penalty (1e10) the instant any bin’s expected count went non-positive, discarding whatever NLL had already been accumulated from other bins in this call. Because scipy.optimize.minimize(…, method=’L-BFGS-B’) estimates gradients via finite differences, a locally-constant NLL surface reads as “gradient ~ 0, already converged,” which can permanently trap the optimizer if its starting guess lands in a mu <= 0 region (e.g. the bounds midpoint, for a user model that is unphysical there). Instead, the contribution from an unphysical bin is now a continuous extension of the real Poisson NLL term (evaluated at a tiny positive floor) plus a quadratic barrier that grows with how far mu_i actually is below that floor, so the gradient w.r.t. mu_i stays informative and points back toward mu_i > 0. This bin-local contribution is added to the running total rather than short-circuiting the whole function. The trigger is mu_i <= mu_floor rather than mu_i <= 0 specifically so there is no razor-thin discontinuity in the tiny sliver (0, mu_floor): the real Poisson term diverges as mu_i -> 0+, while the floor-based term is constant (evaluated at mu_floor, not at the actual mu_i), so triggering only at mu_i <= 0 would make the NLL drop right at that boundary instead of continuing to climb. For all mu_i > mu_floor (i.e. all physically meaningful values – mu_floor is astronomically small), behavior is practically identical to before this change (matches to within 1e-12 relative/absolute tolerance; see tests/test_smoothing.py).

pyfc.binned.conditional_fit_grid_1d(test_val, fix_idx, n_params, N_obs, cond_grid_points, S_sumw2, B_sumw2, use_finite_mc, compute_rates_func)

Performs a conditional maximum likelihood fit with one parameter fixed (Profiling).

Statistical Theory: To construct the numerator of the Profile Likelihood Ratio for a 1D confidence interval, the likelihood is maximized (NLL minimized) while keeping the parameter of interest fixed to a specific test value. The remaining (nuisance) parameters are allowed to vary (profiled out).

Expression: NLL_cond = min( NLL(theta_test, nu) ) over all nuisance parameters nu.

Parameters:

test_valfloat

The fixed value for the parameter of interest.

fix_idxint

The index of the parameter to fix.

n_paramsint

Total number of parameters in the model.

N_obs, S_sumw2, B_sumw2array_like, any shape (all matching)

Data counts and variances.

cond_grid_pointsarray_like, 2D

Pre-computed grid of the (N-1) free nuisance parameters to scan over.

use_finite_mcbool

Flag to use finite MC likelihood formulation.

compute_rates_funccallable

User-provided physics mapping function.

Returns:

min_nllfloat

The minimum negative log-likelihood given the fixed parameter.

best_paramsarray_like, 1D

The full parameter combination (including the fixed one) that produced min_nll.

pyfc.binned.conditional_fit_grid_2d(test_vA, test_vB, fix_A, fix_B, n_params, N_obs, cond_grid_points, S_sumw2, B_sumw2, use_finite_mc, compute_rates_func)

Performs a conditional maximum likelihood fit with two parameters fixed.

Statistical Theory: Analogous to the 1D case, this computes the profiled likelihood for a 2D joint confidence region. Two parameters of interest are fixed to a grid point, and the NLL is minimized over all remaining (N-2) nuisance parameters.

Parameters:

test_vA, test_vBfloat

The fixed values for the two parameters of interest.

fix_A, fix_Bint

The indices of the parameters being fixed.

n_paramsint

Total number of parameters in the model.

N_obs, S_sumw2, B_sumw2array_like, any shape (all matching)

Data counts and variances.

cond_grid_pointsarray_like, 2D

Pre-computed grid of the (N-2) free nuisance parameters to scan over.

use_finite_mcbool

Flag to use finite MC likelihood formulation.

compute_rates_funccallable

User-provided physics mapping function.

Returns:

min_nllfloat

The minimum negative log-likelihood given the two fixed parameters.

best_paramsarray_like, 1D

The full parameter combination that produced min_nll.

pyfc.binned.generate_and_fit_toys_grid_1d(test_val, fix_idx, true_params, n_params, full_grid_points, cond_grid_points, n_toys, S_sumw2, B_sumw2, use_finite_mc, compute_rates_func)

Generates Poisson toys and computes the 1D test statistic distribution (Feldman-Cousins).

Statistical Theory & Assumptions: To ensure proper frequentist coverage (as per the Feldman-Cousins unified approach), we cannot assume the Profile Likelihood Ratio strictly follows an asymptotic chi-square distribution (Wilks’ Theorem), especially near physical boundaries. Instead, we build the empirical distribution of the test statistic (t) by running parametric bootstraps (toys).

For each toy: 1. A synthetic dataset (toy_N) is sampled from Poisson(mu_true). 2. The unconditional global minimum NLL is found for the toy. 3. The conditional minimum NLL (fixing the parameter of interest) is found. 4. The test statistic is calculated: t = NLL_cond - NLL_uncond. (Due to numerical imprecision in grid searches, t is bounded at a minimum of 0.0).

Parameters:

test_valfloat

The fixed value for the parameter of interest.

fix_idxint

The index of the fixed parameter.

true_paramsarray_like, 1D

The parameter combination used to generate the true expected counts for the toys.

n_paramsint

Total number of parameters.

S_sumw2, B_sumw2array_like, any shape

Variances (matching whatever shape compute_rates_func’s mu produces).

full_grid_pointsarray_like, 2D

Grid for the unconditional fit.

cond_grid_pointsarray_like, 2D

Grid for the conditional (nuisance parameter) fit.

n_toysint

The number of pseudo-experiments to generate.

use_finite_mcbool

Flag to use finite MC likelihood during toy fitting.

compute_rates_funccallable

User-provided physics mapping function.

Returns:

t_statisticsarray_like, 1D

Array of length n_toys containing the calculated test statistic for each toy.

pyfc.binned.generate_and_fit_toys_grid_2d(test_vA, test_vB, fix_A, fix_B, true_params, n_params, full_grid_points, cond_grid_points, n_toys, S_sumw2, B_sumw2, use_finite_mc, compute_rates_func)

Generates Poisson toys and computes the 2D test statistic distribution.

Statistical Theory: Similar to the 1D case, this function calculates empirical distributions of the Profile Likelihood Ratio test statistic for 2D parameter joint combinations. t = NLL_cond_2D - NLL_uncond.

Parameters:

test_vA, test_vBfloat

The fixed values for the two parameters of interest.

fix_A, fix_Bint

The indices of the parameters being fixed.

true_paramsarray_like, 1D

The parameter combination used to generate expected counts for the toys.

n_paramsint

Total number of parameters.

S_sumw2, B_sumw2array_like, any shape

Variances (matching whatever shape compute_rates_func’s mu produces).

full_grid_pointsarray_like, 2D

Grid for the unconditional fit.

cond_grid_pointsarray_like, 2D

Grid for the conditional (nuisance parameter) fit.

n_toysint

The number of pseudo-experiments to generate.

use_finite_mcbool

Flag to use finite MC likelihood during toy fitting.

compute_rates_funccallable

User-provided physics mapping function.

Returns:

t_statisticsarray_like, 1D

Array of length n_toys containing the calculated test statistic for each toy.

pyfc.binned.unconditional_fit_grid(N_obs, full_grid_points, S_sumw2, B_sumw2, use_finite_mc, compute_rates_func)

Performs an unconditional maximum likelihood fit via exhaustive grid search.

Statistical Theory: To construct the denominator of the Profile Likelihood Ratio (PLR) test statistic, we must find the global maximum of the likelihood function (minimum of the NLL) allowing all parameters in the model to vary freely across their physical bounds.

Expression: NLL_best = min( NLL(theta) ) for all theta in parameter space.

Parameters:

N_obsarray_like, any shape

The observed data or toy counts.

full_grid_pointsarray_like, 2D

A pre-computed matrix where each row represents a complete N-dimensional parameter combination to evaluate.

S_sumw2, B_sumw2array_like, same shape as N_obs

Sum of squared MC weights (sumw2) per bin, for signal and background respectively.

use_finite_mcbool

Flag to use finite MC likelihood formulation.

compute_rates_funccallable

User-provided physics mapping function.

Returns:

min_nllfloat

The absolute minimum negative log-likelihood found in the grid.

best_paramsarray_like, 1D

The parameter combination that produced min_nll.

pyfc.config module

Configuration and Argument Parsing Module

This file defines the configuration management system for the PyFC framework. It establishes a hierarchical parameter loading mechanism (Hardcoded Defaults -> JSON Configuration File -> Command Line Arguments) to control the execution of the Feldman-Cousins confidence interval construction.

While this file does not execute the mathematical routines, it defines the statistical, algorithmic, and computational parameters that govern the underlying profile likelihood ratio tests and Monte Carlo toy generation.

Created: v0.1.0 (July 24, 2026) Author: Mauricio Bustamante (mbustamante@gmail.com)

This file was released as part of the PyFC code, stored at https://github.com/mbustama/FeldmanCousins, which exists under a GNU GPL v3 License.

pyfc.config.generate_sample_config(filename='../config/example_fc_config.json')

Generates a sample JSON configuration file containing default execution parameters.

Statistical & Algorithmic Context: The default configuration initializes a standard binned likelihood analysis calculating confidence intervals at 68% and 90% Confidence Levels (CL). It enables the finite Monte Carlo correction (assuming template statistical uncertainties require a Poisson-Gamma mixture likelihood) and utilizes an adaptive Monte Carlo toy generation strategy to optimize computational resources when empirical p-values are far from the critical threshold alpha.

Parameters:

filenamestr, optional

The destination file path where the JSON configuration will be written. Defaults to “../config/example_fc_config.json”.

Returns:

None

The function writes a file to disk and prints a confirmation message.

pyfc.config.parse_arguments()

Parses execution parameters using a strict hierarchy: Defaults -> JSON Config File -> CLI Arguments.

Statistical Theory & Configuration Parameters: - likelihood_type: Defines the probability density function (PDF) form.

‘binned’ uses Poisson/Negative Binomial counting statistics; ‘unbinned’ uses continuous event likelihoods.

• cl: The target Confidence Level(s) (1 - alpha), defining the required frequentist coverage probability (e.g., 0.90 for 90% coverage).

• n_toys: The baseline number of parametric bootstrap pseudo-experiments used to build the empirical Profile Likelihood Ratio (PLR) distribution.

• use_finite_mc_correction_binned: If True, convolves the standard Poisson likelihood with a Gamma prior to account for limited statistics in MC templates.

• strategy: The optimization algorithm used to minimize the Negative Log-Likelihood (NLL). Options include exhaustive ‘grid’, gradient-based ‘scipy’, or nested sampling ‘ultranest’.

Parameters:

None (Reads directly from command line via sys.argv)

Returns:

configdict

A dictionary containing the fully resolved configuration parameters necessary to execute the Feldman-Cousins orchestrator.

pyfc.generate_config module

Configuration Generator Module (Interactive CLI)

This file contains a standalone interactive command-line interface (CLI) designed to generate a properly formatted JSON configuration file (fc_config.json) for the PyFC (Feldman-Cousins) orchestrator. It handles user input collection, type casting, rule validation, and default parameter fallback.

While this script does not perform mathematical calculations directly, it defines the statistical hyper-parameters (such as Confidence Levels, likelihood formulations, and Monte Carlo toy counts) that dictate the rigorousness and underlying assumptions of the frequentist confidence interval construction.

Usage (Command Line Interface):

To launch the interactive configuration generator, run the script directly from your terminal:

$ python generate_config.py

Alternatively, if the PyFC package is installed in your Python environment:

$ python -m pyfc.generate_config

Interactive Instructions: - You will be prompted sequentially for various configuration values. - Type your desired value and press [Enter]. - To accept the default value displayed in brackets (e.g., [Default: binned]),

simply press [Enter] without typing anything.

• For boolean (yes/no) questions, accepted inputs include ‘y’, ‘yes’, ‘t’, ‘true’, ‘1’ for True, and ‘n’, ‘no’, ‘f’, ‘false’, ‘0’ for False (case-insensitive).

• For list inputs (like Confidence Levels), provide values separated by commas (e.g., 0.68, 0.90, 0.95).

Created: v0.1.0 (July 24, 2026) Author: Mauricio Bustamante (mbustamante@gmail.com)

This file was released as part of the PyFC code, stored at https://github.com/mbustama/FeldmanCousins, which exists under a GNU GPL v3 License.

pyfc.generate_config.get_input(prompt_text, default_val, cast_func, choices=None, validator=None, error_msg='Invalid input.')

Core interactive prompt loop for CLI data entry.

Handles rendering the prompt, intercepting empty inputs for default fallback, parsing the input via the provided casting function, and executing all validation rules (like bounds checking for statistical parameters).

Parameters:

prompt_textstr

The question or prompt presented to the user.

default_valany

The fallback value used if the user submits an empty response (presses Enter).

cast_funccallable

A function (like int, str, or parse_bool) used to transform the string input into the desired target data type.

choiceslist, optional

A specific list of allowed values. If provided, input must exist in this list.

validatorcallable, optional

A function that takes the parsed value and returns True if valid, False otherwise. Used for mathematical bounds (e.g., ensuring n_toys > 0).

error_msgstr, optional

The message displayed if the validator function fails.

Returns:

parsed_valany

The fully validated, type-casted user input (or the default value).

pyfc.generate_config.main()

Main entry point for the configuration generator.

Statistical Theory & Configuration Impact: This sequence builds the fc_config.json matrix. It defines: - PDF formulation (Binned vs Unbinned) which changes the NLL definition. - Confidence Levels (CL) which sets the coverage integration targets. - Monte Carlo sampling size (n_toys), controlling the resolution of the empirical

test statistic PDF. Higher toys reduce statistical noise in the p-value calculation at the expense of computation time.

• Likelihood optimization strategies (SciPy, UltraNest, etc.) used to find the global and conditional minima of the parameter space during profiling.

• Advanced topological settings (adaptive toys, grid sparsification) which employ heuristics to skip unnecessary toy generations in regions definitively outside the targeted confidence bounds.

Parameters:

None

Returns:

None

Outputs a serialized JSON file containing the validated configuration.

pyfc.generate_config.parse_bool(value)

Safely converts string inputs to boolean values.

Statistical Context: Boolean flags in the configuration control critical methodological choices, such as whether to apply the finite Monte Carlo correction (which shifts the likelihood from a standard Poisson distribution to a Negative Binomial Poisson-Gamma mixture) or whether to use adaptive sampling for the toys.

Parameters:

valuestr or bool

The raw input provided by the user (e.g., ‘y’, ‘n’, ‘true’, ‘1’).

Returns:

bool

The parsed boolean True or False.

Raises:

ValueError:

If the input string cannot be mapped to a valid boolean state.

pyfc.generate_config.parse_float_list(value)

Converts a comma-separated string to a list of floats.

Statistical Context: This is primarily used to parse the Confidence Levels (CL). In the Feldman-Cousins approach, the confidence level (1 - alpha) determines the critical value of the Profile Likelihood Ratio test statistic. The framework will compute intervals ensuring exact frequentist coverage at these specified probabilities (e.g., 0.68 for 1-sigma, 0.90 for 90% CL).

Parameters:

valuestr or list

The raw input string containing comma-separated numbers (e.g., “0.68, 0.90”).

Returns:

list of float

A list of floating-point numbers parsed from the input.

pyfc.generate_config.parse_str_list(value)

Converts a comma-separated string to a list of strings.

Parameters:

valuestr or list

The raw input string containing comma-separated names.

Returns:

list of str

A list of stripped string values (typically representing physical parameter names).

pyfc.optimizers module

Continuous Optimization Module

This file contains the continuous optimization routines used to minimize the Negative Log-Likelihood (NLL) functions for both binned and unbinned models. It provides interfaces for gradient-based local optimization via SciPy (L-BFGS-B) and robust global optimization via UltraNest (Nested Sampling). These optimizers are essential for profiling out nuisance parameters and finding the global and conditional likelihood maxima required for the Profile Likelihood Ratio test statistic in the Feldman-Cousins construction.

A note on joint/simplex-constrained parameters:

Every bounds_list entry here is an INDEPENDENT per-parameter box constraint only – param_i in [lo_i, hi_i], with no notion of a relationship between two different parameters. Neither L-BFGS-B/SLSQP nor UltraNest’s prior transform (which maps a unit hypercube onto bounds_list one coordinate at a time) have any native concept of a joint constraint like a simplex (a + b <= 1) or a sphere (a^2 + b^2 <= 1). If your physical model has such a constraint, do not work around it with a hand-rolled penalty function multiplying your rate – that reproduces the exact flat-gradient trap described in calc_nll’s docstring, just user-side instead of library-side. Two purpose-built mechanisms exist instead, and can be used together:

• bounds_func (accepted by every function below): tightens a free parameter’s own box bound based on whatever OTHER parameter(s) the current scan step has fixed. Use this when the constraint only ever couples the scan’s currently-FIXED test parameter(s) to a free nuisance parameter – it’s cheap, since it keeps the optimizer’s search box (and its bounds-midpoint starting guess) always physical, with no boundary- adjacent evaluations needed at all.

• constraints (accepted by unconditional_fit_scipy, the two conditional_fit_*_scipy functions, and the two conditional_fit_*_ ultranest functions besides unconditional_fit_ultranest): a list of scipy.optimize.LinearConstraint/NonlinearConstraint objects, expressed in the FULL n_params space. Use this when the constraint couples multiple parameters that are SIMULTANEOUSLY free (e.g. neither is ever the scan’s fixed test parameter) – a case bounds_func cannot express. On the scipy side this switches method from ‘L-BFGS-B’ to ‘SLSQP’ (the only two of PyFC’s methods that support constraints together with bounds); on the UltraNest side, a violated constraint just yields a very large negative log-likelihood at that point, since nested sampling doesn’t use gradients and isn’t vulnerable to the flat-region trap that motivates the scipy-side machinery.

See the README section “Handling Joint/Simplex-Constrained Parameters” for a full worked example (a neutrino flavor-fraction fit with f_e + f_mu <= 1).

Created: v0.1.0 (July 24, 2026) Author: Mauricio Bustamante (mbustamante@gmail.com)

This file was released as part of the PyFC code, stored at https://github.com/mbustama/FeldmanCousins, which exists under a GNU GPL v3 License.

pyfc.optimizers.conditional_fit_1d_scipy(test_val, fix_idx, n_params, data, bounds_list, compute_rates_func, seed=None, pdf_components=None, likelihood_type='binned', S_sumw2=None, B_sumw2=None, use_finite_mc=False, bounds_func=None, constraints=None, scipy_method=None, n_restarts=1)

Performs a conditional maximum likelihood fit (profiling 1 parameter) using SciPy.

Statistical Theory: To evaluate the numerator of the Profile Likelihood Ratio for 1D intervals, one parameter of interest is fixed to a specific ‘test_val’, and the NLL is minimized over the remaining (N-1) nuisance parameters.

Mathematical Expression: NLL_cond = min_{nu} NLL(theta_fixed, nu | data) where theta_fixed is the parameter of interest and nu are the nuisance parameters.

Parameters:

test_valfloat

The value at which to fix the parameter of interest.

fix_idxint

The index of the parameter being fixed.

n_paramsint

Total number of parameters.

data, bounds_listvarious

Dataset (binned counts of any shape, or unbinned events) and boundary constraints.

compute_rates_funccallable

User-provided mapping function matching parameters to physical expectations.

seedarray_like, 1D, optional

Initial guess for the free parameters.

pdf_componentslist of callable, optional

Probability density functions, one per model component. Required when likelihood_type=”unbinned”; unused for “binned”.

likelihood_type, S_sumw2, B_sumw2, use_finite_mcvarious

Likelihood formulation configurations.

bounds_funccallable, optional

Advanced hook for expressing a joint/simplex constraint between the currently-fixed test_val and the free nuisance parameters (see module docstring). If provided, called as bounds_func({fix_idx: test_val}, free_indices, bounds_list) and must return one (lo, hi) tuple per entry of free_indices, in that order. This lets the box itself be tightened so it is always physical (e.g. f_mu <= 1 - f_e_test), which also fixes the default bounds-midpoint starting guess so it can never land in forbidden territory. If None (default), each free parameter uses its static bounds_list[i] unmodified – identical to the behavior before bounds_func existed.

constraintslist of scipy.optimize.LinearConstraint/NonlinearConstraint, optional

Joint constraints among the FULL n_params vector (see module docstring and _project_linear_constraint). Since fix_idx is not part of this function’s optimization vector, each constraint is projected down to the free-parameter subspace before being passed to SciPy, substituting test_val as a constant offset. When non-empty, method defaults to ‘SLSQP’ instead of ‘L-BFGS-B’. When None/empty (default), behavior and method are UNCHANGED from before this parameter existed. Use constraints (rather than bounds_func) for constraints among multiple simultaneously-FREE nuisance parameters, which bounds_func cannot express.

scipy_methodstr, optional

Explicit override for the scipy.optimize.minimize method.

n_restartsint, optional

Number of distinct starting points to try, keeping whichever converges to the lowest NLL (see _minimize_with_restarts). Default 1 preserves the original single-fit behavior (a single fit from x0), except that res.success is now always checked; a non-converged result triggers one cheap perturbed retry and a warning if that also fails.

Returns:

min_nllfloat

The profiled minimum negative log-likelihood.

best_parray_like, 1D

The full parameter array (size N) including the fixed test_val and the optimized nuisance parameters.

pyfc.optimizers.conditional_fit_1d_ultranest(test_val, fix_idx, n_params, data, bounds_list, compute_rates_func, verbose=1, pdf_components=None, likelihood_type='binned', S_sumw2=None, B_sumw2=None, use_finite_mc=False, bounds_func=None, constraints=None)

Performs a conditional maximum likelihood fit (profiling 1 parameter) using UltraNest.

Statistical Theory: Executes Reactive Nested Sampling in the reduced (N-1) dimensional subspace of nuisance parameters to robustly find the conditional maximum likelihood. This prevents false confidence interval exclusions that can occur if a gradient optimizer falls into a local minimum when profiling.

Parameters:

test_valfloat

The value at which to fix the parameter of interest.

fix_idxint

The index of the parameter being fixed.

n_paramsint

Total number of parameters in the model.

data, bounds_listvarious

Dataset (binned counts of any shape, or unbinned events) and boundary constraints.

compute_rates_funccallable

User-provided mapping function matching parameters to physical expectations.

verboseint

Verbosity level.

pdf_componentslist of callable, optional

Probability density functions, one per model component. Required when likelihood_type=”unbinned”; unused for “binned”.

likelihood_type, S_sumw2, B_sumw2, use_finite_mcvarious

Likelihood formulation configurations.

bounds_funccallable, optional

Advanced hook for expressing a joint/simplex constraint between the currently-fixed test_val and the free nuisance parameters (see module docstring). If provided, called as bounds_func({fix_idx: test_val}, free_indices, bounds_list) and must return one (lo, hi) tuple per entry of free_indices, in that order. If None (default), each free parameter uses its static bounds_list[i] unmodified – identical to the behavior before bounds_func existed.

constraintslist of scipy.optimize.LinearConstraint/NonlinearConstraint, optional

Joint constraints among the full n_params vector, evaluated directly on the fully-assembled parameter point (no projection needed – nested sampling doesn’t use gradients, so a violated constraint simply makes log_likelihood return -1e10 at that point). If None/empty (default), behavior is unchanged.

Returns:

min_nllfloat

The profiled minimum negative log-likelihood.

best_parray_like, 1D

The full parameter array (size N) including the fixed test_val and the optimized nuisance parameters.

pyfc.optimizers.conditional_fit_2d_scipy(test_vA, test_vB, fix_A, fix_B, n_params, data, bounds_list, compute_rates_func, seed=None, pdf_components=None, likelihood_type='binned', S_sumw2=None, B_sumw2=None, use_finite_mc=False, bounds_func=None, constraints=None, scipy_method=None, n_restarts=1)

Performs a conditional maximum likelihood fit (profiling 2 parameters) using SciPy.

Statistical Theory: Analogous to the 1D profiling function, but computes the joint conditional minimum for a 2D confidence region by fixing two parameters and optimizing the remaining (N-2) nuisance parameters.

Parameters:

test_vA, test_vBfloat

The values at which to fix the two parameters of interest.

fix_A, fix_Bint

The indices of the parameters being fixed.

n_paramsint

Total number of parameters.

data, bounds_listvarious

Dataset (binned counts of any shape, or unbinned events) and boundary constraints.

compute_rates_funccallable

User-provided mapping function matching parameters to physical expectations.

seedarray_like, 1D, optional

Initial guess for the free parameters.

pdf_componentslist of callable, optional

Probability density functions, one per model component. Required when likelihood_type=”unbinned”; unused for “binned”.

likelihood_type, S_sumw2, B_sumw2, use_finite_mcvarious

Likelihood formulation configurations.

bounds_funccallable, optional

Advanced hook for expressing a joint/simplex constraint between the two currently-fixed test parameters and the free nuisance parameters (see module docstring). If provided, called as bounds_func({fix_A: test_vA, fix_B: test_vB}, free_indices, bounds_list) and must return one (lo, hi) tuple per entry of free_indices, in that order. If None (default), each free parameter uses its static bounds_list[i] unmodified – identical to the behavior before bounds_func existed.

constraintslist of scipy.optimize.LinearConstraint/NonlinearConstraint, optional

Joint constraints among the FULL n_params vector, projected down to the free-parameter subspace (substituting test_vA/test_vB as a constant offset) before being passed to SciPy. When non-empty, method defaults to ‘SLSQP’ instead of ‘L-BFGS-B’. When None/empty (default), behavior and method are UNCHANGED from before this parameter existed.

scipy_methodstr, optional

Explicit override for the scipy.optimize.minimize method.

n_restartsint, optional

Number of distinct starting points to try, keeping whichever converges to the lowest NLL (see _minimize_with_restarts). Default 1 preserves the original single-fit behavior (a single fit from x0), except that res.success is now always checked; a non-converged result triggers one cheap perturbed retry and a warning if that also fails.

Returns:

min_nllfloat

The profiled minimum negative log-likelihood.

best_parray_like, 1D

The full parameter array (size N) including the two fixed values and the optimized nuisance parameters.

pyfc.optimizers.conditional_fit_2d_ultranest(test_vA, test_vB, fix_A, fix_B, n_params, data, bounds_list, compute_rates_func, verbose=1, pdf_components=None, likelihood_type='binned', S_sumw2=None, B_sumw2=None, use_finite_mc=False, bounds_func=None, constraints=None)

Performs a conditional maximum likelihood fit (profiling 2 parameters) using UltraNest.

Statistical Theory: Executes Reactive Nested Sampling in the (N-2) dimensional nuisance parameter subspace to construct a robust 2D confidence region profile.

Parameters:

test_vA, test_vBfloat

The values at which to fix the two parameters of interest.

fix_A, fix_Bint

The indices of the parameters being fixed.

n_paramsint

Total number of parameters in the model.

data, bounds_listvarious

Dataset (binned counts of any shape, or unbinned events) and boundary constraints.

compute_rates_funccallable

User-provided mapping function matching parameters to physical expectations.

verboseint

Verbosity level.

pdf_componentslist of callable, optional

Probability density functions, one per model component. Required when likelihood_type=”unbinned”; unused for “binned”.

likelihood_type, S_sumw2, B_sumw2, use_finite_mcvarious

Likelihood formulation configurations.

bounds_funccallable, optional

Advanced hook for expressing a joint/simplex constraint between the two currently-fixed test parameters and the free nuisance parameters (see module docstring). If provided, called as bounds_func({fix_A: test_vA, fix_B: test_vB}, free_indices, bounds_list) and must return one (lo, hi) tuple per entry of free_indices, in that order. If None (default), each free parameter uses its static bounds_list[i] unmodified – identical to the behavior before bounds_func existed.

constraintslist of scipy.optimize.LinearConstraint/NonlinearConstraint, optional

Joint constraints among the full n_params vector, evaluated directly on the fully-assembled parameter point (no projection needed). If None/empty (default), behavior is unchanged.

Returns:

min_nllfloat

The profiled minimum negative log-likelihood.

best_parray_like, 1D

The full parameter array (size N) including the two fixed values and the optimized nuisance parameters.

pyfc.optimizers.unconditional_fit_scipy(data, n_params, bounds_list, compute_rates_func, seed=None, pdf_components=None, likelihood_type='binned', S_sumw2=None, B_sumw2=None, use_finite_mc=False, bounds_func=None, constraints=None, scipy_method=None, n_restarts=1)

Performs an unconditional global maximum likelihood fit using SciPy’s L-BFGS-B.

Statistical Theory & Assumptions: To evaluate the denominator of the Profile Likelihood Ratio (PLR), we must find the parameter combination that yields the absolute minimum Negative Log-Likelihood (NLL).

Mathematical Expression: NLL_uncond = min_{theta} NLL(theta | data) where theta represents the full set of N free parameters in the model.

The L-BFGS-B algorithm is a quasi-Newton method that approximates the Broyden-Fletcher-Goldfarb-Shanno (BFGS) update to the Hessian matrix, scaled for bounded parameter spaces. It assumes the likelihood surface is sufficiently smooth and twice-differentiable. Because it is a local optimizer, it can be sensitive to the initial seed.

Parameters:

dataarray_like

The observed dataset (binned counts of any shape, or unbinned events).

n_paramsint

Total number of parameters in the model.

bounds_listlist of tuples

A list of (min, max) bounds for each parameter.

compute_rates_funccallable

User-provided mapping function matching parameters to physical expectations.

seedarray_like, 1D, optional

Initial parameter guess. If None, the midpoint of the bounds is used.

pdf_componentslist of callable, optional

Probability density functions, one per model component (e.g. signal, background). Required when likelihood_type=”unbinned”; each is evaluated exactly once on data and the resulting arrays are reused across every NLL evaluation in this fit. Unused for “binned”.

likelihood_typestr

Either “binned” or “unbinned”.

S_sumw2, B_sumw2array_like, optional

Sum of squared MC weights (sumw2) per bin, for finite MC binned likelihoods.

use_finite_mcbool

Flag to enable the Poisson-Gamma mixture likelihood for binned data.

bounds_funccallable, optional

Advanced hook for joint/simplex parameter constraints (see module docstring). If provided, called as bounds_func({}, list(range(n_params)), bounds_list) – an empty fixed_values dict since nothing is fixed in an unconditional fit – and must return one (lo, hi) tuple per parameter, in index order. If None (default), bounds_list is used unmodified.

constraintslist of scipy.optimize.LinearConstraint/NonlinearConstraint, optional

Joint constraints among the full n_params vector (see module docstring). When non-empty, method defaults to ‘SLSQP’ instead of ‘L-BFGS-B’ (required, since L-BFGS-B does not support constraints). When None/empty (default), behavior and method are UNCHANGED from before this parameter existed.

scipy_methodstr, optional

Explicit override for the scipy.optimize.minimize method (e.g. ‘trust-constr’ for better-conditioned but slower constrained fits). Takes precedence over the constraints-based default.

n_restartsint, optional

Number of distinct starting points to try, keeping whichever converges to the lowest NLL (see _minimize_with_restarts). Default 1 preserves the original single-fit behavior (a single fit from x0), except that res.success is now always checked; a non-converged result triggers one cheap perturbed retry and a warning if that also fails.

Returns:

min_nllfloat

The minimal negative log-likelihood found by the optimizer.

best_paramsarray_like, 1D

The parameter combination corresponding to min_nll.

pyfc.optimizers.unconditional_fit_ultranest(data, n_params, bounds_list, compute_rates_func, verbose=1, pdf_components=None, likelihood_type='binned', S_sumw2=None, B_sumw2=None, use_finite_mc=False, bounds_func=None, constraints=None)

Performs an unconditional global maximum likelihood fit using UltraNest.

Statistical Theory & Assumptions: UltraNest relies on Reactive Nested Sampling, an algorithm traditionally used for Bayesian evidence computation and posterior sampling. However, because it exhaustively maps the likelihood surface across the entire prior volume without relying on gradients, it acts as a highly robust global optimizer, virtually immune to local minima traps.

The algorithm evaluates the likelihood by converting a uniformly sampled unit hypercube [0, 1]^N into the physical parameter bounds via a prior transform. The returned value is the Maximum Likelihood Estimate (MLE) discovered during the integration phase.

Parameters:

dataarray_like

The observed dataset (binned counts of any shape, or unbinned events).

n_paramsint

Total number of parameters in the model.

bounds_listlist of tuples

A list of (min, max) bounds for each parameter.

compute_rates_funccallable

User-provided mapping function matching parameters to physical expectations.

verboseint

Verbosity level controlling UltraNest output (2 for full tracking).

pdf_componentslist of callable, optional

Probability density functions, one per model component. Required when likelihood_type=”unbinned”; unused for “binned”.

likelihood_type, S_sumw2, B_sumw2, use_finite_mcvarious

Likelihood formulation configurations.

bounds_funccallable, optional

Advanced hook for joint/simplex parameter constraints (see module docstring). If provided, called as bounds_func({}, list(range(n_params)), bounds_list) – an empty fixed_values dict since nothing is fixed in an unconditional fit – and must return one (lo, hi) tuple per parameter, in index order. If None (default), bounds_list is used unmodified.

constraintslist of scipy.optimize.LinearConstraint/NonlinearConstraint, optional

Joint constraints among the full n_params vector. Since nested sampling doesn’t use gradients, it doesn’t suffer the flat-gradient trap that motivates the scipy-side SLSQP/projection machinery: a violated constraint simply makes log_likelihood return a very large negative value (-1e10) at that point, no projection needed. If None/empty (default), behavior is unchanged.

Returns:

min_nllfloat

The absolute minimum negative log-likelihood (-1 * max_logL).

best_paramsarray_like, 1D

The parameter combination corresponding to the maximum log-likelihood point.

pyfc.orchestrator module

Feldman-Cousins Orchestrator Module

This file acts as the central hub and main driver for the PyFC package. It coordinates the execution of the Feldman-Cousins confidence interval construction, managing the interplay between the data, physical models, likelihood optimizers, and Monte Carlo toy generators. It handles checkpointing, parameter space scanning (including an optimized 2D grid sparsification algorithm), and data archiving.

Usage (Command Line Interface):

This script can be executed directly from the terminal to run a demonstration of the PyFC pipeline using automatically generated mock data (binned or unbinned).

To run the demonstration using the current configuration (or fc_config.json):

$ python orchestrator.py

Alternatively, if the PyFC package is installed in your Python environment:

$ python -m pyfc.orchestrator

You can override configuration parameters directly via command-line arguments. For example, to run an unbinned analysis with 500 toys using the SciPy optimizer:

$ python orchestrator.py –likelihood_type unbinned –n_toys 500 –strategy scipy

Created: v0.1.0 (July 24, 2026) Author: Mauricio Bustamante (mbustamante@gmail.com)

This file was released as part of the PyFC code, stored at https://github.com/mbustama/FeldmanCousins, which exists under a GNU GPL v3 License.

class pyfc.orchestrator.NumpyEncoder(*, skipkeys=False, ensure_ascii=True, check_circular=True, allow_nan=True, sort_keys=False, indent=None, separators=None, default=None)

Bases: JSONEncoder

Custom JSON encoder to safely serialize NumPy data types.

Standard json libraries cannot handle np.ndarray or specific NumPy numerical types (like np.float64). This intercepts those objects during the json.dump process and converts them to standard Python equivalents.

default(obj)

Implement this method in a subclass such that it returns a serializable object for o, or calls the base implementation (to raise a TypeError).

For example, to support arbitrary iterators, you could implement default like this:

def default(self, o):
    try:
        iterable = iter(o)
    except TypeError:
        pass
    else:
        return list(iterable)
    # Let the base class default method raise the TypeError
    return JSONEncoder.default(self, o)

pyfc.orchestrator.compute_fc_intervals(data, grids, compute_rates_func=None, generate_toy_func=None, pdf_components=None, cl=None, n_toys=500, strategy='scipy', num_cores=None, verbose=1, adaptive_toys=True, toy_batch_size=200, sparsify_grid=False, warm_start=True, likelihood_type='binned', S_mc_pool=None, B_mc_pool=None, output_file=None, save_log=True, save_directory='output/example_fc_output', use_finite_mc_correction_binned=True, S_sumw2=None, B_sumw2=None, compute_1D_intervals=True, compute_2D_intervals=True, param_names=None, smooth_1d=True, smooth_2d=True, bounds_func=None, constraints=None, scipy_method=None, n_restarts=1, neighbor_seeding=True)

Main execution pipeline for the Feldman-Cousins unified approach.

Statistical Theory: The Feldman-Cousins method constructs confidence intervals (or regions) that transition smoothly between upper limits and two-sided bounds. It relies on the Profile Likelihood Ratio as the test statistic to determine ordering.

For a given parameter point $ heta_{ ext{test}}$, the data test statistic is: $$ t_{ ext{data}} = -2 ln

rac{mathcal{L}( heta_{ ext{test}}, hat{hat{oldsymbol{ u}}} | ext{data})}{mathcal{L}(hat{oldsymbol{ heta}} | ext{data})} $$

where $hat{oldsymbol{ heta}}$ are the unconditional Maximum Likelihood Estimate (MLE) parameters, and $hat{hat{oldsymbol{

u}}}$ are the nuisance parameters profiled (maximized)

while fixing $ heta = heta_{ ext{test}}$. (Note: Since calc_nll returns $-2lnmathcal{L}$, this simplifies in code to cond_nll - data_uncond_nll).

Because the asymptotic distribution of $t$ (Wilks’ theorem) may fail near physical boundaries, we evaluate the exact critical threshold $t_{ ext{critical}}$ at a required confidence level $lpha$ by generating and fitting Monte Carlo pseudo-experiments (toys) generated under the null hypothesis $( heta_{ ext{test}}, hat{hat{oldsymbol{

u}}})$.

The point $ heta_{ ext{test}}$ is included in the final confidence set if: $$ t_{ ext{data}} leq t_{ ext{critical}}(lpha) $$

dataarray_like

The observed data: binned counts (any N-dimensional shape, e.g. a genuine 2D (E, cos_theta) histogram) or unbinned events.

gridslist of np.ndarray

A list of arrays, each defining the evaluation points for a parameter.

compute_rates_funccallable

User-provided mapping function matching parameters to physical expectations. For likelihood_type=”binned”: (params, S_sumw2, B_sumw2) -> (mu, sigma2), where mu/sigma2 must have the same shape as data; any fixed template arrays should be referenced via closure/module-global rather than passed in. For likelihood_type=”unbinned”: (params, probs) -> (expected_total, p_events), where probs is the list of pre-evaluated pdf_components densities.

generate_toy_funccallable, optional

User-provided unbinned parametric bootstrap function.

pdf_componentslist of callable, optional

Probability density functions, one per model component (e.g. signal, background, …no longer limited to exactly two). Required for likelihood_type=”unbinned”; ignored (with a warning) for “binned”. Each is evaluated exactly once per fit as pdf(data) and the resulting arrays are reused across every subsequent NLL evaluation in that fit.

cllist of float, optional

Confidence levels to compute (e.g., [0.68, 0.90]).

n_toysint, optional

Number of pseudo-experiments to generate per point.

strategystr, optional

Optimizer strategy: “grid”, “scipy”, “ultranest”, or “hybrid”.

num_coresint, optional

Number of parallel threads for toy generation. Defaults to None, which is forwarded as-is to the underlying ThreadPoolExecutor/ProcessPoolExecutor/numba.set_num_threads calls; all three treat None as “use every available hardware thread” (equivalent to os.cpu_count()). Pass an explicit integer to cap this, e.g. to match a Slurm/PBS core allocation.

verboseint, optional

0 = Silent, 1 = Normal, 2 = Debug.

adaptive_toysbool, optional

When True and strategy in (“scipy”, “ultranest”, “hybrid”), stops generating toys for a grid point early once its accept/reject verdict is statistically settled (a strict 99.9% Wilson confidence interval on the running p-value estimate lies entirely on one side of the target significance, and at least toys.ADAPTIVE_MIN_TOYS toys have been generated) – see toys._adaptive_stop_decision. This only meaningfully saves time for points far from the accept/ reject boundary; points near it will and should run the full n_toys. No effect for strategy=”grid” (its toy generators are a numba prange-parallelized batch operation, not a natural fit for a sequential early-exit check).

toy_batch_sizeint, optional

For strategy in (“scipy”, “ultranest”, “hybrid”), toys are submitted/collected from the executor in batches of this size instead of all n_toys at once, bounding how many toy datasets / in-flight results are alive at once (matters most for likelihood_type=”unbinned”; binned toys are small fixed-size per-bin arrays with no comparable memory concern). Also the granularity at which adaptive_toys checks its stopping condition. Same total n_toys are generated either way when adaptive_toys=False – this only changes memory/dispatch pattern, never the statistics. No effect for strategy=”grid”.

warm_startbool

Algorithmic enhancement to reduce execution time via checkpointing.

sparsify_gridbool, optional

For 2D scans, coarsely samples the grid and interpolates the t_critical surface, then refines only cells adjacent to the interpolated accept/reject boundary, instead of evaluating every cell. Defaults to False. Known limitation: refinement is a single, non-iterative pass with a 1-cell-wide halo around ~5-per-axis coarse nodes, so for grids much larger than ~20x20 per axis it under-counts the true accepted region (deep interior/exterior cells far from any coarse node are never evaluated and are force-excluded as a result). See the comment above the “Sparsification Array Tracing” block in this function’s body for details.

likelihood_typestr

“binned” or “unbinned”.

S_mc_pool, B_mc_poolarray_like, optional

Pool of events to bootstrap for unbinned toy generation.

output_file, save_log, save_directorystr/bool

I/O file structures. save_directory (default “output/example_fc_output”, relative to the current working directory) is where checkpoints and final results are written; see the README’s “Outputs, Plots, and Checkpointing” section. output_file defaults to None, meaning no final .npz/.json archive is written – only the in-memory results dict (this function’s first return value) is populated; pass a string prefix (e.g. “fc_results”) to also write {save_directory}/{output_file}.npz and .json.

use_finite_mc_correction_binnedbool, optional

Toggle for Poisson-Gamma mixture likelihood.

S_sumw2, B_sumw2array_like, optional

Sum of squared MC weights (sumw2) per bin, for the finite MC correction.

compute_1D_intervals, compute_2D_intervalsbool

Switches for calculating 1D profiles or 2D joint contours.

param_nameslist of str, optional

Names used for plot labels. Defaults to None, which auto-generates [“param1”, “param2”, …] matching len(grids) – i.e. it always matches the actual number of parameters in this specific model, unlike hardcoding a fixed-length example list would.

smooth_1d, smooth_2dbool, optional

Toggles interpolation smoothing for final plot outputs.

bounds_funccallable, optional

Advanced hook that lets the box bounds of the profiled (free) parameters depend on whatever parameter(s) are currently fixed by the scan, e.g. to express a joint constraint like f_e + f_mu <= 1 without a user-side penalty function. Forwarded to whichever conditional_fit_*/unconditional_fit_* optimizer is active (scipy or ultranest strategies only; ignored for strategy=”grid”). See the optimizers.py module docstring for the full signature contract and a worked example, and the README section “Handling joint/simplex- constrained parameters” for an end-to-end walkthrough. Limitation: this only handles constraints between a currently-FIXED test parameter and a free nuisance parameter – it cannot express a constraint between two parameters that are BOTH free at the same time (use constraints for that instead).

constraintslist of scipy.optimize.LinearConstraint/NonlinearConstraint, optional

Joint constraints among the FULL n_params vector, for expressing constraints between multiple simultaneously-free nuisance parameters (which bounds_func cannot express). Forwarded to whichever scipy/ ultranest optimizer is active; ignored for strategy=”grid”. For the scipy path, method switches from ‘L-BFGS-B’ to ‘SLSQP’ whenever constraints are non-empty (required, since L-BFGS-B doesn’t support constraints); for the ultranest path, a violated constraint simply yields a very large negative log-likelihood at that point. When None/empty (default), behavior and method are UNCHANGED from before this parameter existed. See the optimizers.py module docstring and the README section “Handling joint/simplex-constrained parameters”.

scipy_methodstr, optional

Explicit override for the scipy.optimize.minimize method (e.g. ‘trust-constr’ for better-conditioned but slower constrained fits). Takes precedence over the constraints-based default.

n_restartsint, optional

Forwarded to the scipy fit functions for the DATA fits (Phase 0’s unconditional fit and the Phase 1/2 conditional fits) – NOT the MC toy fits, which are already well-seeded from the conditional MLE (true_params). Number of distinct starting points to try per fit, keeping whichever converges to the lowest NLL (see optimizers._minimize_with_restarts). Default 1 preserves pre- FIX-4 behavior (a single fit per point), except that res.success is now always checked, with a cheap perturbed retry and warning on non-convergence – this applies even at the default n_restarts=1, for every scipy fit call including the toy fits.

neighbor_seedingbool, optional

When True (default) and strategy=”scipy”, the DATA fit (not the MC toys, which are already seeded from the conditional MLE) at each 1D/ 2D scan grid point is seeded from an adjacent, already-evaluated grid point’s profiled parameters, instead of always starting fresh from the bounds midpoint – these loops already iterate the grid in order, so neighboring points’ solutions are informative starting guesses. Whenever a neighbor seed is used, the effective number of restarts for that specific call is raised to at least 2 (the neighbor-seeded start plus a fresh bounds-midpoint start, keeping whichever is better), so a bad neighbor optimum can’t silently cascade forward across many subsequent grid points. Set False to disable and always start from the bounds midpoint (the original behavior for the data fit).

resultsdict

A massive dictionary containing all t_data, thresholds, accepted masks, and best-fit geometries evaluated across the parameter space.

figmatplotlib.figure.Figure or None

The Matplotlib Figure object generated by generate_corner_plot, or None if plotting was skipped.

pyfc.orchestrator.tqdm(iterable, *args, **kwargs)

pyfc.plotting module

Plotting and Visualization Module

This module provides the visualization routines for the PyFC package, generating publication-quality corner plots that map the multi-dimensional parameter space. It visualizes both the 1D Profile Likelihood Ratios (with integrated Monte Carlo thresholds) and the 2D joint confidence contours using matplotlib.

Created: v0.1.0 (July 24, 2026) Author: Mauricio Bustamante (mbustamante@gmail.com)

This file was released as part of the PyFC code, stored at https://github.com/mbustama/FeldmanCousins, which exists under a GNU GPL v3 License.

pyfc.plotting.generate_corner_plot(results, config)

Generates a lower-triangle corner plot of the Feldman-Cousins confidence intervals.

Statistical & Visual Context: The diagonal axes display the 1D Profile Likelihood Ratio (PLR). The solid black curve represents the data test statistic ($t_{ ext{data}}$), while the dashed colored lines represent the exact critical thresholds ($t_{ ext{critical}}$) derived from Monte Carlo toys at specified Confidence Levels (CL). The shaded regions denote the accepted parameter intervals where $t_{ ext{data}} leq t_{ ext{critical}}$.

The off-diagonal axes display the 2D joint contours. The contours are drawn exactly where the differential surface $z = t_{ ext{data}} - t_{ ext{critical}}$ crosses zero, providing a smooth boundary of the rigorous frequentist confidence region.

Parameters:

resultsdict

The comprehensive results dictionary produced by compute_fc_intervals, containing evaluated test statistics, parameter grids, and critical thresholds.

configdict

Visualization configuration dictionary containing: - n_params (int): Number of physical parameters. - param_names (list of str): Axis labels. - cl (list of float): Confidence levels to plot. - smooth_1d, smooth_2d (bool): Toggles for Gaussian interpolation smoothing. - save_directory (str): Output path.

Returns:

figmatplotlib.figure.Figure, or None

The generated figure, also saved to disk as fc_corner_plot.pdf. The caller owns the returned figure and is responsible for closing it (e.g. plt.close(fig)) if generating many plots in a loop – this function does not close it itself. Returns None instead if matplotlib isn’t installed (MATPLOTLIB_AVAILABLE is False).

pyfc.toys module

Monte Carlo Toy Generation and Fitting Module

This module manages the core frequentist computation of the PyFC package: the generation and parallel evaluation of pseudo-experiments (toys). To establish exact coverage in the Feldman-Cousins framework, the distribution of the test statistic must be empirically derived at each point in the parameter grid by simulating data under the null hypothesis, and fitting it both unconditionally and conditionally.

This file specifically handles the execution when using continuous optimizers (SciPy, UltraNest) and manages efficient parallelization by routing workloads to either Thread or Process pools depending on Python’s Global Interpreter Lock (GIL) limitations for the requested likelihood type.

Created: v0.1.0 (July 24, 2026) Author: Mauricio Bustamante (mbustamante@gmail.com)

This file was released as part of the PyFC code, stored at https://github.com/mbustama/FeldmanCousins, which exists under a GNU GPL v3 License.

pyfc.toys.generate_and_fit_toys_python(true_params, n_params, fit_mode, fix_idx, fix_A, fix_B, t_vA, t_vB, bounds_list, n_toys, strategy, num_cores=None, verbose=1, pdf_components=None, likelihood_type='binned', S_mc_pool=None, B_mc_pool=None, S_sumw2=None, B_sumw2=None, use_finite_mc=False, compute_rates_func=None, generate_toy_func=None, bounds_func=None, constraints=None, scipy_method=None, toy_batch_size=None, adaptive_toys=False, t_data=None, alpha=None)

Handles threaded generation and fitting of MC toys for 1D profiling and 2D contours.

Statistical Context: To evaluate if a grid point $t_{test}$ belongs in the confidence interval, we generate n_toys pseudo-experiments assuming the null hypothesis ($ heta = t_{test}$, with nuisance parameters at their conditional maximum likelihood values true_params). Each toy is then fit unconditionally and conditionally to build the empirical distribution of the Profile Likelihood Ratio.

Parallelization Context: - Unbinned PDFs often rely on pure Python functions (like scipy.stats.norm.pdf).

These are strictly bound by the Global Interpreter Lock (GIL). Multithreading provides zero speedup. Thus, we branch to a ProcessPoolExecutor to spawn distinct processes.

• Binned operations (vectorized NumPy math) naturally release the GIL. Thus, we can use a lighter-weight ThreadPoolExecutor and avoid the heavy overhead of process spawning and memory IPC.

Parameters:

true_paramsarray_like

The physical parameters (test values + profiled nuisance values) used to generate toys.

n_paramsint

Total number of parameters in the model.

fit_modestr

Either “1d” or “2d”.

fix_idx, fix_A, fix_Bint or None

Indices of the fixed parameters depending on the fit mode.

t_vA, t_vBfloat

The test values to fix the parameters at during the conditional fit.

bounds_listvarious

Boundary constraints.

n_toysint

Number of pseudo-experiments to generate and evaluate.

strategystr

Optimizer selection (“scipy”, “ultranest”, “hybrid”).

num_coresint

Number of parallel workers.

verboseint

Logging level.

pdf_componentslist of callable, optional

Probability density functions, one per model component. Required when likelihood_type=”unbinned”; unused for “binned”.

likelihood_typestr

“binned” or “unbinned”.

S_mc_pool, B_mc_poolarray_like, optional

Source pools for bootstrapping unbinned events.

S_sumw2, B_sumw2array_like, optional

Sum of squared MC weights (sumw2) per bin, for finite MC corrections.

use_finite_mcbool

Toggle for Poisson-Gamma mixture likelihoods.

compute_rates_funccallable

User-provided expected rates mapping function.

generate_toy_funccallable

User-provided unbinned toy bootstrapper.

bounds_funccallable, optional

Advanced hook for joint/simplex parameter constraints, forwarded unchanged to the underlying *_scipy/*_ultranest fit calls (see optimizers.py module docstring). If None (default), behavior is unchanged from before this parameter existed.

constraintslist of scipy.optimize.LinearConstraint/NonlinearConstraint, optional

Joint constraints, forwarded unchanged to the underlying *_scipy/*_ultranest fit calls. If None/empty (default), behavior is unchanged.

scipy_methodstr, optional

Explicit scipy.optimize.minimize method override, forwarded to the scipy fit calls.

toy_batch_sizeint, optional

Chunk size for dispatching toys to the executor: instead of submitting all n_toys tasks in one executor.map call, they are submitted/collected in batches of this size, bounding how many toy datasets / in-flight result objects are alive at once (matters most for likelihood_type=”unbinned”, where each toy is a variable-size event array passed through ProcessPoolExecutor IPC; binned toys are small fixed-size per-bin arrays and don’t have this memory concern). None or <= 0 means “one batch of the full n_toys” (previous behavior, unchanged numerics either way – this only affects how toys are grouped for dispatch/early-stop checks, never how many are generated when adaptive_toys=False).

adaptive_toysbool, optional

If True, stop generating toys for this point early once the accept/reject verdict at alpha is statistically settled (see _adaptive_stop_decision). Checked once per toy_batch_size batch, so a smaller toy_batch_size allows finer-grained (but not free – still bounded by ADAPTIVE_MIN_TOYS) early stopping. Requires t_data and alpha to be provided; silently has no effect otherwise (falls back to always generating all n_toys).

t_datafloat, optional

The observed data’s test statistic at this grid point, required for adaptive_toys to have any effect.

alphafloat, optional

Target significance (1 - CL) for the adaptive_toys stopping decision. When multiple confidence levels are requested, callers should pass the smallest alpha (largest CL) – see _adaptive_stop_decision’s docstring.

Returns:

t_statsnp.ndarray

Array of the computed test statistics – length n_toys, unless adaptive_toys stopped early, in which case it is shorter. Callers must index/quantile using len(t_stats), not the original n_toys.

pyfc.unbinned module

Unbinned Likelihood and Grid Optimization Module

This file implements the Extended Unbinned Maximum Likelihood (EUML) routines for the PyFC package. Unlike binned analyses which rely on histograms, the unbinned method evaluates the probability density directly on an event-by-event basis. This is statistically powerful for small datasets where binning would cause unacceptable information loss.

The module also includes brute-force grid search optimizers for the unbinned likelihood, which serve as robust alternatives to gradient-based minimizers when the likelihood surface is highly irregular.

Created: v0.1.0 (July 24, 2026) Author: Mauricio Bustamante (mbustamante@gmail.com)

This file was released as part of the PyFC code, stored at https://github.com/mbustama/FeldmanCousins, which exists under a GNU GPL v3 License.

pyfc.unbinned.calc_nll_unbinned(params, len_obs, probs, compute_rates_func)

Computes the Extended Unbinned Negative Log-Likelihood (NLL).

Mathematical Expression: -ln(L) = N_expected - sum(ln(lambda(x_i)))

Note: The factorial term ln(N_obs!) is dropped from the standard EUML formula because it is a constant that depends only on the data, and therefore strictly cancels out when computing delta-NLL test statistics.

Parameters:

paramsarray_like

The physical parameters.

len_obsint

The total number of observed events (N_obs).

probslist of np.ndarray

The pre-evaluated PDFs (one per pdf_components[i]) for the observed events.

compute_rates_funccallable

User-provided function (params, probs) -> (expected_total, p_events) mapping parameters and the pre-evaluated component densities to the expected total event count and the per-event mixture density.

Returns:

float

The calculated NLL value. For events where the model predicts an unphysical or near-zero density (p_events[k] <= 1e-12), that event’s contribution is a smooth, continuous extension of -log(p) rather than a flat penalty (see note below).

Note on unphysical (p_events <= p_floor = 1e-12) events:

Earlier versions of this function returned a flat penalty (1e10) the instant ANY single event had a non-positive predicted density, short- circuiting the entire likelihood. Because scipy.optimize.minimize(…, method=’L-BFGS-B’) estimates gradients via finite differences, a locally -constant NLL surface reads as “gradient ~ 0, already converged,” which can permanently trap the optimizer. Instead, each offending event’s density is clipped to a tiny positive floor for the -log(p) term (a continuous extension of the real likelihood), and a quadratic barrier proportional to how far below that floor p_events[k] actually is is added per-event, so the gradient stays informative instead of flat. This is applied elementwise (not collapsed via np.any), so multiple unphysical events each contribute independently. The trigger is p_events <= p_floor rather than p_events <= 0 specifically so there is no razor-thin discontinuity in the tiny sliver (0, p_floor): the real -log(p) term diverges as p_events -> 0+, while the floor-based term is constant (evaluated at p_floor, not at the actual p_events[k]), so triggering only at p_events <= 0 would make the NLL drop right at that boundary instead of continuing to climb. For all p_events[k] > p_floor (i.e. all physically meaningful values – p_floor is astronomically small), behavior is practically identical to before this change (matches to within 1e-12 relative/absolute tolerance; see tests/test_smoothing.py).

pyfc.unbinned.conditional_fit_grid_unbinned_1d(test_val, fix_idx, n_params, obs_events, pdf_components, cond_grid_points, compute_rates_func)

Performs a brute-force 1D conditional profile scan for the unbinned likelihood. Fixes one parameter of interest and scans the remaining grid.

Parameters:

test_valfloat

The fixed value for the parameter of interest.

fix_idxint

The index of the parameter to fix.

n_paramsint

Total number of parameters.

obs_eventsarray_like

The observed data events.

pdf_componentslist of callable

Probability density functions, one per model component.

cond_grid_pointsarray_like, 2D

The sub-grid of nuisance parameters to profile over.

compute_rates_funccallable

User-provided expected rates function.

Returns:

min_nllfloat

The minimum negative log-likelihood given the fixed parameter.

best_paramsarray_like, 1D

The parameter combination that produced min_nll.

pyfc.unbinned.conditional_fit_grid_unbinned_2d(test_vA, test_vB, fix_A, fix_B, n_params, obs_events, pdf_components, cond_grid_points, compute_rates_func)

Performs a brute-force 2D conditional profile scan for the unbinned likelihood. Fixes two parameters of interest to map joint contours.

Parameters:

test_vA, test_vBfloat

The fixed values for the two parameters of interest.

fix_A, fix_Bint

The indices of the fixed parameters.

n_paramsint

Total number of parameters.

obs_eventsarray_like

The observed data events.

pdf_componentslist of callable

Probability density functions, one per model component.

cond_grid_pointsarray_like, 2D

The sub-grid of nuisance parameters to profile over.

compute_rates_funccallable

User-provided expected rates function.

Returns:

min_nllfloat

The minimum negative log-likelihood given the two fixed parameters.

best_paramsarray_like, 1D

The parameter combination that produced min_nll.

pyfc.unbinned.generate_and_fit_toys_grid_unbinned_1d(test_val, fix_idx, true_params, n_params, pdf_components, full_grid_points, cond_grid_points, n_toys, S_mc_pool, B_mc_pool, compute_rates_func, generate_toy_func)

Sequentially generates and evaluates 1D unbinned toys entirely using grid search.

Used when gradient optimizers fail and external parallelization is disabled. It evaluates the Profile Likelihood Ratio (cond_nll - uncond_nll) for each simulated unbinned dataset.

Parameters:

test_valfloat

The fixed value for the parameter of interest.

fix_idxint

The index of the fixed parameter.

true_paramsarray_like, 1D

The physical parameters dictating the expected event yields.

n_paramsint

Total number of parameters in the model.

pdf_componentslist of callable

Probability density functions, one per model component.

full_grid_pointsarray_like, 2D

Grid for the unconditional fit.

cond_grid_pointsarray_like, 2D

Grid for the conditional (nuisance parameter) fit.

n_toysint

Number of pseudo-experiments to generate.

S_mc_pool, B_mc_poolnp.ndarray

Pre-generated continuous events for bootstrapping.

compute_rates_funccallable

User-provided physics mapping function.

generate_toy_funccallable

User-provided function to bootstrap the toy data.

Returns:

t_statisticsarray_like, 1D

Array of length n_toys containing the test statistic for each toy.

pyfc.unbinned.generate_and_fit_toys_grid_unbinned_2d(test_vA, test_vB, fix_A, fix_B, true_params, n_params, pdf_components, full_grid_points, cond_grid_points, n_toys, S_mc_pool, B_mc_pool, compute_rates_func, generate_toy_func)

Sequentially generates and evaluates 2D unbinned toys entirely using grid search.

Parameters:

test_vA, test_vBfloat

The fixed values for the two parameters of interest.

fix_A, fix_Bint

The indices of the parameters being fixed.

true_paramsarray_like, 1D

The physical parameters dictating the expected event yields.

n_paramsint

Total number of parameters in the model.

pdf_componentslist of callable

Probability density functions, one per model component.

full_grid_pointsarray_like, 2D

Grid for the unconditional fit.

cond_grid_pointsarray_like, 2D

Grid for the conditional (nuisance parameter) fit.

n_toysint

Number of pseudo-experiments to generate.

S_mc_pool, B_mc_poolnp.ndarray

Pre-generated continuous events for bootstrapping.

compute_rates_funccallable

User-provided physics mapping function.

generate_toy_funccallable

User-provided function to bootstrap the toy data.

Returns:

t_statisticsarray_like, 1D

Array of length n_toys containing the test statistic for each toy.

pyfc.unbinned.unconditional_fit_grid_unbinned(obs_events, pdf_components, full_grid_points, compute_rates_func)

Performs a brute-force global scan over the parameter space to find the unconditional MLE.

Because evaluating a PDF function pdf(obs_events) for every parameter combination can be incredibly slow in Python, this function evaluates each PDF exactly ONCE for the observed dataset, and reuses those static arrays across the entire grid loop.

Parameters:

obs_eventsarray_like

The observed data events.

pdf_componentslist of callable

Probability density functions, one per model component (e.g. signal, background).

full_grid_pointsarray_like, 2D

Matrix containing all parameter combinations to scan.

compute_rates_funccallable

User-provided expected rates function.

Returns:

min_nllfloat

The absolute minimum negative log-likelihood found.

best_paramsarray_like, 1D

The parameter combination corresponding to min_nll.

Module contents

Feldman-Cousins Frequentist Analysis Framework

Created: v0.1.0 (July 24, 2026) Author: Mauricio Bustamante (mbustamante@gmail.com)

pyfc.compute_fc_intervals(data, grids, compute_rates_func=None, generate_toy_func=None, pdf_components=None, cl=None, n_toys=500, strategy='scipy', num_cores=None, verbose=1, adaptive_toys=True, toy_batch_size=200, sparsify_grid=False, warm_start=True, likelihood_type='binned', S_mc_pool=None, B_mc_pool=None, output_file=None, save_log=True, save_directory='output/example_fc_output', use_finite_mc_correction_binned=True, S_sumw2=None, B_sumw2=None, compute_1D_intervals=True, compute_2D_intervals=True, param_names=None, smooth_1d=True, smooth_2d=True, bounds_func=None, constraints=None, scipy_method=None, n_restarts=1, neighbor_seeding=True)

Main execution pipeline for the Feldman-Cousins unified approach.

Statistical Theory: The Feldman-Cousins method constructs confidence intervals (or regions) that transition smoothly between upper limits and two-sided bounds. It relies on the Profile Likelihood Ratio as the test statistic to determine ordering.

For a given parameter point $ heta_{ ext{test}}$, the data test statistic is: $$ t_{ ext{data}} = -2 ln

rac{mathcal{L}( heta_{ ext{test}}, hat{hat{oldsymbol{ u}}} | ext{data})}{mathcal{L}(hat{oldsymbol{ heta}} | ext{data})} $$

where $hat{oldsymbol{ heta}}$ are the unconditional Maximum Likelihood Estimate (MLE) parameters, and $hat{hat{oldsymbol{

u}}}$ are the nuisance parameters profiled (maximized)

while fixing $ heta = heta_{ ext{test}}$. (Note: Since calc_nll returns $-2lnmathcal{L}$, this simplifies in code to cond_nll - data_uncond_nll).

Because the asymptotic distribution of $t$ (Wilks’ theorem) may fail near physical boundaries, we evaluate the exact critical threshold $t_{ ext{critical}}$ at a required confidence level $lpha$ by generating and fitting Monte Carlo pseudo-experiments (toys) generated under the null hypothesis $( heta_{ ext{test}}, hat{hat{oldsymbol{

u}}})$.

The point $ heta_{ ext{test}}$ is included in the final confidence set if: $$ t_{ ext{data}} leq t_{ ext{critical}}(lpha) $$

dataarray_like

The observed data: binned counts (any N-dimensional shape, e.g. a genuine 2D (E, cos_theta) histogram) or unbinned events.

gridslist of np.ndarray

A list of arrays, each defining the evaluation points for a parameter.

compute_rates_funccallable

User-provided mapping function matching parameters to physical expectations. For likelihood_type=”binned”: (params, S_sumw2, B_sumw2) -> (mu, sigma2), where mu/sigma2 must have the same shape as data; any fixed template arrays should be referenced via closure/module-global rather than passed in. For likelihood_type=”unbinned”: (params, probs) -> (expected_total, p_events), where probs is the list of pre-evaluated pdf_components densities.

generate_toy_funccallable, optional

User-provided unbinned parametric bootstrap function.

pdf_componentslist of callable, optional

Probability density functions, one per model component (e.g. signal, background, …no longer limited to exactly two). Required for likelihood_type=”unbinned”; ignored (with a warning) for “binned”. Each is evaluated exactly once per fit as pdf(data) and the resulting arrays are reused across every subsequent NLL evaluation in that fit.

cllist of float, optional

Confidence levels to compute (e.g., [0.68, 0.90]).

n_toysint, optional

Number of pseudo-experiments to generate per point.

strategystr, optional

Optimizer strategy: “grid”, “scipy”, “ultranest”, or “hybrid”.

num_coresint, optional

Number of parallel threads for toy generation. Defaults to None, which is forwarded as-is to the underlying ThreadPoolExecutor/ProcessPoolExecutor/numba.set_num_threads calls; all three treat None as “use every available hardware thread” (equivalent to os.cpu_count()). Pass an explicit integer to cap this, e.g. to match a Slurm/PBS core allocation.

verboseint, optional

0 = Silent, 1 = Normal, 2 = Debug.

adaptive_toysbool, optional

When True and strategy in (“scipy”, “ultranest”, “hybrid”), stops generating toys for a grid point early once its accept/reject verdict is statistically settled (a strict 99.9% Wilson confidence interval on the running p-value estimate lies entirely on one side of the target significance, and at least toys.ADAPTIVE_MIN_TOYS toys have been generated) – see toys._adaptive_stop_decision. This only meaningfully saves time for points far from the accept/ reject boundary; points near it will and should run the full n_toys. No effect for strategy=”grid” (its toy generators are a numba prange-parallelized batch operation, not a natural fit for a sequential early-exit check).

toy_batch_sizeint, optional

For strategy in (“scipy”, “ultranest”, “hybrid”), toys are submitted/collected from the executor in batches of this size instead of all n_toys at once, bounding how many toy datasets / in-flight results are alive at once (matters most for likelihood_type=”unbinned”; binned toys are small fixed-size per-bin arrays with no comparable memory concern). Also the granularity at which adaptive_toys checks its stopping condition. Same total n_toys are generated either way when adaptive_toys=False – this only changes memory/dispatch pattern, never the statistics. No effect for strategy=”grid”.

warm_startbool

Algorithmic enhancement to reduce execution time via checkpointing.

sparsify_gridbool, optional

For 2D scans, coarsely samples the grid and interpolates the t_critical surface, then refines only cells adjacent to the interpolated accept/reject boundary, instead of evaluating every cell. Defaults to False. Known limitation: refinement is a single, non-iterative pass with a 1-cell-wide halo around ~5-per-axis coarse nodes, so for grids much larger than ~20x20 per axis it under-counts the true accepted region (deep interior/exterior cells far from any coarse node are never evaluated and are force-excluded as a result). See the comment above the “Sparsification Array Tracing” block in this function’s body for details.

likelihood_typestr

“binned” or “unbinned”.

S_mc_pool, B_mc_poolarray_like, optional

Pool of events to bootstrap for unbinned toy generation.

output_file, save_log, save_directorystr/bool

I/O file structures. save_directory (default “output/example_fc_output”, relative to the current working directory) is where checkpoints and final results are written; see the README’s “Outputs, Plots, and Checkpointing” section. output_file defaults to None, meaning no final .npz/.json archive is written – only the in-memory results dict (this function’s first return value) is populated; pass a string prefix (e.g. “fc_results”) to also write {save_directory}/{output_file}.npz and .json.

use_finite_mc_correction_binnedbool, optional

Toggle for Poisson-Gamma mixture likelihood.

S_sumw2, B_sumw2array_like, optional

Sum of squared MC weights (sumw2) per bin, for the finite MC correction.

compute_1D_intervals, compute_2D_intervalsbool

Switches for calculating 1D profiles or 2D joint contours.

param_nameslist of str, optional

Names used for plot labels. Defaults to None, which auto-generates [“param1”, “param2”, …] matching len(grids) – i.e. it always matches the actual number of parameters in this specific model, unlike hardcoding a fixed-length example list would.

smooth_1d, smooth_2dbool, optional

Toggles interpolation smoothing for final plot outputs.

bounds_funccallable, optional

Advanced hook that lets the box bounds of the profiled (free) parameters depend on whatever parameter(s) are currently fixed by the scan, e.g. to express a joint constraint like f_e + f_mu <= 1 without a user-side penalty function. Forwarded to whichever conditional_fit_*/unconditional_fit_* optimizer is active (scipy or ultranest strategies only; ignored for strategy=”grid”). See the optimizers.py module docstring for the full signature contract and a worked example, and the README section “Handling joint/simplex- constrained parameters” for an end-to-end walkthrough. Limitation: this only handles constraints between a currently-FIXED test parameter and a free nuisance parameter – it cannot express a constraint between two parameters that are BOTH free at the same time (use constraints for that instead).

constraintslist of scipy.optimize.LinearConstraint/NonlinearConstraint, optional

Joint constraints among the FULL n_params vector, for expressing constraints between multiple simultaneously-free nuisance parameters (which bounds_func cannot express). Forwarded to whichever scipy/ ultranest optimizer is active; ignored for strategy=”grid”. For the scipy path, method switches from ‘L-BFGS-B’ to ‘SLSQP’ whenever constraints are non-empty (required, since L-BFGS-B doesn’t support constraints); for the ultranest path, a violated constraint simply yields a very large negative log-likelihood at that point. When None/empty (default), behavior and method are UNCHANGED from before this parameter existed. See the optimizers.py module docstring and the README section “Handling joint/simplex-constrained parameters”.

scipy_methodstr, optional

Explicit override for the scipy.optimize.minimize method (e.g. ‘trust-constr’ for better-conditioned but slower constrained fits). Takes precedence over the constraints-based default.

n_restartsint, optional

Forwarded to the scipy fit functions for the DATA fits (Phase 0’s unconditional fit and the Phase 1/2 conditional fits) – NOT the MC toy fits, which are already well-seeded from the conditional MLE (true_params). Number of distinct starting points to try per fit, keeping whichever converges to the lowest NLL (see optimizers._minimize_with_restarts). Default 1 preserves pre- FIX-4 behavior (a single fit per point), except that res.success is now always checked, with a cheap perturbed retry and warning on non-convergence – this applies even at the default n_restarts=1, for every scipy fit call including the toy fits.

neighbor_seedingbool, optional

When True (default) and strategy=”scipy”, the DATA fit (not the MC toys, which are already seeded from the conditional MLE) at each 1D/ 2D scan grid point is seeded from an adjacent, already-evaluated grid point’s profiled parameters, instead of always starting fresh from the bounds midpoint – these loops already iterate the grid in order, so neighboring points’ solutions are informative starting guesses. Whenever a neighbor seed is used, the effective number of restarts for that specific call is raised to at least 2 (the neighbor-seeded start plus a fresh bounds-midpoint start, keeping whichever is better), so a bad neighbor optimum can’t silently cascade forward across many subsequent grid points. Set False to disable and always start from the bounds midpoint (the original behavior for the data fit).

resultsdict

A massive dictionary containing all t_data, thresholds, accepted masks, and best-fit geometries evaluated across the parameter space.

figmatplotlib.figure.Figure or None

The Matplotlib Figure object generated by generate_corner_plot, or None if plotting was skipped.

pyfc.generate_config()

Main entry point for the configuration generator.

Statistical Theory & Configuration Impact: This sequence builds the fc_config.json matrix. It defines: - PDF formulation (Binned vs Unbinned) which changes the NLL definition. - Confidence Levels (CL) which sets the coverage integration targets. - Monte Carlo sampling size (n_toys), controlling the resolution of the empirical

test statistic PDF. Higher toys reduce statistical noise in the p-value calculation at the expense of computation time.

• Likelihood optimization strategies (SciPy, UltraNest, etc.) used to find the global and conditional minima of the parameter space during profiling.

• Advanced topological settings (adaptive toys, grid sparsification) which employ heuristics to skip unnecessary toy generations in regions definitively outside the targeted confidence bounds.

Parameters:

None

Returns:

None

Outputs a serialized JSON file containing the validated configuration.

pyfc.generate_corner_plot(results, config)

Generates a lower-triangle corner plot of the Feldman-Cousins confidence intervals.

Statistical & Visual Context: The diagonal axes display the 1D Profile Likelihood Ratio (PLR). The solid black curve represents the data test statistic ($t_{ ext{data}}$), while the dashed colored lines represent the exact critical thresholds ($t_{ ext{critical}}$) derived from Monte Carlo toys at specified Confidence Levels (CL). The shaded regions denote the accepted parameter intervals where $t_{ ext{data}} leq t_{ ext{critical}}$.

The off-diagonal axes display the 2D joint contours. The contours are drawn exactly where the differential surface $z = t_{ ext{data}} - t_{ ext{critical}}$ crosses zero, providing a smooth boundary of the rigorous frequentist confidence region.

Parameters:

resultsdict

The comprehensive results dictionary produced by compute_fc_intervals, containing evaluated test statistics, parameter grids, and critical thresholds.

configdict

Visualization configuration dictionary containing: - n_params (int): Number of physical parameters. - param_names (list of str): Axis labels. - cl (list of float): Confidence levels to plot. - smooth_1d, smooth_2d (bool): Toggles for Gaussian interpolation smoothing. - save_directory (str): Output path.

Returns:

figmatplotlib.figure.Figure, or None

The generated figure, also saved to disk as fc_corner_plot.pdf. The caller owns the returned figure and is responsible for closing it (e.g. plt.close(fig)) if generating many plots in a loop – this function does not close it itself. Returns None instead if matplotlib isn’t installed (MATPLOTLIB_AVAILABLE is False).