Model Fitting

Top-level compiled model.

dict_to_model(model_dict, list_free_params; nB_workspace=100) -> FlatModel

Compile a flat parameter dict into a FlatModel ready for evaluation.

Arguments

modeldict : Dict{String,Any} with "component,param" keys. Values may be Float64 (fixed or free) or String (expression). Special string key "component,profile" triggers the Hankel pathway. listfreeparams : Vector{String} of keys to optimize (must be numeric in modeldict). nB_workspace : Pre-allocated baseline grid size for HankelWorkspace buffers. Should be ≥ the number of baselines you will evaluate at.

Returns

A FlatModel whose evalmodel / evalmodel_grad methods can be called in a loop with zero Dict allocation.

Example

model_dict = Dict(
    "star,ud"      => 0.8,
    "star,f"       => 0.6,
    "ring,profile" => "exp(-(\$R / \$scale)^2 / 2)",
    "ring,scale"   => 1.5,
    "ring,udout"   => 12.0,
    "ring,f"       => "1 - \$star,f",
)
list_free_params = ["star,ud", "star,f", "ring,scale"]
model = dict_to_model(model_dict, list_free_params)

Backward-compatible alias for dict_to_model.

model_to_vis(model, x, uv; wl=nothing, mjd=nothing, n=0)

Alias for eval_model. Compute complex visibilities from a parametric model, for consistency with the model_to_obs / model_to_chi2 family.

eval_model(model, x, uv; wl=nothing, mjd=nothing, n=0) -> Vector{Complex{T}}

Evaluate the model visibility at the given parameter vector x and baselines uv.

Arguments

model : FlatModel from dicttomodel x : current free-parameter values, length = length(model.listfreeparams) uv : 2×N matrix of (u,v) spatial frequencies in cycles/rad wl : wavelength per UV point (metres), length N — enables $WL in expressions mjd : MJD per UV point, length N — enables $MJD in expressions n : Bessel order for Hankel components (default 0)

Returns

V : complex visibility vector, length N

eval_model_grad(model, x, uv; wl=nothing, mjd=nothing, n=0) -> (V, J)

Evaluate model visibility and its Jacobian w.r.t. x.

Returns

V : complex visibility vector (length nB) J : Jacobian, Matrix{Complex{T}} of shape (nB, length(x)); T follows x and uv

Uses ForwardDiff for the resolver and analytic components; for Hankel components uses the cached-K chain rule via hankelvisfwd! + pullback. This is more efficient than naive ForwardDiff through the full Hankel loop when nB is large, since Bessel evaluations are shared across parameters.

display_model(model_dict, list_free_params; lb, ub)

Print a human-readable summary of the model setup, grouped by component. Shows each parameter's value (or expression), whether it is free or fixed, and its bounds. Warns about common mistakes (value out of bounds, missing bounds, flux fractions that don't sum to 1).

Call this before fit_model / fit_model_ultranest to verify your setup.

default_bounds(model_dict, list_free_params; data = nothing, max_size = nothing) -> (lb, ub)

Suggest lower and upper bounds for each free parameter, as two Dict{String,Float64} ready to pass to any fitter's lb/ub keyword.

Bounds are chosen from the parameter's suffix — the part after the first comma, so "star,ud" is looked up as "ud". Azimuthal-mode keys ("c,az amp1", "c,az projang1") are recognised by prefix. Anything unrecognised, including bare global names, gets (-Inf, Inf); callers that need finite bounds (nested sampling) should check for that.

Angular sizes and offsets are the only bounds worth tuning, and they are the ones the data can tune. Pass data and the ceiling becomes 2 λ/B_min from the actual uv coverage — the array cannot constrain a source larger than the scale its shortest baseline senses. Without data the ceiling is DEFAULT_MAX_SIZE_MAS; max_size overrides both.

lb, ub = default_bounds(model_dict, free)              # 100 mas ceiling
lb, ub = default_bounds(model_dict, free; data = data) # sized to this uv coverage

Everything else is physical rather than instrumental — flux fractions in [0,1], inclination in [0,90], position angles in [-180,180] — and does not depend on the data.

See also display_model, which validates values against whatever bounds you end up with, and max_angular_scale.

Default upper bound on an angular size, in mas, when no data is supplied.

Long-baseline interferometry senses scales between roughly λ/2Bmax and λ/Bmin. Measured on three real datasets, that is 0.5–1.2 mas at the fine end and 18–41 mas at the coarse end: a source much larger than λ/Bmin is fully resolved out and its size is not constrained at all. 100 mas therefore leaves generous headroom over anything an array can measure, while staying a real bound. Pass data to `defaultbounds` to size this from the actual uv coverage instead.

max_angular_scale(data) -> Float64

Largest angular scale the uv coverage is sensitive to, in mas: λ/B_min, from the shortest baseline present. A source larger than this is fully resolved and its size is unconstrained, which makes it the natural ceiling for a size bound.

ModelConstraint(lhs, op, rhs; tol = DEFAULT_CONSTRAINT_TOL)

A relation between model parameters, enforced during fitting.

lhs and rhs are parameter names ("disk,diamout"), expressions over parameter names ("2 * star,ud"), or — for rhs — a plain number. op is one of "<", "<=", ">", ">=", "=", given as a string or a symbol. tol is how much violation counts as none; see DEFAULT_CONSTRAINT_TOL.

fit_model hands these to NLopt as real nonlinear constraints, so they hold at the optimum. fit_model_lsqfit and fit_model_ultranest have no such machinery and apply a one-sided quadratic penalty instead, which is soft: a steep enough χ² can overrule it, and tol then sets the stiffness rather than a feasibility tolerance.

ModelConstraint("disk,diamout", ">", "disk,diamin")   # a ring must have an outside
ModelConstraint("star,f", "+", 0.0)                   # ── invalid: `+` is not a relation
ModelConstraint("star,f + disk,f", "=", 1.0)          # fluxes sum to one
ModelConstraint("star,ud", "<", 4.0; tol = 0.01)      # deliberately loose

Globals — model keys with no comma — must be written with a \$, as "\$PA", exactly as in a model expression. A bare name with no comma and no arithmetic is taken to be one and gets its \$ added.

See also fit_model and read_model_file.

parse_constraints(specs) -> Vector{ModelConstraint}

Accept constraints in any of the forms a caller is likely to have them in.

A ModelConstraint passes through. A tuple is read positionally as (lhs, op, rhs) or (lhs, op, rhs, tol) — the PMOIRED prior layout, so a converted PMOIRED model needs no reshaping. A Dict or NamedTuple is read by the keys param/lhs, op, value/rhs and tol, which is the shape a TOML file parses to.

check_constraints(constraints, model_dict; verb = true) -> Vector{Bool}

Report which constraints the model satisfies as it stands, without fitting anything.

Useful before a long run: a starting model that violates a constraint is not an error — the penalty will push it back — but it is worth knowing about, because it may equally mean the constraint was written backwards.

Default tolerance for a constraint, in the units of the quantity being constrained.

It means two related things, one per enforcement mechanism. To NLopt it is the feasibility tolerance: how far outside the constraint a point may sit and still count as satisfied. To the penalty used by the least-squares and nested-sampling fitters it is the stiffness: a violation of tol costs one unit of penalty.

1e-3 is a thousandth of a milliarcsecond on an angular size, and a tenth of a percent on a flux fraction — below anything interferometry resolves, in both readings.

read_model_file(path) -> (; model, free, lb, ub, constraints, priors, name)

Read a TOML model file into the arguments a fitter takes.

The result destructures straight into any of the fitters:

m = read_model_file("binary.toml")
res = fit_model(m.model, m.free, data; m.lb, m.ub, m.constraints, m.priors)

Every section is optional. A file with only [model] reads back as a model with no free parameters, which is what model_to_image and model_to_sed want. Bounds absent from the file are not filled in from default_bounds — an absent bound means ±Inf, and silently substituting a guess would make two fits from the same file differ by which version of the defaults was in the package.

See also write_model_file, ModelConstraint and default_bounds.

write_model_file(path, model_dict; free, lb, ub, constraints, priors, name)

Write a model and its fit settings to a TOML file readable by read_model_file.

write_model_file("binary.toml", model; free, lb, ub,
                 constraints = [ModelConstraint("disk,diamout", ">", "disk,diamin")])

Keys are sorted, so the same model always produces the same bytes and two files can be diffed. Bounds that are infinite on both sides are omitted rather than written as [-inf, +inf]: they constrain nothing, and every fitter already treats a missing bound that way.

A FitResult is not what gets written — the fitted values are, if the caller passes them in model_dict. Writing the result of a fit back out is therefore write_model_file(path, merge(model, Dict(zip(res.list_free_params, res.x_opt))); free = res.list_free_params, ...).

fit_model(model_dict, list_free_params, data; kwargs...) -> FitResult

Fit a parametric model to interferometric data using NLopt.

Arguments

• model_dict::Dict{String} — flat parameter dictionary (values are numbers or expression strings)
• list_free_params::Vector{String} — names of the free parameters to optimize
• data::OIdata — interferometric data

Keywords

• lb, ub — Dict{String,Float64} of lower/upper bounds per parameter (default: ±Inf)
• weights — observable weights [V2, T3amp, T3phi, visamp, visphi, flux, diffphase]
• priors — Vector of (expr_str, target, sigma) tuples for Gaussian penalties
• constraints — relations between parameters that bounds cannot express, as ModelConstraints or (lhs, op, rhs[, tol]) tuples. Handed to NLopt as real nonlinear constraints, so they hold at the optimum; a method that cannot take them is wrapped in :AUGLAG rather than replaced. The reported chi2 is the data χ², unaffected
• method — NLopt algorithm (default :LD_LBFGS; use :LN_NELDERMEAD for gradient-free)
• maxeval — maximum function evaluations (default 2000)
• ftol_rel, xtol_rel — convergence tolerances
• vonmises — use von Mises statistic for T3phi
• nB_workspace — Hankel workspace size (default: nuv from data)
• verb — print per-evaluation chi2 breakdown

fit_model(model::FlatModel, x0, data; kwargs...) -> FitResult

Fit a pre-compiled FlatModel starting from parameter vector x0.

This method skips the dict_to_model step, which is useful when fitting the same model multiple times (e.g. bootstrap, grid search).

Bounds lb and ub can be Dict{String,Float64} (keyed by parameter name) or Vector{Float64} (ordered to match model.list_free_params).

fit_model_lsqfit(model_dict, list_free_params, data; kwargs...) -> LsqFitResult

Fit a parametric model to interferometric data using Levenberg-Marquardt (LsqFit.jl). Returns parameter covariance and 1σ uncertainties.

Uses the analytic Jacobian from residuals_flat_jac (Wirtinger chain rule through the complex visibility Jacobian), so no finite-difference overhead.

Arguments

• model_dict::Dict{String} — flat parameter dictionary
• list_free_params::Vector{String} — names of the free parameters to optimize
• data::OIdata — interferometric data

Keywords

• lb, ub — Dict{String,Float64} of lower/upper bounds per parameter
• weights — observable weights [V2, T3amp, T3phi, visamp, visphi, flux, diffphase]
• vonmises — use von Mises statistic for T3phi
• nB_workspace — Hankel workspace size
• maxIter — maximum iterations (default 200)
• constraints — relations between parameters that bounds cannot express, as ModelConstraints or (lhs, op, rhs[, tol]) tuples. They enter as extra rows of the residual vector, which is exactly where Levenberg-Marquardt can act on them — and so also enter the Jacobian, meaning covar and stderror account for them. The reported chi2 excludes them, so it stays comparable with what every other fitter reports
• verb — print per-evaluation chi2 breakdown

fit_model_lsqfit(model::FlatModel, x0, data; kwargs...) -> LsqFitResult

Levenberg-Marquardt fit on a pre-compiled FlatModel starting from x0. Bounds lb and ub can be Dicts or Vectors.

fit_model_nested(model_dict, list_free_params, data; kwargs...) -> NestedResult
fit_model_nested(model::FlatModel, data; kwargs...)             -> NestedResult

Fit a parametric model by Bayesian nested sampling, returning the posterior and the evidence as well as a best-fit point.

Unlike the gradient fitters this needs no starting guess — it samples the whole prior volume — but it does need finite bounds on every free parameter, since those bounds are the prior.

using OITOOLS, Nautilus
r = fit_model_nested(model_dict, ["star,ud"], data;
                     lb = Dict("star,ud" => 0.5), ub = Dict("star,ud" => 10.0))
r.logz, r.logzerr        # the evidence, for comparing models
r.posterior              # equally weighted samples, (n_samples, n_params)

backend selects the sampler and defaults to nested_backend(); see set_nested_backend! to change it for a session.

Keywords, both backends

Keywords, :nautilus

Keywords, :ultranest

Evidence from two different samplers is worth more than evidence from one: logz should agree within logzerr, and NestedResult.backend records which produced any given number.

fit_model_ultranest(args...; kwargs...) -> NestedResult

fit_model_nested pinned to the UltraNest backend. Needs using PythonCall and the Python ultranest package (declared in CondaPkg.toml).

nested_backend() -> Symbol

Which nested sampler fit_model_nested will use: :nautilus, :ultranest, or :none when no extension is loaded.

set_nested_backend!(b::Symbol) -> Symbol

Choose the nested sampler: :nautilus or :ultranest.

The backend's package must already be loaded, and this says so immediately rather than letting a long run end in a MethodError.

chi2_map(model_dict, list_free_params, data, p1, p2; kwargs...) -> Chi2Map

Evaluate χ² on a grid over the two free parameters p1 and p2, holding every other free parameter at its starting value.

m = read_model_file("LD_power.toml")
map = chi2_map(m.model, m.free, data, "star,ldpow", "star,alpha"; m.lb, m.ub, weights = [1.0,0,0])
res = FitResult(map)          # the best grid point, as any other fitter would report it

An axis with no finite bound and no explicit range is an error rather than a guess: a grid needs an extent, and inventing one would quietly decide the answer.

p1 and p2 must both be free. Mapping a fixed parameter would be a perfectly reasonable thing to want, but x[i] ↔ list_free_params[i] is what makes the result a FitResult at all, so free it first.

delta_chi2_levels

Δχ² for the 68.3, 95.4 and 99.7% confidence regions of two jointly estimated parameters: 2.30, 6.17, 11.8.

These are the two-parameter values, not the one-parameter 1, 4, 9. Contouring a two-parameter map at 1σ = 1 draws a region about a third too small, which is the usual way a χ² map is misread.

model_warnings(model_dict, list_free_params; lb, ub) -> Vector{String}

Everything wrong with a model that can be seen without fitting it, one sentence each.

Four checks, and every one of them describes a setup that runs happily and answers the wrong question: a starting value outside its own bounds, bounds that exclude everything, and flux fractions that do not sum to one. The last is the quiet one – a second component added with f = 1 doubles the model's flux, and the only symptom is a χ² in the millions.

Shared with display_model, which prints these to a terminal, so a GUI and a script cannot disagree about what is wrong with the same model.

Flux fractions given as expressions are not checked: their values depend on the resolver, so summing the literals would be summing the wrong things.

model_to_obs(model, x, data) -> (v2, t3amp, t3phi, visamp, visphi)

Evaluate a FlatModel at parameter vector x and return the model observables matching the data structure.

Returns a NamedTuple with fields v2, t3amp, t3phi (degrees), visamp, visphi (degrees). Empty vectors are returned for observable types not present in data.

model_to_residuals(model, x, data)

Compute normalised residuals (model - data) / error for each observable type. Returns a NamedTuple with fields v2, t3amp, t3phi, visamp, visphi. Phase residuals (t3phi, visphi) are wrapped to [-180, 180] before dividing by error.

model_to_chi2(model, x, data; weights=[1,1,1,0,0,0,0], verb=false, vonmises=false)

Alias for chi2_flat. Compute the weighted chi-squared of a parametric model against data, for consistency with the model_to_obs / model_to_vis family.

model_to_chi2_fg(model, x, data; weights=[1,1,1,0,0,0,0], verb=false, vonmises=false)

Alias for chi2_flat_fg. Compute chi-squared and its gradient w.r.t. x, for consistency with the model_to_obs / model_to_vis family.

model_to_image(model, x; nx=256, pixsize=0.1, oversample=1, normalize=true, wl=nothing)
    -> Matrix{Float64}

Synthesise an image from a FlatModel at parameter vector x.

The image is computed by evaluating the model visibility on the 2-D FFT frequency grid and applying an inverse real FFT (irfft), which exploits the Hermitian symmetry V(-u,-v) = conj(V(u,v)) to halve the number of visibility evaluations.

Arguments

• model::FlatModel — compiled model from dict_to_model
• x::AbstractVector — current free-parameter values

Keywords

• nx — image size in pixels (default 256)
• pixsize — pixel size in mas (default 0.1)
• oversample — oversampling factor (default 1)
• normalize — normalize image to unit sum (default true)
• wl — wavelength in microns (scalar); enables $WL in spectrum expressions

Returns

img — nx × nx Matrix{Float64}

model_to_sed(model, x, wl_grid) -> (total, components)

Evaluate the spectral energy distribution of a FlatModel by computing V(u=0, v=0) at each wavelength — the zero-baseline visibility equals the total flux.

Per-component fluxes are obtained by resolving the f / spectrum parameters at each wavelength.

Arguments

• model::FlatModel — compiled model from dict_to_model
• x::AbstractVector — current free-parameter values
• wl_grid — wavelength array in METRES, the unit $WL carries in the resolver and uv_lam carries in the data. A grid in microns evaluates every $WL expression at a wavelength a million times too large and returns fluxes that are wrong without erroring.

Returns

• total::Vector{Float64} — total flux at each wavelength (= V(0,0))
• components::Dict{String,Vector{Float64}} — per-component flux

model_to_flux(model, x; wl=nothing)

Total flux of the model at zero baseline: real(V(0,0)).

For polychromatic models, pass wl as a scalar or 1-element vector. See also model_to_sed for the full spectral energy distribution.

bootstrap_fit(model_dict, list_free_params, data; kwargs...) -> BootstrapResult

Estimate parameter uncertainties by nonparametric block bootstrap: refit the model to nboot replicates in which the blocks of data (by default one per MJD and telescope configuration) are resampled, and summarise the scatter of the best-fit parameters.

Unlike the analytic covariance of fit_model_lsqfit, this does not assume the quoted error bars are correct or uncorrelated: correlated calibration errors and mis-stated error bars show up as extra scatter between replicates. It does assume the blocks are numerous and exchangeable — with a single snapshot (a handful of blocks) the estimate is unreliable, and perturb_data based Monte Carlo or the analytic errors are the better tool.

Arguments

• model_dict::Dict{String} — flat parameter dictionary
• list_free_params::Vector{String} — names of the free parameters
• data::OIdata — interferometric data

Keywords

• nboot — number of replicates (default 200; use ≥1000 for stable percentiles)
• mode — :replacement (default), :halfsample or :weights; see block_counts and block_weights. :pmoired reproduces PMOIRED's scheme and is biased low by about √2 — comparison only
• granularity — :config (default), :epoch or :point; see data_blocks
• fitter — :lsqfit (default, Levenberg-Marquardt with analytic Jacobian) or :nlopt
• lb, ub, weights, vonmises, nB_workspace — passed through to the fitter
• method, maxeval — NLopt settings (fitter=:nlopt)
• maxIter — LsqFit iterations (fitter=:lsqfit)
• start_scatter — perturb each replicate's starting point by this many σ of the full-data fit (default 0; PMOIRED uses the equivalent of 1). Non-zero values mix optimiser scatter into the uncertainty and are meant for comparison studies only.
• sigma_clipping — reject replicates beyond this many σ from the median (default nothing; PMOIRED's plots default to 4.5)
• chi2r_max — reject replicates whose reduced χ² exceeds this value
• seed — RNG seed; replicate i uses Xoshiro(seed + i) so results are reproducible regardless of thread scheduling
• threaded — run replicates on all available Julia threads (default true)
• verb — progress reporting (default true)

Example

data  = readoifits("data/AlphaCenA.oifits")[1,1]
model = Dict{String,Any}("star,ldlin" => 8.0, "star,u" => 0.3, "star,f" => 1.0)
boot  = bootstrap_fit(model, ["star,ldlin", "star,u"], data;
                      lb=Dict("star,ldlin"=>5.0, "star,u"=>0.0),
                      ub=Dict("star,ldlin"=>12.0, "star,u"=>1.0),
                      weights=[1.0, 0.0, 0.0], nboot=500)

bootstrap_driver(fitfun, x_opt, list_free_params; kwargs...) -> BootstrapResult

Run nboot bootstrap replicates of an arbitrary fit and summarise the scatter of the resulting parameters.

This is the model-agnostic half of bootstrap_fit: it owns the replicate loop, the seeding, the threading, the failure handling, the clipping and the statistics, and knows nothing about what is being fitted or how the data are resampled. Packages with their own estimators (or their own notion of a resampling unit) can reuse the statistics by calling this directly.

Arguments

• fitfun(state, rng) -> (x, chi2r) — performs one replicate: it resamples the data and fits them, returning the best-fit parameters and the reduced χ². The driver never resamples; that is deliberate, so that callers are free to resample structures OITOOLS knows nothing about. A third element may be returned, (x, chi2r, extra), which is stored in info (see below).
• x_opt — best fit to the unresampled data, recorded in the result.
• list_free_params — parameter names; length sets the expected width of x.

Keywords

• nboot — number of replicates (default 200)
• worker_init(w::Int) -> state — build per-worker state, called once per worker and serially, before any task is spawned (compiling or allocating concurrently is exactly what this exists to avoid). state is nothing when worker_init is not given.
• nworkers — number of replicate-level workers; 0 (default) means min(nthreads(), nboot). Set it to 1 when fitfun threads internally, so the two levels do not oversubscribe.
• info — optional AbstractVector of length nboot; when given and fitfun returns a third element, info[i] receives it. Written at the replicate's own index, never in completion order.
• sigma_clipping, chi2r_max — reject replicates beyond this many σ from the median, or above this reduced χ²
• seed — replicate i uses Xoshiro(seed + i), so results do not depend on thread scheduling
• threaded — use several workers (default true; ignored if nworkers is set)
• verb — progress reporting
• mode, granularity, nblocks — metadata only, copied into the result

Example

# a fit with no data at all: recover the scatter of a known Gaussian
x_true = [1.0, 2.0]
r = bootstrap_driver((state, rng) -> (x_true .+ 0.1 .* randn(rng, 2), 1.0),
                     x_true, ["a", "b"]; nboot = 500, seed = 1, verb = false)
r.sigma   # ≈ [0.1, 0.1]

data_blocks(data::OIdata; granularity=:config, mjd_digits=5) -> DataBlocks

Partition data into resampling blocks.

Granularity

• :config (default) — one block per (MJD, telescope configuration), where the configuration is the baseline for V²/VIS, the triangle for T3 and the telescope for FLUX. All wavelength channels of a block stay together, and observables sharing a configuration at the same MJD (e.g. |V| and V² on the same baseline) stay together too. This is PMOIRED's "spectral vector".
• :epoch — one block per MJD: an entire observation, all baselines and all observables, is kept or dropped as a unit. Coarser, more conservative; needs many epochs to be usable.
• :point — one block per data point. The classical i.i.d. bootstrap; it destroys the spectral and per-configuration correlation structure and will therefore underestimate uncertainties on real data. Provided for comparison studies.

mjd_digits sets the rounding used to group MJDs into a common epoch (default 5 decimals ≈ 0.86 s, matching PMOIRED).

resample_blocks(data, blocks; mode=:replacement, rng=Random.default_rng()) -> OIdata
resample_blocks(data; granularity=:config, mode=:replacement, ...) -> OIdata

Bootstrap replicate of data: the observable values and their error bars are left untouched, and the blocks of blocks are drawn according to mode (see block_counts). A block drawn twice contributes its data twice; a block not drawn is absent.

The uv table is carried over unchanged — duplicated data points simply share the same uv row — so the model evaluation cost per replicate is unchanged.

block_counts(nblocks, mode; rng) -> Vector{Int}

Multiplicity drawn for each block under the given resampling mode. Exposed mainly for testing and for the calibration study of the different schemes.

Modes

• :replacement (default) — draw nblocks blocks with replacement. This is the textbook nonparametric bootstrap; multiplicities are ≈ Poisson(1), i.e. mean 1 and variance 1.
• :halfsample — keep a random half of the blocks, each once. Balanced repeated replication: correctly calibrated (the finite-population correction cancels the factor 2 from halving the data), and cheaper since each replicate fit uses half the data. In the validation study under demos/bootstrap_validation this was the best-calibrated scheme whenever blocks were numerous (ratio 0.94-1.06); with very few blocks it turns conservative (1.24-1.56), since each replicate then fits very little data.
• :pmoired — PMOIRED's scheme: build two independent copies of the dataset, drop a random half of the blocks in each, and fit the union. Multiplicities are 0/1/2 with probabilities 1/4, 1/2, 1/4, so their variance is 1/2 rather than 1.

resample_blocks additionally accepts mode=:weights, the multiplier (Bayesian) bootstrap, which draws continuous weights rather than integer multiplicities — see block_weights.

apply_block_counts(data, blocks, counts) -> OIdata

Build the OIdata in which each block of blocks appears counts[i] times.

block_weights(nblocks; rng) -> Vector{Float64}

Random weights for the multiplier ("Bayesian") bootstrap: nblocks draws from Dirichlet(1, …, 1) rescaled to mean 1, i.e. mean 1 and variance 1 − 1/nblocks, the same first two moments as the multiplicity of a draw with replacement.

Unlike block_counts, these are continuous: no block is ever dropped entirely and none is duplicated, so a replicate costs exactly one dataset instead of one-and-a-bit. See apply_block_weights.

apply_block_weights(data, blocks, w) -> OIdata

Multiplier bootstrap replicate: instead of duplicating and dropping blocks, give block i the weight w[i] in the χ² by dividing its error bars by sqrt(w[i]).

For a weighted least-squares fit this is the standard multiplier (or Bayesian) bootstrap: with mean(w) = 1 and var(w) = 1 it has the same asymptotic behaviour as resampling blocks with replacement, at the cost of one dataset per replicate and with no discrete jumps as blocks enter and leave.

The reduced χ² of such a replicate is not meaningful; the parameter scatter is.

perturb_data(data::OIdata; rng=Random.default_rng()) -> OIdata

Return a copy of data with independent Gaussian noise added to every observable, each draw scaled by that point's error bar. uv coordinates, flags and the number of data points are unchanged.

Useful for building noisy realisations of a dataset — simulation, end-to-end tests, sensitivity checks.

resample_data(data::OIdata)

Deprecated alias for perturb_data.

Despite the name, this function never resampled anything: it adds Gaussian noise drawn from the error bars. For a genuine bootstrap (resampling which observations are used) see bootstrap_fit and resample_blocks.

BootstrapResult

Outcome of bootstrap_fit.

Fields

• samples — nboot × npar matrix of per-replicate best-fit parameters
• list_free_params — parameter names, matching the columns of samples
• x_opt — best fit to the full dataset
• median, sigma, sigma_minus, sigma_plus — median and 16/84 percentile half-widths of the (masked) bootstrap distribution
• covar, correlation — parameter covariance and correlation matrices
• chi2r — reduced χ² of each replicate
• mask — replicates kept after sigma/χ² clipping
• nblocks, mode, granularity — resampling setup
• nfailed — replicates whose fit threw an error

NestedResult

Result of fit_model_nested, whichever sampler produced it.

backend is not decoration. Two independent nested samplers agreeing on logz is evidence; one of them quoted without saying which is not, because the estimator, the bounding strategy and the stopping rule all differ between them.

UltraNestResult

Alias for NestedResult, so code written against the UltraNest-only API keeps working.

FitResult

What every point-estimate fitter returns: fit_model, fit_model_lsqfit and the grid search.

ndof counts data points and is not reduced by the number of free parameters, so it is not a degrees-of-freedom in the model-comparison sense; compute AIC or BIC from chi2 and the parameter count rather than reading them off chi2r.

See also LsqFitResult, UltraNestResult and Chi2Map.

Chi2Map

A χ² surface over two free parameters, with everything needed to draw it and to report the best point on it.

Parameters not on the two axes are held at their starting values throughout: a map is a slice through the χ² surface, not a profile likelihood, and reading it as one overstates how well the two plotted parameters are determined.

See also chi2_map and delta_chi2_levels.

DataBlocks

Partition of an OIdata into resampling units. Each block holds the row indices it owns in each observable table (OI_VIS, OI_VIS2, OI_T3, OI_FLUX). Built by data_blocks.