fit Function

fit fits steady-state or transient GM/GRSM models to data for a single gene (or coupled units), writes results through finalize, and returns fit results and diagnostics.

For coupled transcribing units, arguments transitions, G, R, S, insertstep, and trace-related settings become tuples of the single-unit type (e.g. two units with G = (2, 3)).

Syntax

fits, stats, measures, data, model, options = fit(; kwargs...)

The keyword form is the usual entry point. Positional overloads exist for advanced / legacy callers; optional inference keywords on the positional path are forwarded into make_structures (see the _MAKE_STRUCTURES_OPTION_KW constant in fit.jl).

Inference methods

Posterior / variational inference is selected with inference_method (also accepted: plain symbols :mh, :nuts, :advi, or the aliases INFERENCE_MH, INFERENCE_NUTS, INFERENCE_ADVI):

Shared budget keywords are harmonized in load_options when building option structs:

• samplesteps: MH — posterior samples to collect; NUTS — n_samples; ADVI — maxiter unless you set maxiter explicitly in the run dict.
• warmupsteps: MH — discarded warmup (shares wall time with sampling via maxtime); NUTS — n_adapts unless n_adapts is set explicitly; ADVI — not the same semantic object (use n_mc for ELBO Monte Carlo draws).

Other cross-method options stored on MHOptions, NUTSOptions, and ADVIOptions in common.jl:

• device: :cpu or :gpu (GPU paths may error if unsupported for a method).
• parallel (alias parallelism): :single, :threaded, or :distributed — used with nchains for multi-chain NUTS/ADVI dispatch (run_inference); MH multi-chain still uses the existing distributed MH chain runner when nchains > 1.
• gradient: method-specific; e.g. :finite, :ForwardDiff, :Zygote for NUTS/ADVI; :none default for MH. String values in TOML/JSON-style dicts are coerced when possible.

For NUTS, the package default is gradient = :finite. This uses central finite differences on the optimized primal likelihood and is the recommended production path for trace, dwell-time, coupled, and other nontrivial likelihoods unless a local benchmark shows :ForwardDiff or :Zygote wins for that specific model. :ForwardDiff and :Zygote remain useful opt-in validation/benchmark paths.

For ADVI, gradient = :finite and gradient = :ForwardDiff both differentiate the ELBO with ForwardDiff. Trace ADVI with plain gradient = :Zygote intentionally falls back to that safer path unless zygote_trace = true.

make_structures merges _current_run_spec[] (from fit(; key=...)) with explicit samplesteps / warmupsteps / maxtime / temp and any extra kwargs..., then calls load_options on that dict. You can also call load_options directly in scripts or tests.

Arguments (high level)

Basic model parameters

• G::Int, R::Int, S::Int, insertstep::Int, transitions: Model topology (see in-source fit docstring).
• coupling, grid, hierarchical: Coupled / grid / hierarchical layouts.
• Shared-parameter multi-fit: hierarchical may also be a keyed shared-parameter named tuple with kind=:shared, levels::HierarchySpec, datasets::Vector{DatasetSpec}, rate_families, transition_families, and emission_families. This combines datasets from multiple folders while sharing some parameter families exactly and splitting others. It does not introduce learned hyperparameters. See Shared Parameter Trace Fits. Shared fits write combined informational files for the run, but operational rates_*.txt files are split by shared-rate group and keep ordinary rate labels.

Data parameters

• datatype: String or symbol for a single legacy datatype, e.g. "rna", :trace, "rnadwelltime", "tracejoint". A tuple or vector such as (:rna, :dwelltime) requests the CombinedData path.
• datapath: Path, vector/tuple of paths, or for CombinedData preferably a modality-keyed NamedTuple, e.g. (rna = "smFISH", dwelltime = ["ON.csv", "OFF.csv"]).
• datacond, cell, gene, nalleles, trace_specs, dwell_specs, …

Fitting / inference parameters

• nchains::Int: Number of parallel chains for fit(nchains, data, model, …) dispatch (run_inference); for MH this matches the existing pooled-chain behavior; for NUTS/ADVI, multi-chain runs are merged when nchains > 1 (see src/inference_common.jl in the package source).
• maxtime: Primary wall budget for MH (warmup + sampling); numeric seconds or strings like "90m", "2h" (see maxtime_seconds). Interpretation for NUTS/ADVI is method-specific (many NUTS/ADVI controls are in the option structs from load_options).
• samplesteps, warmupsteps, temp: As above; temp is MH temperature; NUTS/ADVI paths use a neutral value for finalize when not applicable.
• propcv: MH proposal CV / covariance loading (see Package overview).

Priors, indices, outputs

• priormean, priorcv, noisepriors, fittedparam, fixedeffects, onstates, decayrate, …
• resultfolder, label, writesamples, burst, optimize, …

For RNA histograms, rna_truncation=:legacy reproduces the v0.7.8 histogram support rule; use :none for all bins. Keep this setting fixed when comparing runs. A negative decayrate or negative decay placeholder in priormean is resolved from gene/cell halflife metadata and converted to inverse minutes before initialization, with 1.0 as the unavailable-metadata fallback.

burst=true computes a posterior burst-size distribution during finalization. optimize=true runs a separate LBFGS optimization initialized at the best sampled point. Neither analysis is recreated later by dataframe assembly.

Do not confuse these summaries:

• row 1 of rates_*.txt is the largest-likelihood point encountered during sampling;
• posterior mean and median are rows 2 and 3;
• optimized_*.txt contains the separate numerical optimizer result;
• burst_*.txt summarizes burst size evaluated over posterior samples.

v2.0 API changes

• CombinedData: Multimodal fits use datatype = (:rna, :trace) or datatype = (:rna, :dwelltime). The order is canonicalized, likelihoods are evaluated per modality, scalar likelihoods are summed, and WAIC pointwise predictions are concatenated in canonical modality order.
• Retired legacy input keywords: infolder and inlabel are no longer part of the public API. Use root, datapath, label, and resultfolder. Older key-based run specs may still contain infolder; fit(; key=...) ignores it during migration.
• Trace and dwell metadata: Prefer trace_specs and dwell_specs. Legacy traceinfo and dttype entries may be consumed from old run specs, then are dropped from newly written run specs.

Run specification and key-based naming

• key = nothing: When a string, fit loads info_<key>.jld2 (companion to the marker TOML) and merges with explicit keywords (keywords win). Outputs use that stem. See Run specification (info TOML).

Returns

The top-level fit(; …) return tuple is:

fits, stats, measures, data, model, options = fit(; ...)

• fits: Fit — posterior samples (MH/NUTS) or variational mean column (ADVI); ll, WAIC-related fields, acceptance summaries depend on the method.
• stats, measures: Stats, Measures — comparable structs across methods where meaningful (ADVI may use single-point proxies for some diagnostics; see the run_advi_fit docstring in the source).
• data, model, options: The concrete Options subtype is MHOptions, NUTSOptions, or ADVIOptions depending on inference_method.

Examples

Basic RNA histogram (MH default)

fits, stats, measures, data, model, options = fit(
    G = 2,
    R = 0,
    transitions = ([1,2], [2,1]),
    datatype = "rna",
    datapath = "data/HCT116_testdata/",
    gene = "MYC",
    datacond = "MOCK",
)

NUTS (same keyword surface; different options type)

fits, stats, measures, data, model, options = fit(
    G = 2,
    R = 0,
    transitions = ([1,2], [2,1]),
    datatype = "rna",
    datapath = "data/HCT116_testdata/",
    gene = "MYC",
    datacond = "MOCK",
    inference_method = :nuts,
    samplesteps = 500,
    warmupsteps = 250,
    parallel = :single,
    gradient = :finite,
)

Omitting gradient in a NUTS fit is equivalent to gradient = :finite.

Combined RNA + dwell-time data

fits, stats, measures, data, model, options = fit(
    datatype = (:rna, :dwelltime),
    datapath = (
        rna = "HBEC_smFISH",
        dwelltime = [
            "dwelltime/CANX_ON.csv",
            "dwelltime/CANX_OFF.csv",
            "dwelltime/CANX_ONG.csv",
        ],
    ),
    gene = "CANX",
    cell = "HBEC",
    datacond = "",
    resultfolder = "FISH",
    dwell_specs = [
        (
            unit = 1,
            onstates = [Int[], Int[], [2, 3]],
            dttype = ["ON", "OFF", "ONG"],
        ),
    ],
)

For new scripts, the keyed datapath form is preferred over the legacy positional "rnadwelltime" layout. See CombinedData API.

Trace fit with multiple MH chains

fits, stats, measures, data, model, options = fit(
    G = 3,
    R = 2,
    S = 2,
    insertstep = 1,
    transitions = ([1,2], [2,1], [2,3], [3,1]),
    datatype = "trace",
    datapath = "data/testtraces",
    cell = "TEST",
    gene = "test",
    datacond = "testtrace",
    trace_specs = [(unit = 1, interval = 1.0, start = 1.0, t_end = -1.0, zeromedian = true, active_fraction = 1.0, background = 0.0)],
    noisepriors = [40., 20., 200., 10.],
    nchains = 4,
)

Trace observation likelihoods use Gaussian densities (prob_Gaussian) for the observed fluorescence values. Because a continuous density has units of inverse fluorescence and can be greater than 1 when the noise scale is small, trace log-likelihood contributions can be positive. This is expected and does not imply a probability greater than 1; it only means the density per fluorescence unit is high near the observed value. If a future workflow needs likelihoods interpreted as finite-bin probability masses, a measurement-resolution correction can be applied after the HMM likelihood, e.g. logL + n_observations * log(delta_y), with delta_y = 1.0 preserving the current behavior.

ExGaussian Background

probfn=prob_ExGaussian opts into Gaussian reporter noise plus a right-skewed exGaussian background. Existing Gaussian and mixture models remain unchanged. Supply five noisepriors entries: [background_mean, background_gaussian_sd, reporter_mean, reporter_sd, tail_scale]. An independent exponential contribution with mean tail_scale is added and centered by subtracting that same mean. Thus backgroundmean retains its meaning. For reporter count n, the mean is `backgroundmean+nreportermeanand variance isbackgroundgaussian_sd^2+nreportersd^2+tailscale^2`. Do not use total background SD as the Gaussian component SD or the tail variance is counted twice. The model changes the marginal distribution, not temporal noise correlations.

For a 3301 model with nine kinetic rates, the five noise slots are 10:14. To hold a calibrated tail fixed, omit slot 14 from both fittedparam and the individual-noise indices in hierarchical; for example, fittedparam=vcat(1:8,10:13), hierarchical=(2,collect(10:13),()) fits kinetics (except decay) and the original four noise quantities, with tailscale fixed. Zero tailscale gives exactly Gaussian emissions. When fitting tailscale, initialize it strictly positive: its sampling transform is logarithmic. Use complete rate vectors of the new size: do not append parameters to an old hierarchical file without rebuilding each complete block. Correlation contexts and coupled `simulatetracevector/simulatetrialsinfer these five noise parameters when counts are omitted, including mixed Gaussian/exGaussian units. Saved per-unit noise priors take precedence in the correlation file writers. Low-levelsimulatorcalls still require the appropriate explicitnoiseparams` count for non-Gaussian emission families.

Quick-start helpers

For steady-state RNA histograms, fit_rna wraps fit with defaults:

result = fit_rna(rna_histogram(counts); G=2, name="MYC", halflife=600.0)
describe_fit(result)

• fit_rna(data; G, R, S, insertstep, name, halflife, outdir, samplesteps, nchains, kwargs...): data is a histogram vector or a file path. Inputs are staged under outdir/data; results go to outdir/results/<name>. halflife is in minutes; if omitted the decay rate is 1.0 per minute. Extra keywords pass to fit. Returns the fit outputs plus rates, outdir and resultdir.
• rna_histogram(counts): per-cell RNA counts to a histogram.
• fit_traces(data; G, R, S, insertstep, name, interval, halflife, noisepriors, zeromedian, datacol, outdir, samplesteps, warmupsteps, nchains, maxtime, kwargs...): live-cell traces from a folder of files, a vector of file paths, or a vector of intensity vectors. interval is the frame interval in minutes and R must be at least 1. Files are copied to outdir/data/traces. warmupsteps defaults to samplesteps (40,000), because trace fits mix poorly without a long adaptive warmup. With zeromedian = false (default) the raw intensities are used and noisepriors default to an estimate from the data; with zeromedian = true traces are median-subtracted and scaled and noisepriors default to [0.0, 0.1, 1.0, 0.1], which biases the fit if the real noise differs. Returns the same fields as fit_rna plus ntraces.
• fit_summary(resultdir): DataFrame of ML, mean and median rates (per minute).
• describe_fit(result): prints the table and output location.
• demo(): short fit of the bundled HCT116_testdata/MYC_MOCK.txt.
• demo_traces(): simulates traces from a known model, fits them and prints fitted and true rates. A short run is only a rough check of the workflow.

See Choose a model and data for model choice, units and file formats.

Notes

1. Rate order — G transitions, R transitions, S if present, decay, noise parameters (see package overview).

2. Convergence — R-hat, ESS, etc. apply to sample-based chains; interpret ADVI diagnostics separately.

3. MH proposal covariance and warmup — See Package overview.

4. Key-based workflows — Use key="..." for reproducible cluster runs and stage_write_run_specs; include inference_method / parallel / gradient in the saved dict when needed.

5. Cluster scripts — makeswarm and related helpers can emit fit(; key=..., inference_method=:nuts, …) overrides; positional gene/coupled scripts append ; kw=... suffixes for the same keywords. See Cluster and batch workflows.

6. Migration from older inputs — Replace infolder with datapath; replace inlabel usage with label; keep output routing in resultfolder. Prefer trace_specs / dwell_specs over traceinfo / dttype.