fstats

FSTATS is a modern Fortran statistical library. The public API is collected in the fstats module, so applications can generally start with use fstats.

Status

Feature Matrix

Category	Functionality	Status	Documentation	Examples
Descriptive statistics	Means, variance, standard deviation, medians, quantiles, trimmed means, covariance, and pooled variance	✅	✅	N/A
Probability distributions	Normal, log-normal, Student's t, F, chi-squared, binomial, Poisson, and multivariate normal distributions	✅	✅	Distribution
Statistical inference	Confidence intervals, t-tests, F-tests, Bartlett's test, Levene's test, and sample-size calculations	✅	✅	ANOVA, model comparison
ANOVA	One-factor and two-factor analysis of variance	✅	✅	ANOVA
Regression	Polynomial linear least squares, regression statistics, R-squared metrics, correlations, numerical Jacobians, and nonlinear Levenberg-Marquardt least squares	✅	✅	Linear, nonlinear, GLM
DOE design and analysis	Full and fractional factorial designs (coded on [-1, 1]), central-composite designs, Latin hypercube sampling, polynomial model fitting (linear, 2-way and 3-way with interaction and squared terms), diagnostics, design efficiency, and confidence/prediction intervals using the t-distribution	✅	✅	Design efficiency, prediction, full factorial
DOE model comparison and ANOVA	Model comparison and model ANOVA	✅	✅	Model comparison
Response-surface optimization	Maximize, minimize, or target a fitted response surface within factor bounds; analytic stationary-point solution (with maximum/minimum/saddle classification) for quadratic models, and a multi-start projected-gradient search otherwise. Results are local optima	✅	✅	Example
Measurement systems analysis	Crossed, nested, and expanded gauge R&R studies with variance components, contribution, study variation, tolerance, and distinct categories	✅	✅	N/A
Sampling and resampling	Bootstrap, random and rejection sampling, Box-Muller, and multivariate-normal sampling	✅	✅	Bootstrap, rejection, Box-Muller
Markov chain Monte Carlo	Chains, target distributions, proposals, samplers, and model evaluation	✅	✅	Regression, example 2, example 3
Signal and numerical methods	Allan variance, LOWESS smoothing, and linear, polynomial, spline, and Hermite interpolation	✅	✅	Allan, interpolation, LOWESS
Missing-data handling	Missing-data imputation	✅	✅	N/A
Special functions	Beta and gamma functions, regularized and incomplete forms, and the digamma function	✅	✅	N/A
Data containers and I/O	Typed data columns and tables with delimited-file input and output	✅	✅	N/A

Documentation links to the Fortran API reference. Status: ✅ available; ⚠️ partial implementation. N/A indicates no standalone example.

Documentation

The generated API documentation is available here.

Building with CMake

CMake 3.24 or newer, a Fortran 2018 compiler, Git, and an OpenMP-capable Fortran compiler/runtime are required. CMake searches for compatible installations of LINALG, COLLECTIONS, and FSTRING; if they are not found, it fetches them automatically. See External Libraries for details.

Configure and build the library:

cmake -S . -B build -DCMAKE_BUILD_TYPE=Release
cmake --build build --parallel

Tests and examples are disabled by default. Enable either or both at configure time:

cmake -S . -B build -DBUILD_TESTING=ON -DBUILD_FSTATS_EXAMPLES=ON
cmake --build build --parallel
ctest --test-dir build --output-on-failure

Install the library, module files, CMake package files, and pkg-config metadata with:

cmake --install build --prefix ./install

FPM can also be used to build this library using the provided fpm.toml.

fpm build --profile release
fpm test --profile release

FPM resolves the linalg, collections, fstring, and openmp dependencies declared by this project. The test suite also uses the fortran_test_helper development dependency. To use FSTATS as a dependency in another FPM project, add the following to that project's fpm.toml:

[dependencies]
fstats = { git = "https://github.com/jchristopherson/fstats", tag = "v1.7.0" }

Then import the API in Fortran with use fstats and build normally with fpm build.

Parallel MCMC

sample_chains runs independent Metropolis-Hastings chains with OpenMP. Supply equally sized arrays of samplers, proposals, and targets. Each target must be initialized independently, rather than copied from an object whose collections or pointer components might share mutable state. The observations are shared read-only; models and sampler callbacks must also be thread-safe.

For a user-defined custom_target extending mcmc_target, with xdata and ydata already defined, declare and initialize the chain objects as follows:

type(mcmc_sampler) :: samplers(4)
type(mcmc_proposal) :: proposals(4)
type(custom_target) :: targets(4)
type(normal_distribution) :: prior
integer :: chain_index

prior%mean_value = 0.0d0
prior%standard_deviation = 1.0d0
do chain_index = 1, size(targets)
	! Add each parameter's prior independently for this chain.
	call targets(chain_index)%add_parameter(prior)
	call targets(chain_index)%add_parameter(prior)
end do
call sample_chains(samplers, xdata, ydata, proposals, targets, niter=10000)
! Retrieve each result with samplers(chain_index)%get_chain().

The parameter count must match the custom model. Each proposal adapts its own scale. As with sample, repeated calls append samples and initialize a new run from the prior means; they do not resume from the last stored state. Omitting niter uses 10,000 iterations per chain. Arrays of polymorphic objects must have a common dynamic type within each array.

Likelihood evaluation uses scaled, two-pass SIMD residual reductions without temporary residual or log-density arrays. Outside an existing OpenMP parallel region, it uses worker threads when the observation count reaches the target's public likelihood_parallel_threshold (default 10,000). This default is a starting point, not a benchmark-derived optimum: tune it for the workload, set it to zero to request threading at every size, or to huge(1_int32) to disable threading for ordinary dataset sizes. Nested worker teams are suppressed so independent chains do not also launch likelihood teams. Reduction order can change floating point results and, consequently, sampled trajectories.

Set OMP_NUM_THREADS to control concurrency, and avoid oversubscription from threaded BLAS or custom models. Parallel chains require a compiler runtime with thread-safe intrinsic random_number, such as GNU Fortran. Random-stream assignment and reproducibility depend on that runtime and thread scheduling; this API does not provide per-chain seeds or thread-count-independent results. Without OpenMP enabled, the same API runs serially.

Distribution Modules

Each concrete distribution has its own module and source file in src/distributions. The module name is fstats_ followed by the type name, for example:

use fstats_distributions, only : distribution
use fstats_normal_distribution, only : normal_distribution

fstats_distributions contains the abstract distribution and multivariate_distribution types, their interfaces, and shared supporting routines. It no longer exports concrete distribution types. Replace direct imports of concrete types from that module with their individual modules, or continue to use use fstats to access all distributions unchanged.

Each built-in distribution implements its own log-density and tail bindings; the abstract bases supply generic compatibility fallbacks for custom laws.

Numerical Behavior

Distributions expose log_pdf, survival, log_cdf, and log_survival in addition to pdf and cdf. Prefer log densities for posterior calculations and direct survival functions for small upper-tail probabilities. Built-in laws use analytic log densities and direct beta/gamma or normal tails. Discrete PDFs are probability masses at integer-valued points; their CDFs include all masses at or below the supplied point, and survival means P(X > x).

Density and tail evaluations with invalid parameters or NaN inputs return NaN. Outside support, PDFs return zero and log densities return negative infinity. A support-endpoint density singularity can legitimately return positive infinity. Custom distributions have compatibility fallbacks that take logs of their PDFs or CDFs; override these methods for accurate extreme tails. Subclasses that change a built-in probability law must also override its log-density/tail methods.

log_regularized_beta, log_regularized_gamma_lower, and log_regularized_gamma_upper avoid underflow in special-function probabilities. The regularized routines accept an optional max_iterations limit and return NaN on invalid input or failure to converge. stable_log1p and stable_expm1 retain small differences lost by log(1 + x) and exp(x) - 1.

MCMC stores log variance and evaluates the default Gaussian likelihood without exponentiating it. likelihood still accepts variance, but nonpositive or nonfinite variance now returns NaN rather than being silently clamped. likelihood_log_variance accepts log variance directly. log_posterior uses evaluate_log_variance_prior, whose density includes the variance-transform Jacobian. Custom targets changing the likelihood or variance prior must override these log-coordinate bindings; overriding only the variance-space methods no longer changes sampling. Initial states must have finite log posterior; invalid or impossible proposals are rejected. Correcting the Jacobian changes sampled variance distributions relative to earlier versions.

Descriptive statistics use scaled/compensated reductions. Empty statistics, sample moments with fewer than two observations, and nonfinite observations return NaN; constant-data correlation also returns NaN. Quantiles require 0 <= q <= 1, trimming requires 0 <= p < 0.5, and pooled variance requires at least two observations per group. These contracts replace some earlier zero results or silent clamping. Standard deviation can remain finite even when variance overflows; truly unrepresentable moments return infinity.

Multivariate-normal log densities use Cholesky solves and log determinants. Regression covariance uses an SVD of the design matrix, not its normal equations. regression_covariance optionally reports numerical rank and condition number; nonlinear fits expose the same diagnostics in convergence_info when covariance is requested. Rank-deficient covariance is a pseudocovariance: discarded directions are not evidence of zero parameter uncertainty. Nonlinear covariance requires additional evaluations of the final Jacobian, counted in the reported function evaluations. Ordinary regression coefficients retain their existing parameterization; scaling polynomial predictors remains advisable.

Experimental Design Modules

All experimental design types are declared in fstats_experimental_design. Procedures are organized into the following modules in src/design:

• fstats_doe_designs: design generation, size calculations, and efficiency assessment.
• fstats_doe_coding: natural-to-coded and coded-to-natural factor conversions.
• fstats_doe_models: model fitting and evaluation.
• fstats_doe_diagnostics: goodness of fit, residual analysis, model comparison, and model ANOVA.
• fstats_doe_prediction: predictions and uncertainty intervals.
• fstats_doe_response_surface: response-surface optimization.

use fstats continues to expose all these types and procedures. Code that previously imported procedures directly from fstats_experimental_design must instead import them from the corresponding module above, or use fstats. Import shared types from fstats_experimental_design when using individual procedure modules. This organization does not change numerical behavior; existing approximations and placeholder calculations remain unchanged.

External Libraries

FSTATS has the following direct dependencies:

Dependency	Purpose	Build integration
LINALG	Linear algebra operations	CMake uses an installed compatible package or fetches it; FPM resolves it from fpm.toml.
COLLECTIONS	Collection types used by the library	CMake uses an installed compatible package or fetches it; FPM resolves it from fpm.toml.
FSTRING	String support	CMake uses an installed compatible package or fetches it; FPM resolves it from fpm.toml.
OpenMP	Fortran parallel runtime, including parallel MCMC	Required by CMake and declared as an FPM dependency.

The test suite has an additional development dependency, fortran_test_helper, which CMake can also fetch when tests are enabled. An optimized BLAS/LAPACK implementation is recommended for performance with LINALG; it is not a direct FSTATS dependency.