mixle.stats.univariate.discrete.skellam module¶

Skellam distributions for differences of independent Poisson counts.

Data type (int): K = N1 - N2 with N1 ~ Poisson(mu1), N2 ~ Poisson(mu2) independent, so

K ranges over all integers (negative, zero, positive). Its log-mass is

log p(k) = -(sqrt(mu1) - sqrt(mu2))^2 + (k/2) * log(mu1/mu2) + log(I_|k|(2*sqrt(mu1*mu2))),

where I_v is the modified Bessel function of the first kind. The exponentially-scaled ive(v, z) = I_v(z) * exp(-z) is used (log I_v(z) = log(ive(v, z)) + z) so the Bessel term does not overflow for large z; combined with -(mu1+mu2) + z this collapses to the stable -(sqrt(mu1) - sqrt(mu2))^2 constant above.

The MLE has no closed form, but the method of moments is exact and closed-form here: with sample mean m and variance v, mu1 = (v + m)/2 and mu2 = (v - m)/2 (since E[K] = mu1-mu2 and Var[K] = mu1+mu2), which the estimator uses (clamped to keep both rates positive).

Reference: Skellam, ‘The frequency distribution of the difference between two Poisson variates’, JRSS A (1946).

class SkellamDistribution(mu1, mu2, name=None, keys=None)[source]

Bases: SequenceEncodableProbabilityDistribution

Skellam distribution: K = N1 - N2 for independent N1 ~ Poisson(mu1), N2 ~ Poisson(mu2).

Parameters:

mu1 (float)
mu2 (float)
name (str | None)
keys (str | None)

density(x)[source]

Probability mass at integer x (see log_density).

Parameters:: x (int)
Return type:: float

log_density(x)[source]

Stable Skellam log-mass at integer x (-inf for non-integer input).

Parameters:: x (int)
Return type:: float

seq_log_density(x)[source]

Vectorized Skellam log-mass at sequence-encoded integer counts x.

Parameters:: x (ndarray)
Return type:: ndarray

mean()[source]

Mean E[X] = mu1 - mu2.

Return type:: float

variance()[source]

Variance Var[X] = mu1 + mu2.

Return type:: float

entropy()[source]

Shannon entropy in nats, by exact two-sided summation of the standard series.

The Skellam entropy -sum_k p_k log p_k has no closed form. The sum runs over the support between this distribution’s own quantiles at 1e-16 and 1 - 1e-16 (plus a safety margin on each side), beyond which the remaining tail mass on both sides is far below double rounding. Terms where the scaled Bessel term underflows to exactly zero (mass indistinguishable from 0 at double precision) are dropped, matching the p log p -> 0 limit.

Return type:: float

classmethod compute_capabilities()[source]: Describe backend support for generated Skellam scoring kernels.

backend_seq_log_density(x, engine)[source]

Engine-neutral vectorized Skellam log-mass for encoded data (see class backend note).

The series yields the UNSCALED log I_k (not log ive = log I_k - z), so the constant here is -(mu1 + mu2) — the legacy path’s -sqrt_diff_sq plus its implicit -z.

Parameters:

x (Any)
engine (Any)

Return type:

Any

cdf(x)[source]

Cumulative distribution function P(X <= x) (via scipy skellam).

Parameters:: x (float)
Return type:: float

quantile(q)[source]

Inverse CDF F^{-1}(q) (via scipy skellam).

Parameters:: q (float)
Return type:: float

sampler(seed=None)[source]

Return a SkellamSampler for this distribution.

Parameters:: seed (int | None)
Return type:: SkellamSampler

estimator(pseudo_count=None)[source]

Return a SkellamEstimator (method of moments).

Parameters:: pseudo_count (float | None)
Return type:: SkellamEstimator

dist_to_encoder()[source]

Return the encoder for Skellam observations.

Return type:: SkellamDataEncoder

enumerator()[source]

Return an enumerator over all integers in descending probability order.

Return type:: SkellamEnumerator

class SkellamEnumerator(dist)[source]

Bases: DistributionEnumerator

Enumerate the two-sided Skellam support in descending probability order.

The Skellam pmf is log-concave (the difference of two independent Poisson counts, each log-concave on the integers), hence unimodal: enumeration hill-climbs from round(mu1 - mu2) to the mode and walks outward with the standard two-pointer merge over the two tails, emitting the larger side each step. Values whose mass underflows to zero are skipped per the enumerator contract, so the iterator ends when both tails are exhausted at float precision.

Parameters:: dist (SkellamDistribution)

class SkellamSampler(dist, seed=None)[source]

Bases: DistributionSampler

Draw iid Skellam observations as the difference of two independent Poisson draws.

Parameters:

dist (SkellamDistribution)
seed (int | None)

sample(size=None, *, batched=True)[source]

Draw size iid Skellam samples (a single int if size is None).

Parameters:

size (int | None)
batched (bool)

Return type:

int | ndarray

class SkellamAccumulator(name=None, keys=None)[source]

Bases: SequenceEncodableStatisticAccumulator

Accumulate the weighted count, sum, and sum-of-squares needed for the moment fit.

Parameters:

name (str | None)
keys (str | None)

update(x, weight, estimate)[source]

Accumulate count, sum, and squared sum for one integer observation.

Parameters:

x (int)
weight (float)
estimate (SkellamDistribution | None)

Return type:

None

initialize(x, weight, rng)[source]

Initialize statistics from one observation.

Parameters:

x (int)
weight (float)
rng (RandomState | None)

Return type:

None

seq_update(x, weights, estimate)[source]

Accumulate count, sum, and squared sum from encoded observations.

Parameters:

x (ndarray)
weights (ndarray)
estimate (SkellamDistribution | None)

Return type:

None

seq_initialize(x, weights, rng)[source]

Initialize statistics from encoded observations.

Parameters:

x (ndarray)
weights (ndarray)
rng (RandomState | None)

Return type:

None

combine(suff_stat)[source]

Merge another Skellam sufficient-statistic tuple.

Parameters:: suff_stat (tuple[float, float, float])
Return type:: SkellamAccumulator

value()[source]

Return count, sum, and squared sum.

Return type:: tuple[float, float, float]

from_value(x)[source]

Replace accumulator contents from a sufficient-statistic tuple.

Parameters:: x (tuple[float, float, float])
Return type:: SkellamAccumulator

scale(c)[source]

Scale all weight-linear sufficient statistics by c.

Parameters:: c (float)
Return type:: SkellamAccumulator

acc_to_encoder()[source]

Return the encoder used by this accumulator.

Return type:: SkellamDataEncoder

class SkellamAccumulatorFactory(name=None, keys=None)[source]

Bases: StatisticAccumulatorFactory

Factory for SkellamAccumulator.

Parameters:

name (str | None)
keys (str | None)

make()[source]

Create a fresh Skellam accumulator.

Return type:: SkellamAccumulator

class SkellamEstimator(pseudo_count=None, suff_stat=None, name=None, keys=None)[source]

Bases: ParameterEstimator

Estimate (mu1, mu2) by the (exact, closed-form) method of moments.

Parameters:

pseudo_count (float | None)
suff_stat (tuple[float, float] | None)
name (str | None)
keys (str | None)

accumulator_factory()[source]

Return an accumulator factory for Skellam moment statistics.

Return type:: SkellamAccumulatorFactory

estimate(nobs, suff_stat)[source]

Estimate a Skellam from the accumulated (count, sum, sum2) via method of moments.

Parameters:

nobs (float | None)
suff_stat (tuple[float, float, float])

Return type:

SkellamDistribution

class SkellamDataEncoder[source]

Bases: DataSequenceEncoder

Encode sequences of iid Skellam observations (integer data type, any sign).

seq_encode(x)[source]

Encode observations as a floating-point integer-valued array.

Parameters:: x (Sequence[int])
Return type:: ndarray