Skip to content

Distribution Fitting

Chapter summary

This chapter explains how ForecastAI fits parametric demand distributions to forecast quantile output, why each family is chosen, how the fitting algorithm works, and how the fitted parameters flow downstream into the MEIO safety-stock optimizer.

Why fit distributions?

Every inventory policy that goes beyond a fixed reorder point requires a model of demand uncertainty over the replenishment lead time. A point forecast gives the expected demand; the fitted distribution gives the full probability landscape — the probability that demand will be 20 % above forecast, or 50 % above, or anywhere on the tail.

The MEIO optimizer uses the fitted distribution to evaluate, for any candidate committed buffer level \(b\), the expected fraction of demand that will be satisfied from stock — the fill rate integral:

\[ FR(b) = 1 - \frac{E\!\left[\max(D_{LT} - b,\; 0)\right]}{E[D_{LT}]} \]

where \(D_{LT}\) is the total demand over lead time. Without a distribution, this integral cannot be computed; MEIO degenerates to a heuristic.

Connection to forecast quantiles

ForecastAI's statistical and neural forecasting models (AutoARIMA, ETS, NHITS, PatchTST …) produce quantile forecasts at levels \(\{0.10, 0.25, 0.50, 0.75, 0.90, 0.95, 0.99\}\) for each horizon step. Distribution fitting converts these discrete quantile points into a continuous parametric family, enabling:

  • Evaluation of the fill-rate integral at any buffer level (not just the pre-computed quantile levels).
  • Scaling to arbitrary lead times via moment-based extrapolation.
  • Compact storage — six numbers replace a 99-point table.

The six distribution families

Normal (Gaussian)

\[ f(x;\,\mu,\sigma) = \frac{1}{\sigma\sqrt{2\pi}}\, \exp\!\left(-\frac{(x-\mu)^2}{2\sigma^2}\right) \]

Parameters: \(\mu\) (location = mean), \(\sigma\) (scale = standard deviation).

Domain: \((-\infty, +\infty)\), but in practice truncated at zero for demand.

When appropriate: High-volume, symmetric demand; the Central Limit Theorem makes Normal an excellent approximation for aggregated weekly demand across many customers when the coefficient of variation \(CV = \sigma/\mu < 0.5\).

Fitting method (quantile matching):

  1. If the median quantile \(Q_{0.5}\) is available, set \(\hat{\mu} = Q_{0.5}\); otherwise use the mean of all available quantile values.
  2. Estimate \(\hat{\sigma}\) from the interquartile range when both \(Q_{0.25}\) and \(Q_{0.75}\) are available:
\[ \hat{\sigma} = \frac{Q_{0.75} - Q_{0.25}}{1.349} \]

The divisor 1.349 follows from the standard normal identity \(\Phi^{-1}(0.75) - \Phi^{-1}(0.25) = 1.349\). When the IQR quartiles are absent, the spread of the available extreme quantiles is used:

\[ \hat{\sigma} = \frac{Q_{p_{\max}} - Q_{p_{\min}}} {\Phi^{-1}(p_{\max}) - \Phi^{-1}(p_{\min})} \]
  1. Poisson variance floor. For discrete demand with low volume, quantiles from probabilistic forecasts are often nearly identical, which would produce \(\hat{\sigma} \approx 0\). The implementation enforces:
\[ \hat{\sigma} \geq \sqrt{\max(\hat{\mu},\, 0)} \]

This matches the natural lower bound on variability for Poisson-distributed counts (where \(\sigma^2 = \mu\)).

Fill-rate integral (Normal):

\[ FR(b) = 1 - \frac{\sigma_{LT}\cdot G\!\left(\tfrac{b-\mu_{LT}}{\sigma_{LT}}\right)}{\mu_{LT}} \]

where \(G(z) = \phi(z) - z\,[1 - \Phi(z)]\) is the standard normal loss function, \(\phi\) is the standard normal PDF, and \(\Phi\) is the CDF.

Gamma

\[ f(x;\,\alpha,\beta) = \frac{x^{\alpha-1}\,e^{-x/\beta}}{\beta^\alpha\,\Gamma(\alpha)}, \quad x > 0 \]

Parameters: \(\alpha > 0\) (shape), \(\beta > 0\) (scale).

Moments: \(\mu = \alpha\beta\), \(\sigma^2 = \alpha\beta^2\), whence \(\alpha = \mu^2/\sigma^2\) and \(\beta = \sigma^2/\mu\).

Domain: \((0, +\infty)\) — naturally suited to demand (always non-negative).

When appropriate: Positive, right-skewed demand typical of intermittent spare parts with non-zero demand; \(CV \in (0.5, 2)\).

Fitting method (quantile matching via L-BFGS-B):

\[ \min_{\alpha,\beta} \sum_{i}\!\left[Q_{\alpha,\beta}(p_i) - \hat{q}_i\right]^2 \]

where \(Q_{\alpha,\beta}(p)\) is the Gamma CDF inverse (percent point function). Initial guesses are method-of-moments estimates:

\[ \alpha_0 = \frac{\bar{q}^2}{\widehat{Var}(q)}, \quad \beta_0 = \frac{\widehat{Var}(q)}{\bar{q}} \]

Both parameters are constrained to \((0.01, \infty)\) during optimization.

Negative Binomial

\[ P(X = k;\,r,p) = \binom{k+r-1}{k}\,p^r\,(1-p)^k, \quad k = 0,1,2,\ldots \]

Parameters: \(r > 0\) (number of successes / dispersion), \(p \in (0,1)\) (success probability).

Moments: $\(\mu = \frac{r(1-p)}{p}, \quad \sigma^2 = \frac{r(1-p)}{p^2}\)$

Note \(\sigma^2 > \mu\), so the Negative Binomial models overdispersed count demand — variance exceeds the mean. This makes it the go-to family for intermittent demand with lumpy (bursty) orders.

When appropriate: Discrete demand where variance > mean, e.g. spare parts ordered in batches or demand that varies by customer segment.

Fitting method: L-BFGS-B quantile matching with bounds \(r \in (0.01, \infty)\), \(p \in (0.01, 0.99)\). When the sample exhibits overdispersion (\(\widehat{Var} > \bar{q}\)), the method-of-moments initial guesses are:

\[ p_0 = \frac{\bar{q}}{\widehat{Var}}, \quad r_0 = \frac{\bar{q}\,p_0}{1-p_0} \]

Otherwise the algorithm starts from \(p_0 = 0.5\), \(r_0 = \bar{q}\).

Log-Normal

\[ f(x;\,\mu,\sigma) = \frac{1}{x\sigma\sqrt{2\pi}}\, \exp\!\left(-\frac{(\ln x - \mu)^2}{2\sigma^2}\right), \quad x > 0 \]

Parameters: \(\mu\) (mean of \(\ln X\)), \(\sigma\) (std of \(\ln X\)).

Moments of the original scale: $\(E[X] = e^{\mu + \sigma^2/2}, \quad Var(X) = \left(e^{\sigma^2}-1\right)e^{2\mu+\sigma^2}\)$

When appropriate: Demand driven by multiplicative shocks (a product that sells 1 unit most weeks but occasionally sells 50 due to a large order). Heavy right tail; \(CV > 1\).

Fitting method: Log-space moment estimation as initial guess:

\[ \hat{\mu}_{\ln} = \overline{\ln q_i}, \quad \hat{\sigma}_{\ln} = s(\ln q_i) \]

then L-BFGS-B optimizes \((s, \text{scale})\) in the scipy parameterisation \((\text{scale} = e^{\hat{\mu}_{\ln}})\) against the quantile residuals. A small offset of \(10^{-10}\) is added before taking logarithms to handle zero quantile values.

Weibull

\[ f(x;\,k,\lambda) = \frac{k}{\lambda}\left(\frac{x}{\lambda}\right)^{k-1} \exp\!\left(-\left(\frac{x}{\lambda}\right)^k\right), \quad x \geq 0 \]

Parameters: \(k > 0\) (shape), \(\lambda > 0\) (scale).

Moments: $\(\mu = \lambda\,\Gamma\!\left(1+\tfrac{1}{k}\right), \quad \sigma^2 = \lambda^2\!\left[\Gamma\!\left(1+\tfrac{2}{k}\right) - \Gamma\!\left(1+\tfrac{1}{k}\right)^2\right]\)$

When appropriate: Failure-time driven demand, MRO (Maintenance Repair and Overhaul) parts whose consumption is driven by asset operating hours rather than customer pull. When \(k < 1\) the hazard rate decreases (infant mortality); \(k = 1\) gives an exponential (constant hazard); \(k > 1\) gives increasing hazard (wear-out). Aerospace and industrial MRO datasets show a preponderance of Weibull-fitted series in the memory notes for this project.

Fitting method: scipy weibull_min parameterisation with \(\text{loc} = 0\) fixed. L-BFGS-B with bounds \(c \in (0.05, 50)\), \(\lambda \in (0.01, \infty)\), initial guess \(c_0 = 1.5\), \(\lambda_0 = \text{median}(q_i)\).

Poisson

\[ P(X = k;\,\mu) = \frac{\mu^k e^{-\mu}}{k!}, \quad k = 0,1,2,\ldots \]

Parameters: \(\mu > 0\) (rate = mean = variance).

When appropriate: Slow-moving, intermittent discrete demand where \(\sigma^2 \approx \mu\) — the Equi-dispersion regime. Typical for very slow spare parts or items with infrequent demand (ADI > 1.32 but without the lumpiness that would warrant Negative Binomial).

Fitting method: Single-parameter L-BFGS-B minimization:

\[ \min_{\mu > 0}\; \sum_i \left[\text{ppf}_{\text{Poi}}(\mu, p_i) - \hat{q}_i\right]^2 \]

with \(\mu_0 = \bar{q}\).

Quantile matching — the fitting algorithm

All six families share the same fitting strategy: quantile matching (also called the method of percentiles or minimum distance estimation).

Figure: Distribution-fitting decision flow — quantiles are averaged across the horizon, then all six families are fit and scored; sparse data falls back to Normal.

flowchart TD
    Q["Forecast quantiles<br/>Q10..Q99 per horizon step"] --> AVG["Average across horizon<br/>→ single-period quantiles q_p"]
    AVG --> NS{"n observations<br/>>= distribution_threshold (25)?"}
    NS -- no --> NF["Fallback: Normal<br/>mu = Q0.5, sigma = max(sqrt(mu), 1e-3)"]
    NS -- yes --> FIT["Fit all 6 families<br/>Normal, Gamma, NegBin, LogNormal, Weibull, Poisson"]
    FIT --> OUT["fitted_distributions<br/>distribution_type, params, service_level_quantiles"]
    NF --> OUT

Why quantile matching over MLE?

MLE requires individual demand observations. The forecasting pipeline produces quantile forecasts, not individual samples. Quantile matching works directly from the available data without reconstruction.

Reduction to a single-period value

Each forecast row contains quantile arrays of length equal to the horizon (24 weekly steps). The fitter reduces each quantile array to a scalar by averaging across horizon steps:

\[ \hat{q}_p = \frac{1}{H}\sum_{h=1}^{H} Q_p(h) \]

This represents the average-period demand quantile. The MEIO optimizer then scales this distribution to the lead-time window when computing fill-rate integrals.

Fallback for sparse quantile sets

When only a single quantile is available (e.g., Croston / ADIDA only produce a median), the fitter creates a Normal distribution with the Poisson variance floor rather than skipping the series:

\[ \hat{\mu} = Q_{0.5}, \quad \hat{\sigma} = \max\!\left(\sqrt{\hat{\mu}},\; 10^{-3}\right) \]

This guarantees that every forecastable series has a distribution for MEIO.

Scoring and selection

Figure: Goodness-of-fit ranking — each family is scored by mean absolute quantile error; the lowest score wins.

flowchart LR
    F1["Normal<br/>loc, scale"] --> S1["score = mean |ppf(p) - q_p|"]
    F2["Gamma<br/>shape, scale"] --> S2["score"]
    F3["Negative Binomial<br/>n, p"] --> S3["score"]
    F4["LogNormal<br/>s, scale"] --> S4["score"]
    F5["Weibull<br/>c, scale"] --> S5["score"]
    F6["Poisson<br/>mu"] --> S6["score"]
    S1 & S2 & S3 & S4 & S5 & S6 --> R["Rank ascending<br/>lowest score wins"]
    R --> W["fitted_distributions.distribution_type"]

After fitting all six families, each is scored by the mean absolute error between the theoretical quantile function and the observed quantile values:

\[ \text{score}(F) = \frac{1}{|\mathcal{P}|}\sum_{p \in \mathcal{P}} |F^{-1}(p) - \hat{q}_p| \]

The family with the lowest score is selected. This is equivalent to minimizing the average quantile prediction error across the grid of observed quantile levels.

No hypothesis test gate

The current implementation does not apply a KS-test p-value threshold as a minimum acceptance criterion — it always selects the best-fitting family even if the fit is poor. The ks_statistic and ks_pvalue columns are stored in fitted_distributions for diagnostic use but are not yet used as a selection gate. A future release may enforce a minimum p-value.

Database table — fitted_distributions

CREATE TABLE {schema}.fitted_distributions (
    item_id                 BIGINT NOT NULL,
    site_id                 BIGINT NOT NULL,
    method                  TEXT NOT NULL,   -- forecasting method (e.g. 'AutoETS')
    forecast_horizon        INTEGER,
    distribution_type       TEXT,            -- 'normal'|'gamma'|'negative_binomial'|
                                             -- 'lognormal'|'weibull'|'poisson'
    mean                    DOUBLE PRECISION,
    std                     DOUBLE PRECISION,
    params                  JSONB,           -- family-specific parameters (see below)
    ks_statistic            DOUBLE PRECISION,
    ks_pvalue               DOUBLE PRECISION,
    service_level_quantiles JSONB,           -- {0.90: val, 0.95: val, 0.99: val}
    scenario_id             BIGINT NOT NULL DEFAULT 1,
    pipeline_id             BIGINT DEFAULT 0
);

The params JSONB column

Distribution Keys in params
normal loc, scale
gamma shape, scale
negative_binomial n, p
lognormal s, scale
weibull c, scale, loc
poisson mu

The service_level_quantiles column pre-computes demand quantiles at the standard service levels configured in scenario.parameters (default: \(\{0.90, 0.95, 0.99\}\)). These are consumed by the API for fast fill-rate lookups without re-evaluating the distribution.

From fitted distribution to MEIO fill rate

The Rust optimizer (files/MEIO/distributions.rs) evaluates the fill-rate integral analytically for each distribution family. For the Normal family the formula is:

\[ FR(b) = 1 - \frac{\sigma_{LT} \cdot G\!\left(\frac{b - \mu_{LT}}{\sigma_{LT}}\right)}{\mu_{LT}} \]

For Gamma, Negative Binomial, Log-Normal, Weibull, and Poisson, the backorder expectation \(E[\max(D_{LT}-b,0)]\) is computed via numerical integration of the complementary CDF:

\[ E[\max(D_{LT}-b,0)] = \int_b^{\infty} [1 - F(x)]\,dx \]

The math.rs module provides CDF implementations for all six families using accurate polynomial approximations (Normal) or recurrence relations (Poisson, Negative Binomial).

Lead-time scaling

The fitter stores single-period parameters. The optimizer scales to a lead-time window \(L\) (in weeks) using:

  • Normal: \(\mu_{LT} = L\cdot\mu_{1}\), \(\sigma_{LT} = \sqrt{L}\cdot\sigma_{1}\) (i.i.d. assumption).
  • Gamma: \(\alpha_{LT} = L\cdot\alpha_1\), \(\beta_{LT} = \beta_1\) (sum of i.i.d. Gammas is Gamma with additive shape).
  • Poisson: \(\mu_{LT} = L\cdot\mu_1\) (sum of i.i.d. Poissons is Poisson).
  • Log-Normal, Weibull, Negative Binomial: moment matching — compute \(\mu_{LT} = L\mu_1\), \(\sigma^2_{LT} = L\sigma^2_1\), then re-fit parameters of the same family from these moments.

Configuration

Distribution fitting is controlled by the meio section of scenario.parameters (parameter_type = 'meio', is_default = TRUE):

Key Default Description
distributions ['normal','gamma','negative_binomial','lognormal','weibull','poisson'] Families attempted during auto-selection
fitting_method 'quantile_matching' Only supported method currently
service_levels [0.90, 0.95, 0.99] Pre-computed quantile levels stored in service_level_quantiles
distribution_threshold 25 Minimum observation count for parametric fit; series below this threshold fall back to Normal

Operational notes

Re-running distribution fitting

Fitting runs as part of the standard pipeline after forecasting. It can also be re-triggered in isolation via the Process Runner by selecting the "Distributions" step for the active pipeline.

Diagnosing poor fits

Query the fitted_distributions table and inspect ks_statistic for unusually high values (> 0.2 is a common informal threshold). Cross-check with the series' ADI (Average Demand Interval) and zero_ratio from the alert engine to understand whether intermittency is driving the poor fit.

Alert engine integration

The alert engine reads distribution_type from fitted_distributions when evaluating rules such as "Uses Normal distribution for skewed demand". The LOWER() comparison in the LATERAL join ensures case-insensitive matching ('lognormal' not 'LogNormal').

Cross-references: Safety Stock · MEIO Theory · K-Curve