mixle.stats.directional.watson module

Watson distribution – a rotationally symmetric distribution for axial data on the sphere.

Axial data are unit vectors identified with their antipodes (x and -x are the same), e.g. fibre or crystal orientations, where the von Mises-Fisher (which distinguishes x from -x) does not apply. The Watson distribution on S^{p-1} concentrates around an axis mu with shape kappa:

f(x; mu, kappa) = M(1/2, p/2, kappa)^{-1} / omega_p * exp(kappa (mu^T x)^2),

where M is Kummer’s confluent hypergeometric function and omega_p = 2 pi^{p/2} / Gamma(p/2) is the sphere’s surface area. kappa > 0 is bipolar (mass near the axis +/-mu), kappa < 0 girdle (mass on the equator orthogonal to mu); both are antipodally symmetric. It is fit by maximum likelihood: mu is the leading (kappa>0) or trailing (kappa<0) eigenvector of the scatter matrix, and kappa solves E[(mu^T x)^2] = mu^T S mu (a monotone 1-D equation in the Kummer ratio).

Reference: Mardia & Jupp, Directional Statistics (Wiley, 2000).

class WatsonDistribution(mu, kappa, name=None, keys=None)[source]

Bases: SequenceEncodableProbabilityDistribution

Watson distribution on the unit sphere S^{p-1} with axis mu and concentration kappa.

Parameters:
density(x)[source]

Return the density at a single unit vector x.

Parameters:

x (ndarray)

Return type:

float

log_density(x)[source]

Return the log-density at a single unit vector x.

Parameters:

x (ndarray)

Return type:

float

seq_log_density(x)[source]

Vectorized log-density for a stack of unit vectors, shape (N, p).

Parameters:

x (ndarray)

Return type:

ndarray

classmethod compute_capabilities()[source]
backend_seq_log_density(x, engine)[source]

Engine-neutral vectorized log-density for (N, p) unit vectors.

Parameters:
Return type:

Any

sampler(seed=None)[source]

Return a sampler for drawing unit vectors from this distribution.

Parameters:

seed (int | None)

Return type:

WatsonSampler

estimator(pseudo_count=None)[source]

Return a maximum-likelihood estimator (scatter eigenvector + Kummer-ratio kappa solve).

Parameters:

pseudo_count (float | None)

Return type:

WatsonEstimator

dist_to_encoder()[source]

Return the data encoder used by this distribution for vectorized methods.

Return type:

WatsonDataEncoder

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

Bases: DistributionSampler

Draw axes by sampling s = |mu^T x| via a numerical inverse-CDF, then a random orthogonal tangent.

Parameters:

  • dist (WatsonDistribution)

  • seed (int | None)

sample(size=None)[source]

Draw observations.

Combinator samplers (mixture/sequence/…) accept batched. With batched=True (the default) each child stream is drawn in one vectorized call instead of a per-draw Python loop – far faster. Because every child sampler owns an independent RandomState, batching consumes each stream in the same order as the loop, so the draws are identical to the legacy path. batched=False forces that legacy per-draw loop as a guaranteed- stable reference. Leaf samplers are already vectorized and ignore the flag.

Parameters:

size (int | None)

Return type:

ndarray

class WatsonAccumulator(dim, name=None, keys=None)[source]

Bases: SequenceEncodableStatisticAccumulator

Accumulate the weighted scatter matrix S = sum_i w_i x_i x_i^T and total weight.

Parameters:

  • dim (int)

  • name (str | None)

  • keys (str | None)

update(x, weight, estimate)[source]
Parameters:

  • x (ndarray)

  • weight (float)

  • estimate (WatsonDistribution | None)

Return type:

None

initialize(x, weight, rng)[source]
Parameters:
Return type:

None

seq_update(x, weights, estimate)[source]
Parameters:
Return type:

None

seq_initialize(x, weights, rng)[source]
Parameters:
Return type:

None

combine(suff_stat)[source]
Parameters:

suff_stat (tuple[ndarray, float])

Return type:

WatsonAccumulator

value()[source]
Return type:

tuple[ndarray, float]

from_value(x)[source]
Parameters:

x (tuple[ndarray, float])

Return type:

WatsonAccumulator

key_merge(stats_dict)[source]

Pool this accumulator’s statistics into stats_dict under its merge key.

The structural default implements the common single-key pattern: store the accumulator under self.keys the first time the key is seen, else combine into the one already there. Accumulators with several named keys (e.g. an HMM’s init/trans/state keys) or a non-accumulator stats payload override this. A keys of None (the default) is a no-op.

Parameters:

stats_dict (dict[str, Any])

Return type:

None

key_replace(stats_dict)[source]

Replace this accumulator’s statistics from the pooled stats_dict entry (see key_merge).

Parameters:

stats_dict (dict[str, Any])

Return type:

None

acc_to_encoder()[source]
Return type:

WatsonDataEncoder

class WatsonAccumulatorFactory(dim, name=None, keys=None)[source]

Bases: StatisticAccumulatorFactory

Factory for WatsonAccumulator.

Parameters:

  • dim (int)

  • name (str | None)

  • keys (str | None)

make()[source]
Return type:

WatsonAccumulator

class WatsonEstimator(dim, name=None, keys=None)[source]

Bases: ParameterEstimator

Maximum-likelihood estimator: scatter eigenvector for the axis, Kummer-ratio solve for kappa.

Parameters:

  • dim (int)

  • name (str | None)

  • keys (str | None)

accumulator_factory()[source]
Return type:

WatsonAccumulatorFactory

estimate(nobs, suff_stat)[source]
Parameters:
Return type:

WatsonDistribution

class WatsonDataEncoder[source]

Bases: DataSequenceEncoder

Encode a sequence of unit vectors as an (N, p) float array.

seq_encode(x)[source]

Encode the iid observation sequence x for vectorized evaluation.

Parameters:

x (Sequence[ndarray])

Return type:

ndarray