Design

One package

Everything is in package stats, as everything is in scipy.stats: a notebook imports one name. The distributions, the descriptive statistics, the tests, the kernel density estimate and the bootstrap share the input type Data ([]float64 or a go-ndarray *ndarray.Array, read in row-major order) and the error kinds ErrParam, ErrData, ErrFit.

Distributions

A distribution is a plain struct of its parameters, with SciPy’s names and meaning, and Loc and Scale (continuous) or Loc (discrete):

stats.Beta{A: 2, B: 5, Loc: 0, Scale: 1}.PPF(0.975) // scipy.stats.beta.ppf(0.975, 2, 5)

Its methods are pure functions of the value: safe for concurrent use, nothing cached, nothing allocated on the scalar paths. The zero value has Scale 0 and is invalid: the methods then return NaN, as SciPy does for out-of-domain parameters, and Validate (or the NewXxx constructor) says why with an error wrapping ErrParam.

Standard form plus generated wrapper

Each distribution is written once, in standard form (Loc 0, Scale 1), as an unexported type implementing contStd or discStd (continuous.go, discrete.go). The exported methods are generated (internal/gen, from a //stats:continuous Name scipyname a=A ... directive) and apply, identically for every distribution:

• the parameter check (NaN otherwise);
• z = (x - Loc)/Scale and the Jacobian (pdf/Scale, logpdf - log Scale);
• the boundary rules SciPy applies (density 0 outside the support, cdf 0 and 1 at and beyond the ends, ppf(0) and ppf(1) the support’s ends, NaN outside [0, 1]);
• X = Loc + Scale·Y for quantiles, moments (with SciPy’s convention for moments that do not exist) and draws, the product rounded before the sum as in SciPy (no fused multiply-add: the result is the same on every architecture).

So the standard form only has to be right in the interior, and a new distribution is one file of formulas.

Accuracy techniques

• Double-double where a formula cancels or amplifies: the exponent of a density ((a-1)·log x - x - log Γ(a) loses log2 of its terms’ size) is carried in double-double (internal/dd, a copy of go-numerics/special’s) and exponentiated by an exponential accurate to about 2^-70 with the scale factors multiplied in before the one final rounding.
• No large terms that cancel: for large shapes the gamma and beta densities take Stirling’s factors out analytically (x^(a-1)e^-x/Γ(a) = √(a/2π)/x · exp(-a·g((x-a)/a) - corr(a)), g(u) = u - log(1+u) with its own series near 0), as Boost does.
• Complements, never 1 - p: SF, ISF and the log forms are their own functions (incomplete function complements, log1p(-q) when p is close to 1), and the distribution functions of t and F choose whichever incomplete-beta argument is below 1/2.
• The special functions are go-numerics/special’s, measured at or below SciPy’s error everywhere.

Draws

Rvs(g, size...) takes a go-ndarray *random.Generator, which reproduces NumPy’s Generator bit for bit for the same seed, and calls the same generator method SciPy’s _rvs calls (or the uniform-then-quantile fallback), so that d.Rvs(random.DefaultRNG(42), 1000) is d_scipy.rvs(size=1000, random_state=np.random.default_rng(42)). Which distributions are bit-identical, and where not why, is measured in rvs-table.md.

Fits

FitXxx(x, &FitOptions{Fixed: map[string]float64{"Loc": 0}}) is scipy.stats.xxx.fit(x, floc=0). Where SciPy has a closed form the result is SciPy’s to the bit (same pairwise sums). Where SciPy solves likelihood equations the solution here is polished to full precision. Elsewhere the likelihood is maximised from SciPy’s own starting point by Nelder-Mead restarted to convergence with tight tolerances: the test asserts that our log-likelihood is at least SciPy’s (whose fmin stops at four digits).

Arrays

Scalar methods compose with arrays through stats.Map(x, d.CDF) (parallel over GOMAXPROCS, each element computed once, so the result does not depend on the number of goroutines) or ndarray’s own x.Map(d.CDF).

Verification

Three layers, all in go test:

1. ULP fixtures against mpmath (testdata/fixtures, generated by testdata/gen/gen.py): every distribution function in standard form over regions that cover the bulk and both tails at every shape regime, the exact value stored as a double-double, SciPy’s value at the same arguments alongside. Per-region budgets (testdata/budgets) hold on every CI architecture. See accuracy.md.
2. SciPy parity fixtures (testdata/scipy): parameterised instances (shapes, loc, scale) compared with SciPy in the bulk, moments, entropy, median, interval, support, invalid parameters, fits, and the draws for fixed seeds, bit for bit.
3. Contract tests for every distribution: NaN propagation, boundaries, quantile round trips at float64 resolution, log forms consistent with the plain ones.

The SciPy and NumPy values are produced with the x86-64 builds (under Rosetta on the development host): NumPy’s macOS arm64 wheel rounds some libm calls differently, and x86-64 is the build bit-identity is claimed against.