Skip to content

Ordered and Zero-Inflated Families

Ordered-categorical families, conditional and inline, and the zero-inflated, hurdle, and mixture distributions PyTorch does not ship.

ordered

Cutpoint-parameterised ordered-categorical morphisms.

Ordered-categorical distributions assign probability mass to a discrete set of ordered categories {0, 1, ..., K-1} via a latent real-valued eta and K - 1 cutpoints c_0 < c_1 < ... < c_{K-2}. The cumulative probability of landing at or below category k is CDF(c_k - eta), where CDF is either the logistic or the standard normal cumulative distribution function:

P(Y <= k | eta, c) = CDF(c_k - eta)         for 0 <= k < K - 1
P(Y = k | eta, c) = CDF(c_k - eta) - CDF(c_{k-1} - eta)
P(Y = 0 | eta, c) = CDF(c_0 - eta)
P(Y = K - 1 | eta, c) = 1 - CDF(c_{K-2} - eta)

The cutpoints argument is supplied at call time as a 1-D tensor of length K - 1, allowing the same morphism to be reused with different cutpoint vectors across observations or models.

Sampling is non-reparameterisable: the support is discrete.

ConditionalOrderedLogistic 

ConditionalOrderedLogistic(domain: AnySpace, codomain: AnySpace, num_categories: int, hidden_dim: int | Sequence[int] | None = None, param_source: ParamSource | None = None, param_source_option: str | None = None)

Bases: _ConditionalOrdered

Conditional ordered-logistic over num_categories ordered classes.

The latent eta(x) is computed by the parameter source. Given a vector of num_categories - 1 ordered cutpoints supplied at call time, the probability of category k is

P(Y = k) = sigmoid(c_k - eta) - sigmoid(c_{k-1} - eta)

with the conventions sigmoid(c_{-1} - eta) = 0 and sigmoid(c_{K-1} - eta) = 1.

PARAMETER DESCRIPTION
domain

Source space.

TYPE: SetObject or ContinuousSpace

codomain

Target space; semantically the set of ordered categories.

TYPE: SetObject or ContinuousSpace

num_categories

Number of ordered categories K >= 2.

TYPE: int

hidden_dim

Hidden layer width for the neural parameter source.

TYPE: int DEFAULT: None

Source code in src/quivers/continuous/ordered.py 
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
def __init__(
    self,
    domain: AnySpace,
    codomain: AnySpace,
    num_categories: int,
    hidden_dim: int | Sequence[int] | None = None,
    param_source: ParamSource | None = None,
    param_source_option: str | None = None,
) -> None:
    if num_categories < 2:
        raise ValueError(
            f"ordered family requires num_categories >= 2; got {num_categories}"
        )
    super().__init__(domain, codomain)
    self._k = int(num_categories)
    self.param_source = _make_source(
        domain,
        1,
        hidden_dim,
        param_source=param_source,
        param_source_option=param_source_option,
    )

ConditionalOrderedProbit 

ConditionalOrderedProbit(domain: AnySpace, codomain: AnySpace, num_categories: int, hidden_dim: int | Sequence[int] | None = None, param_source: ParamSource | None = None, param_source_option: str | None = None)

Bases: _ConditionalOrdered

Conditional ordered-probit over num_categories ordered classes.

Identical to :class:ConditionalOrderedLogistic but uses the standard normal cumulative distribution function in place of the logistic CDF.

PARAMETER DESCRIPTION
domain

Source space.

TYPE: SetObject or ContinuousSpace

codomain

Target space; semantically the set of ordered categories.

TYPE: SetObject or ContinuousSpace

num_categories

Number of ordered categories K >= 2.

TYPE: int

hidden_dim

Hidden layer width for the neural parameter source.

TYPE: int DEFAULT: None

Source code in src/quivers/continuous/ordered.py 
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
def __init__(
    self,
    domain: AnySpace,
    codomain: AnySpace,
    num_categories: int,
    hidden_dim: int | Sequence[int] | None = None,
    param_source: ParamSource | None = None,
    param_source_option: str | None = None,
) -> None:
    if num_categories < 2:
        raise ValueError(
            f"ordered family requires num_categories >= 2; got {num_categories}"
        )
    super().__init__(domain, codomain)
    self._k = int(num_categories)
    self.param_source = _make_source(
        domain,
        1,
        hidden_dim,
        param_source=param_source,
        param_source_option=param_source_option,
    )

_ordered

Ordered-categorical distribution families that PyTorch does not ship natively (OrderedLogistic, OrderedProbit).

The DSL surfaces these as inline call sites:

observe y : Resp <- OrderedLogistic(eta, cutpoints)

with eta a real predictor (one row per observation) and cutpoints either a globally shared sorted vector of length K - 1 or a per-row sorted matrix of shape (batch, K - 1). Per-row cutpoints land in the ordinal-mixed-model setting where every participant has distinct thresholds; the distribution's broadcast rules handle both shapes uniformly.

OrderedLogistic 

OrderedLogistic(predictor: Tensor, cutpoints: Tensor, validate_args: bool | None = None)

Bases: Distribution

Ordered-logit distribution over K = cutpoints.shape[-1] + 1 ordered categories indexed 0 .. K - 1.

With a real predictor eta and a sorted cutpoint vector c = (c_0 < c_1 < ... < c_{K-2}):

  • P(Y = 0) = sigmoid(c_0 - eta)
  • P(Y = k) = sigmoid(c_k - eta) - sigmoid(c_{k-1} - eta) for 0 < k < K - 1
  • P(Y = K - 1) = 1 - sigmoid(c_{K-2} - eta)

Cutpoint broadcasting:

  • Shared cutpoints c of shape (K - 1,) apply uniformly to every row.
  • Per-row cutpoints of shape (batch, K - 1) apply distinct thresholds per observation (the ordinal-mixed-model case).
  • Arbitrary leading batch shape is supported; the last axis is always the cutpoint axis.

Log probabilities are computed without forming the difference of sigmoids: with a < b the two cutpoint offsets bounding a category, sigmoid(b) - sigmoid(a) equals sigmoid(b) * sigmoid(-a) * (1 - exp(a - b)), so its log is a sum of two log-sigmoids and a log(1 - exp(.)) term, each accurate in both tails. A category far from the predictor thus keeps its exact log probability rather than one floored at the smallest normal float.

Reference: McCullagh 1980.

Source code in src/quivers/continuous/_ordered.py 
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
def __init__(
    self,
    predictor: Tensor,
    cutpoints: Tensor,
    validate_args: bool | None = None,
) -> None:
    if cutpoints.dim() == 0:
        raise ValueError(
            "OrderedLogistic: `cutpoints` must have at least one "
            f"dimension carrying the K-1 thresholds, got "
            f"shape={tuple(cutpoints.shape)}"
        )
    if cutpoints.shape[-1] < 1:
        raise ValueError(
            "OrderedLogistic: `cutpoints` last dimension must "
            f"have size >= 1 (K >= 2 categories), got "
            f"shape={tuple(cutpoints.shape)}"
        )
    self.predictor = predictor
    self.cutpoints = cutpoints
    self._num_categories = int(cutpoints.shape[-1]) + 1
    batch_shape = torch.broadcast_shapes(predictor.shape, cutpoints.shape[:-1])
    super().__init__(
        batch_shape=batch_shape,
        event_shape=torch.Size(()),
        validate_args=validate_args,
    )

_zip_hurdle

Sugar wrappers for zero-inflated, hurdle, and mixture families.

Each class here desugars at construction time to a composition in the measure algebra:

  • ZeroInflatedPoisson(pi, rate) := Mixture([pi, 1-pi], [PointMass(0), Poisson(rate)])
  • HurdlePoisson(pi, rate) := Mixture([pi, 1-pi], [PointMass(0), Restrict(Poisson(rate), low=1)])
  • MixtureNormal(w, loc, scale) := Mixture(w, [Normal(loc_k, scale_k) for k])

Internally each is a Mixture of PointMass and Restrict / base components, which the rewrite pass recognises and lifts to the canonical Bernoulli-branch + conditional-density form Stan / PyMC / NumPyro expect. A single Mixture-with-PointMass emitter covers every family in place of per-family per-base Cartesian-product renderer logic.

ZeroInflatedPoisson 

ZeroInflatedPoisson(zero_prob: Tensor, rate: Tensor, validate_args: bool | None = None)

Bases: Distribution

Mixture([pi, 1-pi], [PointMass(0), Poisson(rate)]).

Total zero mass pi + (1-pi) * exp(-rate) per Lambert 1992; the count component is the unrestricted Poisson, so excess zeros come from both the structural-zero spike and the natural Poisson zeros. Contrast with HurdlePoisson in which the count component is zero-truncated.

Reference: Lambert 1992.

Source code in src/quivers/continuous/_zip_hurdle.py 
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
def __init__(
    self,
    zero_prob: Tensor,
    rate: Tensor,
    validate_args: bool | None = None,
) -> None:
    self.zero_prob = zero_prob
    self.rate = rate
    base = torch.distributions.Poisson(rate)
    self._mixture = _two_component_inflation(zero_prob, base)
    super().__init__(
        batch_shape=self._mixture.batch_shape,
        event_shape=self._mixture.event_shape,
        validate_args=validate_args,
    )

HurdlePoisson 

HurdlePoisson(zero_prob: Tensor, rate: Tensor, validate_args: bool | None = None)

Bases: Distribution

Mixture([pi, 1-pi], [PointMass(0), Restrict(Poisson(rate), low=1)]).

Two-stage hurdle model: a Bernoulli(1 - pi) decides whether the response is zero or positive; conditional on positivity, the response is a zero-truncated Poisson. Because the two components have measure-disjoint supports, the joint log-density factors into a Bernoulli plus a zero-truncated Poisson term, and the rewrite pass emits that factored form to every backend.

Reference: Mullahy 1986.

Source code in src/quivers/continuous/_zip_hurdle.py 
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
def __init__(
    self,
    zero_prob: Tensor,
    rate: Tensor,
    validate_args: bool | None = None,
) -> None:
    self.zero_prob = zero_prob
    self.rate = rate
    base = torch.distributions.Poisson(rate)
    # Wrap the truncated component in `Normalize` so the hurdle's
    # second component is a probability measure on `{1, 2, ...}`
    # before mixing; the resulting Mixture has total mass one and
    # the user-facing `log_prob` returns the proper hurdle
    # density rather than the unnormalised sub-measure form.
    truncated_base = Normalize(Restrict(base, low=1.0))
    self._mixture = _two_component_inflation(zero_prob, truncated_base)
    super().__init__(
        batch_shape=self._mixture.batch_shape,
        event_shape=self._mixture.event_shape,
        validate_args=validate_args,
    )

MixtureNormal 

MixtureNormal(weights: Tensor, loc: Tensor, scale: Tensor, validate_args: bool | None = None)

Bases: Distribution

Mixture(weights, [Normal(loc_k, scale_k) for k]).

Same-family Gaussian mixture with K = weights.shape[-1] components. Mixing weights live in the K-simplex along the last axis; component means and scales also vary along the last axis and broadcast over the leading batch.

For mixing arbitrary families together, construct the Mixture directly with a list of heterogeneous components.

Reference: McLachlan & Peel 2000.

Source code in src/quivers/continuous/_zip_hurdle.py 
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
def __init__(
    self,
    weights: Tensor,
    loc: Tensor,
    scale: Tensor,
    validate_args: bool | None = None,
) -> None:
    if weights.shape[-1] != loc.shape[-1]:
        raise ValueError(
            "MixtureNormal: weights and loc must agree on the "
            f"component axis (last dim); got weights={tuple(weights.shape)}, "
            f"loc={tuple(loc.shape)}"
        )
    if weights.shape[-1] != scale.shape[-1]:
        raise ValueError(
            "MixtureNormal: weights and scale must agree on the "
            f"component axis (last dim); got weights={tuple(weights.shape)}, "
            f"scale={tuple(scale.shape)}"
        )
    self.weights = weights
    self.loc = loc
    self.scale = scale
    self._num_components = int(weights.shape[-1])
    components = [
        torch.distributions.Normal(loc[..., k], scale[..., k])
        for k in range(self._num_components)
    ]
    self._mixture = Mixture(weights, components)
    super().__init__(
        batch_shape=self._mixture.batch_shape,
        event_shape=self._mixture.event_shape,
        validate_args=validate_args,
    )

ZeroOneInflatedBeta 

ZeroOneInflatedBeta(mu: Tensor, phi: Tensor, zoi: Tensor, coi: Tensor, validate_args: bool | None = None)

Bases: Distribution

Zero-one inflated beta (Ospina & Ferrari 2010 / brms ZOIB).

Mixture of point masses at 0 and 1 with a continuous Beta on (0, 1):

  • P(y = 0) = zoi * (1 - coi)
  • P(y = 1) = zoi * coi
  • P(y ∈ (0, 1)) = (1 - zoi) · Beta(μ·φ, (1-μ)·φ)

Parameters follow the mean-precision Beta reparameterisation and the brms zoi / coi inflation probabilities (all on the unit interval). Formula frontend emits ZeroOneInflatedBeta(mu, phi, zoi, coi).

Source code in src/quivers/continuous/_zip_hurdle.py 
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
def __init__(
    self,
    mu: Tensor,
    phi: Tensor,
    zoi: Tensor,
    coi: Tensor,
    validate_args: bool | None = None,
) -> None:
    self.mu = mu
    self.phi = phi
    self.zoi = zoi
    self.coi = coi
    batch_shape = torch.broadcast_shapes(mu.shape, phi.shape, zoi.shape, coi.shape)
    super().__init__(
        batch_shape=batch_shape,
        event_shape=torch.Size(),
        validate_args=validate_args,
    )