mixle.stats.directional.von_mises_fisher module

Create, estimate, and sample from a von Mises-Fisher distribution.

Defines the VonMisesFisherDistribution, VonMisesFisherSampler, VonMisesFisherAccumulatorFactory, VonMisesFisherAccumulator, VonMisesFisherEstimator, and the VonMisesFisherDataEncoder classes for use with mixle.

Data type: Union[Sequence[float], np.ndarray] (a unit-norm vector on the (p-1)-sphere in R^p).

The von Mises-Fisher (vmf) distribution on the (p-1) sphere in R^{p}. Assume x_mat = (X_1,..,X_p) follows a vmf distribution with mean direction vector mu = (mu_1, mu_2, …, mu_p) s.t. ||mu||=1 and concentration parameter kappa > 0. The vmf log-density if given by

log(f(x; mu, kappa)) = log(c_p(kappa)) + kappa * dot(mu, x),

where dot is a dot product and

log(c_p(kappa)) = (p/2-1)log(kappa) - (p/2)*log(2*pi) + log(B_{p/2-1}(kappa)), where

log(B_{p/2-1}(kappa)) = denotes the modified Bessel function of the first kind at order p/2-1.

Numerical notes:

Evaluating log I_v(kappa) directly with scipy.special.iv overflows for large kappa, and the exponentially scaled scipy.special.ive underflows when the order v = p/2 - 1 is large relative to kappa (high dimension with modest concentration). The helper lniv() therefore uses log(ive) + kappa where ive has support and falls back to the uniform large-order asymptotic expansion (Abramowitz & Stegun 9.7.7) implemented in lniv_uniform() when ive underflows. Both the normalizing constant and the Bessel-ratio Newton iteration in VonMisesFisherEstimator.estimate() rely on lniv().

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

lniv_uniform(v, ln_z)

log I_v(z) by the uniform large-order asymptotic (A&S 9.7.7):

I_v(v t) ~ exp(v eta) / (sqrt(2 pi v) (1 + t^2)^{1/4}), eta = sqrt(1 + t^2) + log(t / (1 + sqrt(1 + t^2))).

Valid uniformly in t = z/v for large v, including t -> 0 where it reduces to the small-argument form (z/2)^v / Gamma(v+1) via Stirling.

Parameters:

• v (float) – Order of the modified Bessel function. Must be positive.

• ln_z (float) – Log of the (positive) argument z.

Returns:

Approximate value of log I_v(z) as a float.

lniv(v, ln_z)

Numerically stable log I_v(e^{ln_z}).

Uses the exponentially scaled Bessel function where it has support and the uniform large-order expansion where ive underflows (large v relative to z; ive cannot underflow for v = 0, so that branch always has v > 0).

Parameters:

• v (float) – Order of the modified Bessel function. Must be non-negative.

• ln_z (float) – Log of the argument z. May be -inf (z = 0).

Returns:

log I_v(z) as a float (-inf when z = 0 and v > 0).

class VonMisesFisherDistribution(mu, kappa, name=None, keys=None)

Bases: SequenceEncodableProbabilityDistribution

Von Mises-Fisher distribution on the (p-1)-sphere with mean direction mu and concentration kappa.

Data type: Union[Sequence[float], np.ndarray] (a unit-norm vector in R^p).

Parameters:

• mu (Sequence[float] | ndarray)

• kappa (float)

• name (str | None)

• keys (str | None)

classmethod compute_capabilities()

classmethod compute_declaration()

static backend_legacy_sufficient_statistics(x, params, engine)

Return row-wise legacy sufficient statistics for resident reductions.

Parameters:

• x (Any)

• params (dict[str, Any])

• engine (Any)

Return type:

tuple[Any, …]

static backend_log_density_from_params(x, mu, kappa, log_const, engine)

Engine-neutral von Mises-Fisher log-density from fitted parameters.

Parameters:

• x (Any)

• mu (Any)

• kappa (Any)

• log_const (Any)

• engine (Any)

Return type:

Any

density(x)

Density of von Mises-Fisher distribution at observation x.

See log_density() for details.

Parameters:

x (Union[Sequence[float], np.ndarray]) – Unit-norm vector in R^p.

Returns:

Density at observation x.

Return type:

float

log_density(x)

Log-density of von Mises-Fisher distribution at observation x.

The log-density is given by

log(f(x; mu, kappa)) = log(c_p(kappa)) + kappa * dot(mu, x),

for x on the (p-1)-sphere. When kappa = 0 this reduces to the uniform density on the sphere.

Parameters:

x (Union[Sequence[float], np.ndarray]) – Unit-norm vector in R^p.

Returns:

Log-density at observation x.

Return type:

float

density_cumulative(x)

Exact probability-ordered cumulative G(x) = P(p(Y) >= p(x)) (the HDR mass at x).

A coordinate-wise CDF is undefined on the sphere (no total order), but since the density is monotone in the cosine t = mu . x (p(y) >= p(x) iff mu.y >= mu.x for kappa >= 0), the highest-density-region mass is the upper tail of the cosine marginal, whose density is f(s) proportional to exp(kappa s) (1 - s^2)^((p-3)/2) on [-1, 1]. G is that tail integral (computed by quadrature; the exp(kappa(s-1)) shift keeps it stable for large kappa and cancels in the ratio). Returned to density_rank as method exact-analytic.

Parameters:

x (Sequence[float] | ndarray)

Return type:

float

density_quantile(q)

Inverse of density_cumulative(): a representative unit vector at cumulative-density q.

q is the highest-density-region mass; since the density is monotone in the cosine t = mu . y, the boundary is the cosine t_q with tail mass q, found by bisection on the cosine marginal. The returned representative is a unit vector at that cosine from mu (t_q * mu + sqrt(1 - t_q^2) * perp for a fixed perp orthogonal to mu). Sweeping q enumerates the sphere in descending density (concentric caps about mu).

Parameters:

q (float)

Return type:

ndarray

seq_log_density(x)

Vectorized evaluation of log-density at sequence encoded input x.

Parameters:

x (np.ndarray) – 2-d numpy array of N unit-norm vectors with p columns.

Returns:

Numpy array of log-density (float) of length N.

Return type:

ndarray

backend_seq_log_density(x, engine)

Engine-neutral vectorized log-density for encoded unit-vector observations.

Parameters:

• x (ndarray)

• engine (Any)

Return type:

Any

classmethod backend_stacked_params(dists, engine)

Return stacked parameters for equal-dimensional von Mises-Fisher mixtures.

Parameters:

• dists (Sequence[VonMisesFisherDistribution])

• engine (Any)

Return type:

dict[str, Any]

classmethod backend_stacked_log_density(x, params, engine)

Return an (n, k) matrix of von Mises-Fisher component log densities.

Parameters:

• x (ndarray)

• params (dict[str, Any])

• engine (Any)

Return type:

Any

classmethod backend_stacked_sufficient_statistics(x, weights, params, engine)

Return component-stacked legacy (count, weighted_vector_sum) statistics.

Parameters:

• x (ndarray)

• weights (Any)

• params (dict[str, Any])

• engine (Any)

Return type:

tuple[Any, Any]

sampler(seed=None)

Create a VonMisesFisherSampler object from parameters of VonMisesFisherDistribution instance.

Parameters:

seed (Optional[int]) – Used to set seed in random sampler.

Returns:

VonMisesFisherSampler object.

Return type:

VonMisesFisherSampler

estimator(pseudo_count=None)

Create a VonMisesFisherEstimator object.

Parameters:

pseudo_count (Optional[float]) – Kept for interface consistency (has no effect on estimation).

Returns:

VonMisesFisherEstimator object.

Return type:

VonMisesFisherEstimator

dist_to_encoder()

Returns a VonMisesFisherDataEncoder object for encoding sequences of data.

Return type:

VonMisesFisherDataEncoder

class VonMisesFisherSampler(dist, seed=None)

Bases: DistributionSampler

Sampler for the VonMisesFisherDistribution using Wood’s rejection sampling scheme.

Parameters:

• dist (VonMisesFisherDistribution)

• seed (int | None)

sample(size=None)

Draw iid unit-norm vectors from the von Mises-Fisher distribution.

Parameters:

size (Optional[int]) – Number of samples to draw. If None, a single vector is returned.

Returns:

Numpy array of shape (dim,) if size is None, else of shape (size, dim).

Return type:

ndarray

class VonMisesFisherAccumulator(dim=None, name=None, keys=None)

Bases: SequenceEncodableStatisticAccumulator

Accumulator for the VonMisesFisherDistribution. Tracks the weighted vector sum and total weight.

Parameters:

• dim (int | None)

• name (str | None)

• keys (str | None)

update(x, weight, estimate)

Update sufficient statistics with a weighted observation.

Parameters:

• x (Union[Sequence[float], np.ndarray]) – Unit-norm vector in R^p.

• weight (float) – Weight for observation.

• estimate (Optional[VonMisesFisherDistribution]) – Previous estimate (unused).

Return type:

None

initialize(x, weight, rng)

Initialize sufficient statistics with a weighted observation.

Parameters:

• x (Union[Sequence[float], np.ndarray]) – Unit-norm vector in R^p.

• weight (float) – Weight for observation.

• rng (RandomState) – Random number generator (unused).

Return type:

None

seq_update(x, weights, estimate)

Vectorized update of sufficient statistics from sequence encoded data.

Non-finite or negative weights are dropped from the vector sum.

Parameters:

• x (np.ndarray) – 2-d numpy array of N unit-norm vectors with p columns.

• weights (np.ndarray) – Weights for each of the N observations.

• estimate (Optional[VonMisesFisherDistribution]) – Previous estimate (unused).

Return type:

None

seq_initialize(x, weights, rng)

Vectorized initialization of sufficient statistics from sequence encoded data.

Parameters:

• x (np.ndarray) – 2-d numpy array of N unit-norm vectors with p columns.

• weights (np.ndarray) – Weights for each of the N observations.

• rng (RandomState) – Random number generator (unused).

Return type:

None

combine(suff_stat)

Combine sufficient statistics from another accumulator into this one.

Parameters:

suff_stat (Tuple[float, np.ndarray]) – Tuple of count and weighted vector sum.

Returns:

Self, with aggregated sufficient statistics.

Return type:

VonMisesFisherAccumulator

value()

Returns sufficient statistics as a Tuple of count and weighted vector sum.

Return type:

tuple[float, ndarray]

from_value(x)

Set sufficient statistics of accumulator from value x.

Parameters:

x (Tuple[float, np.ndarray]) – Tuple of count and weighted vector sum.

Return type:

VonMisesFisherAccumulator

key_merge(stats_dict)

Merge sufficient statistics of object instance with suff stats containing matching keys.

Parameters:

stats_dict (Dict[str, Any]) – Dict mapping keys to accumulators with shared sufficient statistics.

Returns:

None.

Return type:

None

key_replace(stats_dict)

Set sufficient statistics of object instance to suff stats with matching keys.

Parameters:

stats_dict (Dict[str, Any]) – Dict mapping keys to accumulators with shared sufficient statistics.

Returns:

None.

Return type:

None

acc_to_encoder()

Returns a VonMisesFisherDataEncoder object for encoding sequences of data.

Return type:

VonMisesFisherDataEncoder

class VonMisesFisherAccumulatorFactory(dim=None, name=None, keys=None)

Bases: StatisticAccumulatorFactory

Factory for creating VonMisesFisherAccumulator objects.

Parameters:

• dim (int | None)

• name (str | None)

• keys (str | None)

make()

Returns a new VonMisesFisherAccumulator object.

Return type:

SequenceEncodableStatisticAccumulator

class VonMisesFisherEstimator(dim=None, pseudo_count=None, name=None, keys=None)

Bases: ParameterEstimator

Estimator for the VonMisesFisherDistribution using the Banerjee et al. approximation for kappa.

Parameters:

• dim (int | None)

• pseudo_count (float | None)

• name (str | None)

• keys (str | None)

accumulator_factory()

Returns a VonMisesFisherAccumulatorFactory for creating VonMisesFisherAccumulator objects.

estimate(nobs, suff_stat)

Estimate a VonMisesFisherDistribution from sufficient statistics.

The mean direction mu is the normalized weighted vector sum. The concentration kappa solves A_p(kappa) = rhat (rhat = ||ssum|| / count), initialized with the closed-form Banerjee et al. approximation and refined with up to three Newton steps. The Bessel-function ratio A_p(kappa) is evaluated through lniv() so large orders p/2 fall back to the uniform large-order asymptotic instead of underflowing. Two guards keep the solution finite: rhat is clamped below 1 (rhat -> 1 sends kappa -> inf), and Newton refinement is skipped within 1e-9 of 1 where A_p’(kappa) -> 0 makes the iteration ill-conditioned while the initializer is already accurate.

Parameters:

• nobs (Optional[float]) – Number of observations (unused).

• suff_stat (Tuple[float, np.ndarray]) – Tuple of count and weighted vector sum.

Returns:

VonMisesFisherDistribution object (uniform on the sphere, kappa = 0, if no data observed).

Return type:

VonMisesFisherDistribution

class VonMisesFisherDataEncoder

Bases: DataSequenceEncoder

Data encoder for sequences of unit-norm vector observations.

seq_encode(x)

Encode a sequence of N unit-norm vectors for vectorized functions.

Parameters:

x (Union[Sequence[float], np.ndarray]) – Sequence of N unit-norm vectors in R^p.

Returns:

2-d numpy array with N rows and p columns.

Return type:

ndarray