How to fit a multimodal negative binomial with autodiff

This guide shows how to fit demand histories that contain low regular values and rare large quantities, such as 0, 1, and 1000. The fitted model mixes a zero-inflated negative binomial for the low mode with a negative binomial for the high mode.

Use this model when a single negative binomial assigns unwanted probability to the quantities between two distinct modes.

Step 1: Define a stable mixture likelihood

Combine the component log-likelihoods with logAddExp. This formulation avoids exponentiating both complete likelihoods before adding them, which can lose the smaller component through numerical underflow.

def autodiff pure logAddExp(a: number, b: number) with
  m = if a > b then a else b
  return m + log(exp(a - m) + exp(b - m))

def autodiff pure bimodalNegativeBinomialLL(
    lowMean: number,
    lowDispersion: number,
    zeroInflation: number,
    lowWeight: number,
    highMean: number,
    highDispersion: number,
    k: number) with
  low = log(lowWeight) + loglikelihood.negativeBinomial(
    lowMean, lowDispersion, zeroInflation, k)
  high = log(1 - lowWeight) + loglikelihood.negativeBinomial(
    highMean, highDispersion, k)
  return logAddExp(low, high)

table Items[item] = with
  [| as item |]
  [| "A" |]

table Weeks[week] = with
  [| as week, as BaseDemand |]
  [| 1, 0 |]
  [| 2, 0 |]
  [| 3, 0 |]
  [| 4, 0 |]
  [| 5, 0 |]
  [| 6, 0 |]
  [| 7, 0 |]
  [| 8, 0 |]
  [| 9, 0 |]
  [| 10, 0 |]
  [| 11, 1 |]
  [| 12, 1 |]
  [| 13, 1 |]
  [| 14, 1 |]
  [| 15, 1 |]
  [| 16, 1000 |]

// Synthetic one-item history for this stand-alone example.
// In production, ItemsWeek contains the item-specific demand histories.
table ItemsWeek = cross(Items, Weeks)
ItemsWeek.Demand = Weeks.BaseDemand

outlierScale = 1000

autodiff Items epochs:500 learningRate:0.03 with
  params Items.lowMean in [0.001 .. 10] auto(0.5, 0)
  params Items.lowDispersion in [1.001 .. 20] auto(2, 0)
  params Items.zeroInflation in [0 .. 0.999] auto(0.5, 0)
  params Items.lowWeight in [0.001 .. 0.999] auto(0.9, 0)
  params Items.highMeanScale in [0.1 .. 2] auto(1, 0)
  params Items.highDispersion in [1.001 .. 20] auto(2, 0)

  highMean = outlierScale * Items.highMeanScale
  itemLogLikelihood = sum(bimodalNegativeBinomialLL(
    Items.lowMean,
    Items.lowDispersion,
    Items.zeroInflation,
    Items.lowWeight,
    highMean,
    Items.highDispersion,
    ItemsWeek.Demand))
  return -itemLogLikelihood

Items.highMean = outlierScale * Items.highMeanScale
Items.Forecast = mixture(
  negativeBinomial(Items.lowMean, Items.lowDispersion, Items.zeroInflation),
  Items.lowWeight,
  negativeBinomial(Items.highMean, Items.highDispersion))

Items.pAt50 = valueAt(cdf(Items.Forecast), 50) -
    valueAt(cdf(Items.Forecast), 49)

show table "Fitted model by item" a4b4 with
  Items.item
  Items.lowMean as "Low component mean"
  Items.lowDispersion as "Low component dispersion"
  Items.zeroInflation as "Zero inflation"
  Items.lowWeight as "Low component weight"
  Items.highMean as "High component mean" // ~1000
  Items.highDispersion as "High component dispersion"
  Items.pAt50 as "P(Demand = 50)" // 0 with this ranvar representation

Replace the synthetic Items, Weeks, and ItemsWeek tables with the item master and item-week history from the target dataset. Keep Items as the observation table of the autodiff block, store the parameters on Items, and sum the weekly log-likelihoods from ItemsWeek. Set outlierScale to the order of magnitude of the high mode instead of fitting a thousand-unit mean directly.

Step 2: Check the fitted modes

Run the script and verify:

• Items.highMean is close to the large observations rather than the low values.
• Items.lowWeight is strictly between 0 and 1.
• Items.pAt50 is negligible compared with the probability around either mode.

The likelihood is not convex. If the two modes collapse into one, initialize Items.lowMean, Items.highMeanScale, and Items.lowWeight from a preliminary split of the history, or run the fit from several initial values.

Keep Items.lowDispersion and Items.highDispersion at or above 1.001. Below this threshold, loglikelihood.negativeBinomial uses its Poisson fallback and the dispersion does not receive a gradient.

The two negative binomial components have non-negative integer support. Their mixture can make the probability between the modes negligible, but it does not create an exact structural gap in the underlying mathematical distribution.

Step 3: Reuse the fitted forecast

Use Items.Forecast outside the autodiff block like any other ranvar. For example, compute item-level quantiles or combine it with independent future periods. The ranvar constructors belong after the fit; use the differentiable loglikelihood.* functions inside the autodiff block.

See the negativeBinomial and mixture references for the materialized ranvar operations.