stats.standard-deviation

def standard_deviation(values: Sequence[float], kind: DeviationKind, decimals: int) -> float

DeviationKind = Literal["population", "sample"]

from fune.stats.standard_deviation import standard_deviation  # stats.standard-deviation@^1

import math
from typing import Sequence

from .stats_standard_deviation_types import DeviationKind


POW10 = [1.0, 10.0, 100.0, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9, 1e10, 1e11, 1e12]


def _round_to(x: float, decimals: int) -> float:
    # Half away from zero on the binary64 value. Python's round() is
    # half-even and would give round(0.5) == 0; this does not.
    scale = POW10[decimals]
    y = x * scale
    r = float(math.floor(y))
    if y - r >= 0.5:
        r += 1.0
    return r / scale + 0.0


def standard_deviation(values: Sequence[float], kind: DeviationKind, decimals: int) -> float:
    """Population or sample standard deviation, two-pass.

    The one-pass "mean of squares minus square of mean" shortcut cancels
    catastrophically on large, close values; two passes do not.
    """
    if isinstance(values, (str, bytes)) or not isinstance(values, (list, tuple)):
        raise TypeError("values must be a list of numbers")
    for v in values:
        if isinstance(v, bool) or not isinstance(v, (int, float)) or not math.isfinite(v):
            raise TypeError("values must be finite numbers, received %r" % (v,))
    if kind not in ("population", "sample"):
        raise ValueError('unknown standard deviation kind "%s"' % (kind,))
    if isinstance(decimals, bool) or not isinstance(decimals, int) or decimals < 0 or decimals > 12:
        raise ValueError("decimals must be a whole number from 0 to 12, received %r" % (decimals,))
    n = len(values)
    if n == 0:
        raise ValueError("values must not be empty")
    if kind == "sample" and n < 2:
        raise ValueError("sample standard deviation needs at least 2 values")

    # Plain left-to-right float additions; sum() on ints would be exact and
    # math.fsum compensated, and either would differ from the other languages.
    total = 0.0
    for v in values:
        total = total + float(v)
    mean = total / n
    squares = 0.0
    for v in values:
        d = float(v) - mean
        squares = squares + d * d
    divisor = n - 1 if kind == "sample" else n
    return _round_to(math.sqrt(squares / divisor), decimals)

Install

fune build

With that line in your source, in a Python project (language python in fune.project), fune build resolves it and nothing else, pins them in fune.lock, downloads only the Python package of each, and builds the code above into your project’s .fune/build, one readable file per capability with a header linking back here. Or pin a range in fune.project and build in one step:

fune add stats.standard-deviation

The manifest, vectors and README with only the Python implementation. Install it without the registry with fune add ./stats.standard-deviation-1.0.0-python.fune, or fetch it from a terminal with fune pull stats.standard-deviation@1.0.0:python.

The whole function, every language, is one file too: stats.standard-deviation-1.0.0.fune, 12,588 bytes, sha256 aea4bddbb12330d26f7bdff755ca24697691ed9aece2b97892c034b13986048d. It installs into a project of any language.

# fune: before stats.standard-deviation

# fune: after stats.standard-deviation

# fune: step stats.standard-deviation after <n|label>

More from the author

**Algorithm.** Two passes: the mean first, then the sum of squared deviations from it. The one-pass textbook shortcut, mean of squares minus square of the mean, cancels catastrophically when the values are large and close together (four readings near one billion come out as garbage or even a negative variance); the vectors include that case.

**Precision.** The result is rounded to `decimals` places, half away from zero, applied to the binary64 value: y = x x 10^decimals, r = floor(y), plus one if y - r >= 0.5, divided back by 10^decimals. So a deviation of exactly 0.5 rounds to 1 at zero places, where Python's `round()` gives 0.

**Why the three languages agree to the bit.** Sums run left to right, and only IEEE-754 +, -, x, /, square root and floor are used. All of these are correctly rounded by the standard (square root included, unlike sin or exp, whose last bit varies between math libraries), so TypeScript, Python and Rust hold the same double before rounding and return the same rounded value.

Sources: NIST/SEMATECH e-Handbook of Statistical Methods, section 1.3.5.6 "Measures of Scale"; B. P. Welford, "Note on a Method for Calculating Corrected Sums of Squares and Products", Technometrics 4(3), 1962, on why the one-pass formula fails.

Files