Likelihood-free inference with nuisance parameters through normalizing flows

Phil Assheton

Source-linked AI summary

Likelihood-free inference with nuisance parameters through normalizing flows

Phil Assheton

arXiv:2609.10534v1Sep 9, 2026stat.MEcs.LGstat.ML

TL;DR

Inference with nuisance parameters seeks pivotal statistics, but such statistics are difficult to construct generally. The paper decomposes a normalizing flow so a learned statistic depends on the interest parameter while nuisance effects are separated, and evaluates its calibration and power across classical problems. It nearly recovers the one-sample t-test, improves worst-case Welch-test size within a constrained variance-ratio range, and outperforms profile likelihood-ratio techniques in small-to-moderate-sample power at much lower computational cost.

Problem
Frequentist inference seeks tests whose size remains invariant across nuisance-parameter configurations, but general problems often lack analytically available pivotal statistics.
Method
The paper decomposes a normalizing flow into a pivoting component conditioned only on the interest parameter and a nuisance component, learning p-values from simulated data.
Results
The method nearly recovers the t-test, improves worst-case Welch-test size within constrained variance-ratio bounds, and has a 1.2 percentage-point average power advantage over bootstrap profile LR for partial biserial correlation.
Takeaways & Limitations
A single trained network can provide rapidly reusable p-values, with competitive calibration and power across the evaluated classical problems.
Takeaways & Limitations
The method is currently constrained to settings with consistent enough local quotient directions; globally, the learned statistic is a training-distribution-dependent compromise.

Abstract

from arXiv · show

We present a simple decomposition of a neural-network-based normalizing flow that naturally uncovers a pivotal statistic (or something close) in the presence of nuisance parameters, based only on a sample generator from the distribution of interest. We show that the statistic is near-pivotal in the sense of minimum average KL-divergence of its $p$-values versus uniform and we argue that it can be expected to have good power when the dimension of the statistic equals the dimension of the parameter. It is able to incorporate prior knowledge about group invariances such as translation and scale. It can discover the one-sample $t$-test almost exactly, outperforms the Welch test in terms of worst-case size over a constrained variance-ratio range and achieves good calibration on partial biserial correlations, while showing higher power (and being much faster) on small-to-moderate samples than profile likelihood-ratio techniques.

1 Introduction

Frequentist inference with nuisance parameters seeks tests whose size is invariant across nuisance configurations, but general approaches rely on approximations or constraints. NeuralCIs addresses this with amortized neural inference trained from simulations, targeting pivotal statistics and p-values without requiring an explicit likelihood.

Similar tests aim to keep test size invariant across all nuisance-parameter configurations.
Profile likelihood ratios and parametric bootstraps provide general solutions, but rely respectively on asymptotic approximations or fitted-model sampling.
Neural networks can learn likelihood-ratio information or approximate likelihoods directly from simulated samples when explicit likelihoods are unavailable.
NeuralCIs trains one network for an entire problem class, amortizing an initial two-to-four-hour training cost across rapid p-value evaluations.
Its central contribution is a normalizing-flow decomposition that yields an on-average near-pivotal statistic with p-values minimizing average KL-divergence to uniform.

2 Normalizing flows

Normalizing flows learn invertible mappings between data and simple target distributions, while this implementation uses back-propagated Jacobians and conditional inputs for flexible simulator-based modeling. Reversing the mapping yields a denormalizing flow used to sample training parameters from unnormalized targets.

A normalizing flow uses an invertible neural change of variables to map complex data distributions toward a standard normal target.
Training minimizes KL divergence by maximizing target log likelihood across samples, with the change-of-variables Jacobian contributing to the objective.
Unlike typical flows, this implementation computes Jacobian determinants by back-propagation and penalizes zero or negative determinants to encourage approximate invertibility.
The Jacobian calculation costs O(n_x^3), but the approach supports conditioning variables and flow decomposition without exotic constrained architectures.
The present approach is computationally expensive and can target only a local diffeomorphism, motivating more scalable methods for higher-dimensional problems.
Adding θ as network inputs produces conditional densities p(x|θ), allowing a different mapping for each parameter value.
Reversing the flow maps standard-normal samples to an unnormalized target and is called a denormalizing flow for parameter sampling.

3 Decomposing the normalizing flow

The proposed decomposition separates interest-parameter information from nuisance information by restricting the pivoting network to the scalar interest function. Training then makes the pivot approximately standard normal, from which calibrated one- and two-tailed p-values are obtained.

The paper’s key contribution is a normalizing-flow decomposition designed to uncover pivotal quantities and their distributions.
The construction assumes the summary dimension equals the parameter dimension, n_x = n_θ, while extensions to n_x ≥ n_θ remain future work.
The flow is split into a pivoting net z_p and a nuisance net z_n, whose outputs are concatenated into z ∈ R^n_x.
Training encourages Z_p|θ to be as close as possible to standard normal on average over the training distribution, measured by KL divergence.
An added constraint discourages orientation flips or folding in z_p and encourages alignment with the nuisance-free component of the ψ-action.
The learned z_p is transformed through the standard normal CDF into one-tailed p-values and then into equal-tailed two-tailed p-values.
For the one-sample t-test example, v = Φ(z_p(m,s,μ_0)) gives a p-value for H_0: μ = μ_0.

4 Formal properties

The formal analysis characterizes the learned statistic as near-pivotal under regularity and dimensionality assumptions, with p-values optimized for average KL proximity to uniform. Power arguments connect the statistic to nuisance removal, local quotient directions, efficient scores, and profile likelihood ratios, while global non-pivotal trade-offs remain possible.

4.1.1 Suitable models: The formal setup assumes dimension matching, nx = nθ, and full-rank transport of probability mass under parameter perturbations.These models are called “full-rank models.”
4.1.2 Standard regularity assumptions: The analysis further assumes a smooth density, regular coordinate patches, finite-variance scores, and full-rank Fisher information.The learned transformation is also assumed one-to-one for the formal change-of-variables arguments, although the implementation does not guarantee global invertibility.
4.2 Near-uniformity (by KL) of p-values (type I error): In the infinite-sample limit, constrained-flow likelihood maximization minimizes average KL divergence between the learned normal variables and a standard normal.A sufficiently flexible nuisance component can Gaussianize the conditional distribution, leaving the pivotal component’s divergence as the optimized term.
4.2 Near-uniformity (by KL) of p-values (type I error): The resulting p-values have minimum average KL divergence from a uniform distribution because KL is invariant under bijective conversion from the standard-normal statistic to p-values.The construction converts the statistic into one- or two-tailed p-values after normalizing transformation.
4.2.1 Suitability of KL optimization: KL optimization provides average near-pivotality and can improve tail calibration, but it may trade severe non-pivotal behavior at some parameter values against better behavior elsewhere.The loss is not tailored to tail-specific calibration or worst-case guarantees, especially for low-probability parameter regions.
4.3.2 Near-transformation models: For near-transformation models, nuisance perturbations leave a local quotient direction containing pure interest-parameter information, making the learned statistic locally pivotal.Globally, the model must stitch together potentially inconsistent local quotient directions, producing a smooth compromise shaped by training emphasis.
4.3.3 The local and asymptotic case: The local efficient-score construction is asymptotically pivotal, making it an ideal target for the learned statistic and yielding the χ2 limit of −2 log Λ.This connects the learned target to classical profile-likelihood behavior in the local asymptotic regime.
4.3.3 The local and asymptotic case: Under local asymptotic normality, efficient scores, profile likelihood ratios, and local maximal invariants share the same local information about the interest parameter.The maximal invariant depends on the data through the efficient score, while the profile likelihood ratio is a function of that invariant.

5 Further considerations

The method incorporates known group invariances and conditioning variables into the flow while carefully selecting training parameters. Architectural changes reduce input complexity and support interactions, with empirical evidence of slightly improved test-size accuracy for the proposed layer design.

Incorporating invariances: Known translation and scale invariances can reduce network complexity while expanding the domain over which valid p-values are produced.The one-sample t-test motivates invariance to location and scale.
Incorporating invariances: A canonicalization rule maps each statistic to a representative orbit point, reducing the transformed input dimension from nx to nx − nG while retaining nθ parameter dimensions.The parameter transformation uses the same data-dependent group element.
Incorporating invariances: The flow computes the Jacobian with respect to the original x through the canonicalization map, allowing a full-rank Jacobian despite dimensionality reduction.This requires nG ≤ nx and relies on gradient transfer through the parameter branch.
Incorporating invariances: The implementation reproduces the t-test example exactly and reports promising numerical results for a case with nx = nθ = 3 and nG = 2.The example demonstrates that Jacobian information can flow through the parameters after canonicalization.
Known values: Additional known values, such as sample size n, can condition the flow so the t-test problem is modeled across multiple sample sizes.These conditioning variables are added as further network inputs.
Parameter sampling: Training uses an inner parameter set covering target data and an outer set covering parameters that can generate overlapping data, making parameter-sample selection central to valid p-values.The procedure approximates Jeffreys-prior sampling using Fisher information estimated from repeated parameter estimates.
Parameter sampling: The estimates-box procedure defines the target region around parameter estimates, but its current construction is described as somewhat ad hoc and slated for refinement.The same data-to-estimate mapping is used to define the bounding box.
Network design: The proposed hidden-layer design combines a Hadamard interaction, layer normalization, and addition to support direct input interactions, stable propagation, and skip-like gradients.In the Behrens–Fisher simulation, standard deviation of test sizes at α = 0.05 was 0.0032 with multiplication versus 0.0035 with ELU.

6 Experimental Methodology

The experiments evaluate NeuralCIs across classic nuisance-parameter problems, testing calibration, power, invariance, and robustness under constrained sampling designs. They compare NeuralCIs with established tests and stress-test canonicalization and non-sufficient or discrete statistics.

Invariance and maximal invariant: The study also examines whether the learned statistic aligns with the maximal invariant using canonicalized and non-canonicalized one-sample t-test models.Canonicalization incorporates location and scale invariance, while the toy experiment assesses convergence toward the maximal invariant.
Evaluation protocol: Power comparisons use paired random samples and size-adjusted thresholds so liberal or conservative calibration does not directly determine the apparent power ranking.The adjustment estimates an alternative-specific p-value threshold rather than using α mechanically.
One-sample t-test: 50000 parameter samples are used to evaluate t-test calibration across µ, σ, and sample-size ranges, with a narrower main-text domain for the non-canonicalized model.The reported non-canonicalized experiments use µ ∈ [−3, 3], σ ∈ [0.333, 3], and n ∈ [3, 100], while the canonicalized model is tested over much broader parameter ranges.
Behrens-Fisher: NeuralCIs are compared with Welch’s test on the Behrens-Fisher problem using shared statistics and explicit size-adjusted power calculations.The model uses unbiased sufficient statistics, known group sample sizes, and location-scale normalization, reducing the statistic input to the estimated variance ratio.
Stress testing: The Behrens-Fisher stress tests sample fresh parameter combinations outside the training range, while retaining a constrained variance-ratio range because of edge weaknesses in parameter sampling.The stated constrained range is σ2/σ1 ∈ [0.333, 3], with expansion to the full inner parameter space left for later work.
Partial biserial correlation: The partial biserial experiment tests correlations and a binarization threshold in a model whose statistics are neither sufficient for the parameter vector nor entirely continuous.This deliberately probes whether the flow-based approach can approximate the discrete case despite its continuous-diffeomorphism foundation.

7 Results

NeuralCIs closely reproduce the one-sample t-test, improve worst-case size relative to Welch, and achieve calibrated error rates with higher power than profile likelihood-ratio methods in small-to-moderate samples. The method also remains well behaved near parameter boundaries, while sample-size comparisons require caution because parameter sampling regions vary.

7.1 t-Test: NeuralCIs reproduced one-tailed t-test p-values extremely closely, with canonicalized differences below 0.0013 in all samples.The canonicalized model stayed within 0.0012 for 99.5% of samples.
7.2 Behrens-Fisher: NeuralCIs had considerably better worst-case size error than Welch, with negligible power difference below 0.1 percentage points.Welch showed a slightly narrower peak, but NeuralCIs had the more favorable worst-case size behavior.
7.3 Partial biserial correlation: At n ∈[20, 100], the χ2 likelihood-ratio test was consistently liberal, while NeuralCIs clustered near the target rate and averaged 1.2 percentage points greater power.This comparison concerns the partial biserial correlation problem.
7.3 Partial biserial correlation: Against bootstrap likelihood ratio, NeuralCIs and the comparator were centered near nominal error rates, while NeuralCIs retained a 1.2-percentage-point size-adjusted power advantage.The comparison used α = 0.01 and α = 0.05; bootstrap estimates were based on only 500 p-values per parameter.
7.3 Partial biserial correlation: NeuralCIs’ advantages over profile likelihood-ratio methods were strongest at sample sizes 20–70, with about a 3-percentage-point average power advantage at sizes 30–50.The advantage narrowed or reversed near n = 100, and interpretation is limited because the parameter sampling region grows with sample size.
7.3 Partial biserial correlation: NeuralCIs maintained error rates close to desired α across broad null-parameter ranges without boundary degeneracy, generating ten billion p-values in fifteen minutes.Profile-likelihood approaches also exhibited non-convergence and boundary issues, especially for bootstrap LR.
7.4 Deviation from the maximal invariant: In the maximal-invariant toy problem, the unexplained-variance proportion became negligible during training and outputs converged to variance 1 with low KS distance.Seven of 250 initializations converged less cleanly, but all improved after 500 additional epochs.

8 Discussion

NeuralCIs performs competitively across classical nuisance-parameter problems, with advantages in worst-case error, power, computation, and boundary behavior. The discussion also identifies current limitations and outlines a route from p-values to confidence intervals.

The Behrens-Fisher method has similar power to Welch while achieving considerably better worst-case false positive rates within pre-specified variance-ratio bounds.
NeuralCIs matches bootstrap profile likelihood-ratio performance for partial biserial correlation at a tiny fraction of the computational cost, with a 1.2-percentage-point average power advantage.The advantage is around 3 percentage points for sample sizes from 30–50.
NeuralCIs generates and performs well on p-value distributions near boundaries, where the profile likelihood-ratio framework required abandoning attempts to do so.
The t-test example nearly recovers the classical pivot, while the maximal-invariant toy problem shows this behavior consistently even when other pivots exist.
Current limitations include edge-of-parameter-space sampling, ad hoc network dimensions, unassessed training variability, and the need for further confidence-interval and statistical-problem evaluations.
A further network could invert the p-value statistic into confidence intervals by mapping statistics and target normal values to compatible interest-parameter boundaries.

9 Conclusion

The paper concludes that decomposed normalizing flows can extract near-pivotal statistics and support competitive frequentist inference on classical nuisance-parameter problems. The fitted network also offers invariance handling and inexpensive subsequent use, but its reasonable-power regime remains dimensionally constrained.

Decomposed normalizing flows extract near-pivotal statistics and achieve competitive accuracy across three classical problems, sometimes with advantages over classical approaches.
The network can incorporate invariance information, potentially improving robustness and extending its reach.
After fitting, the network is much cheaper to use than bootstraps or likelihood-surface optimizers.A major limitation is that reasonable power requires n_x = n_θ.

Declarations

The paper reports independent authorship, open code availability, no external funding, and no competing interests.

The core contribution and codebase were developed independently by the author, with language models used during subsequent theoretical development.
The project code is available at github.com/philassheton/neuralcis.
The work received no external funding.
The author declares no competing interests.

A Parameter sampling details

The parameter-sampling procedure is presented as an engineering section that supports NeuralCIs sampling, while documenting sampling targets, boundary restrictions, and estimator assumptions. The authors also identify implementation choices and known sampling issues for future improvement.

The engineering section can be skipped for a rough understanding of the model and is expected to improve in future work.
The sampling scheme targets future confidence-interval generation rather than directly guarding against inappropriate null values for p-values.
X_target is defined using a parameter-estimate bounding box, retaining statistics whose estimates fall within that box.
The denormalizing-flow sampler is built through MCMC-style walks, neural networks, and successive flow fits targeting inner and outer parameter regions.
The current sampling procedure has an issue for particularly wide-sampled estimates because its epsilon value must depend on the modeled square-root Fisher information.
Alternative estimators may be used instead of maximum-likelihood estimates, but their distributions must transform under group actions when invariances are incorporated.

B Experimental sampling for Partial biserial correlation

The evaluation samples parameters within bounded regions designed to limit boundary hits, while separately sampling power targets for partial biserial correlations. The section also connects the scalar flow coordinate to maximal invariants and profile likelihood ratios under transformation-model assumptions.

Parameter sampling: The evaluation uses a narrow parameter-sampling scheme to avoid regions where baseline likelihood-ratio optimization frequently reaches boundary values.The sampled region differs from the training distribution and is deliberately constrained for smaller sample sizes.
Parameter sampling: Parameter bounds for prevalences and correlations are constructed using a target boundary rate, Fisher’s r-to-z transform, and point-biserial conversion factors.The partial biserial correlation bounds additionally adjust for the sample-size degrees of freedom.
Parameter sampling: Parameters are sampled uniformly between computed minima and maxima, including sample size, while separate partial biserial correlations target power values from 0.05 to 0.90.The power-target construction samples β uniformly from 0.10 to 0.95 and converts it into corresponding alternative correlations.
Transformation-model connection: In a full-rank, dimension-matched transformation model, the scalar flow coordinate zp is linked to the one-dimensional maximal invariant associated with the interest parameter.The remaining coordinates absorb nuisance directions, and the profile likelihood ratio is invariant under the joint group action.
Transformation-model connection: Because the profile likelihood ratio is invariant under the group action, it can be expressed as a scalar function of the maximal invariant.The derivation relies on cancellation of the common Jacobian factor in the likelihood ratio.
Transformation-model connection: For the Gaussian location model, the profile log likelihood ratio is obtained by optimizing the quadratic likelihood over nuisance directions and is a function of the maximal invariant.The translation example illustrates this relationship explicitly.

C.2 Why we think zp will tend towards mψ

The paper argues that zp tends toward the maximal invariant because alternative pivotal solutions require globally coordinated nuisance-dependent cancellations. Such solutions are possible, but random initialization and stochastic-gradient training are unlikely to produce substantial versions of them.

Targeting the maximal invariant: A Gaussianized maximal invariant is a natural target for zp because nuisance transformations move probability within invariant leaves while leaving the distribution across leaves pivotal.The nuisance coordinates zn model within-leaf structure, while p(mψ) remains independent of λ.
Alternative solutions: The loss can admit zp solutions that depend on nuisance coordinates u while remaining standard normal across λ, but these require considerable coordinated cancellations.The paper presents this as a possible loss-optimal alternative rather than the typical learned solution.
Alternative solutions: For monotone zp, nuisance dependence must be offset across regions as λ shifts probability mass through u-space, imposing global structure across all z and λ.The periodic uniform-shift example demonstrates how such cancellation can occur.
Alternative solutions: Non-monotone zp can create additional cancellation through multiple branches, but the diffeomorphism constraint requires coordinated smooth, non-intersecting contour families.A single fold cannot stand alone because the contours must remain continuous and non-crossing.
Conclusion: The paper therefore expects only a very small nuisance-dependent component in zp, because random initialization and training are unlikely to generate the required global coordination.The loss may be indifferent to some structured dependencies, but most random nuisance dependence is expected to be removed during training.
Alternative solutions: The architecture’s interest-direction constraint further requires calibration to balance probability flux as ψ changes, adding another coordination requirement for nuisance-dependent solutions.The flux into level sets must be counterbalanced by shifts in zp values as the data distribution moves.

C.2.3 Consequence of contamination for power

The paper models residual nuisance dependence in zp as independent standard-normal contamination and finds that its effect on noncentrality and power is only second-order in the contamination amplitude.

Contamination model: The likely nuisance-dependent component of zp may be uninformative about ψ or pure noise, motivating an independent standard-normal contamination benchmark.The benchmark treats the contaminating component as carrying no ψ-signal.
Power consequence: The contamination divides the alternative noncentrality by a factor involving 1 + ϵ^2, so the resulting power loss is only second-order in ϵ.Power is described as a smooth function of noncentrality in this benchmark.

D.1 t-Test: the edges of Θinner

The t-test evaluation examines performance beyond the estimates bounding box and diagnoses edge degradation, sampling-scheme weaknesses, likelihood-fit failures, boundary cases, and slower-converging training runs.

Parameter-region evaluation: The extended evaluation tests whether valid p-values persist throughout Θinner beyond the estimates bounding box by sampling Bonferroni-adjusted regions around its boundary.The construction uses separate parameter-wise p-values and retains points allowed by the combined confidence criterion.
Parameter-region evaluation: Larger sample sizes perform relatively well across confidence levels, whereas small-n cases worsen progressively toward the box edge.The paper attributes this pattern to weaknesses in the current parameter-sampling scheme for wider-spread estimates.
Sampling limitations: The current sampling scheme can overweight tightly spread parameters when the bounding or target region is large relative to the estimate distribution.The authors identify this as a potential catastrophic distortion and propose weighting by collision probabilities over regions rather than single points.
Sampling limitations: For small σ values, equal sampling probability underweights parameters whose sampling distributions overlap less of the bounding box, while extending the box can make the correction dominate.This disproportionate effect is strongest at small sample sizes.
Numerical diagnostics: The diagnostics count likelihood non-convergence, near-boundary estimates, and degenerate binary samples across main and size-adjust simulations.A degenerate sample has all binary A* values equal, making the biserial correlation undefined.
Training stability: Further training brought the seven slower-converging runs into the same neat convergence pattern as the other 243 runs.The additional training used a fresh set of 500 epochs of 100 steps and restarted from the initial higher learning rate.