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:
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

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:
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

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)[source]

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

Parameters:

size (int | None)

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:
Return type:

None

seq_update(x, weights, estimate)[source]

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

Parameters:
Return type:

None

seq_initialize(x, weights, rng)[source]

Initialize statistics from encoded observations.

Parameters:
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

key_merge(stats_dict)[source]

Merge keyed statistics into stats_dict when keys are configured.

Parameters:

stats_dict (dict[str, Any])

Return type:

None

key_replace(stats_dict)[source]

Replace this accumulator from keyed statistics when available.

Parameters:

stats_dict (dict[str, Any])

Return type:

None

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(name=None, keys=None)[source]

Bases: ParameterEstimator

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

Parameters:

  • 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:
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