diff --git a/README.md b/README.md index 508a73c..4576e70 100644 --- a/README.md +++ b/README.md @@ -19,7 +19,7 @@ Refer to the [online documentation](https://spherical-stats.readthedocs.io/en/la * Orientation tensor * Parametric distributions with scipy.stats like API: * Modeling axial data: Angular central gaussian distribution (ACG) - * Modeling vector data: Elliptically symmetrical angular gausian distribution (ESAG) + * Modeling vector data: Elliptically symmetrical angular gausian distribution (ESAG), Von Mises-Fisher distribution (VMF) Example usage of the distributions: diff --git a/setup.py b/setup.py index 55efdfb..57d7a41 100644 --- a/setup.py +++ b/setup.py @@ -7,7 +7,7 @@ setup( name='spherical_stats', - version='1.1', + version='1.2', description='Spherical statistics in Python', author='Daniel Schmitz', author_email='danielschmitzsiegen@gmail.com', diff --git a/spherical_stats/_acg.py b/spherical_stats/_acg.py index 48bb1b7..d2a3856 100644 --- a/spherical_stats/_acg.py +++ b/spherical_stats/_acg.py @@ -66,12 +66,21 @@ def fit(vectors, tol): return l class ACG(object): - ''' + r''' Angular Central Gaussian distribution Args: cov_matrix (optional, ndarray (3, 3) ): Covariance matrix of the ACG distribution + The angular central gaussian distribution is an elliptically symmetrical distribution + for axial data. Its PDF is defined as + + .. math:: + + p_{ACG}(\mathbf{x}|\mathbf{\Lambda}) = \frac{1}{4\pi\sqrt{|\Lambda|}}(\mathbf{x}\Lambda^{-1}\mathbf{x})^{-\frac{3}{2}} + + with covariance matrix :math:`\Lambda` and its determinant :math:`|\Lambda|` . + Notes ------- Reference: Tyler, Statistical analysis for the angular central Gaussian distribution on the sphere diff --git a/spherical_stats/_esag.py b/spherical_stats/_esag.py index 8e26579..5e1be3c 100644 --- a/spherical_stats/_esag.py +++ b/spherical_stats/_esag.py @@ -222,22 +222,27 @@ def _fit(vectors, print_summary = False): return optimized_params class ESAG(object): - ''' + r''' Elliptically symmetrical angular Central Gaussian distribution Args: params (optional, ndarray (5, ) ): Parameters of the distribution - ``params`` are the following: :math:`(\mu_0, \mu_1, \mu_2, \gamma_1, \gamma_2)`.\n - The principal orientation vectors is given by the normalized vector :math:`\mu=(\mu_0, \mu_1, \mu_2)/||(\mu_0, \mu_1, \mu_2)||`.\n + ``params`` are the following: :math:`(\mu_0, \mu_1, \mu_2, \gamma_1, \gamma_2)`. + + The principal orientation vectors is given by the normalized vector :math:`\boldsymbol{\mu}=[\mu_0, \mu_1, \mu_2]^T/||[\mu_0, \mu_1, \mu_2]^T||` and the shape of the distribution is controlled by the parameters :math:`\gamma_1` and :math:`\gamma_2`. Notes ------- - Reference: Paine et al. An elliptically symmetric angular Gaussian distribution, - Statistics and Computing volume 28, 689–697 (2018)\n - Unlike the original Matlab implementation the distribution is fitted using the L-BFGS-B algorithm + The formula of the ESAG PDF is quite complicated, developers are referred to the reference below. + + Note that unlike the original Matlab implementation the distribution is fitted using the L-BFGS-B algorithm based on a finite difference approximation of the gradient. So far, this has proven to work succesfully. + + Reference: Paine et al. An elliptically symmetric angular Gaussian distribution, + Statistics and Computing volume 28, 689–697 (2018) + ''' def __init__(self, params = None): diff --git a/spherical_stats/_vmf.py b/spherical_stats/_vmf.py index 2c25f54..826934c 100644 --- a/spherical_stats/_vmf.py +++ b/spherical_stats/_vmf.py @@ -67,19 +67,32 @@ def obj(s): return mu, kappa class VMF: - ''' - Van Mises-Fisher distribution + r""" + Von Mises-Fisher distribution Args: - mu (optional, ndarray (3, ) ): Mean orientation + mu (optional, ndarray (3, ) ): Mean orientation kappa (optional, float): positive concentration parameter - Notes - ------- - References:\n - Mardia, Jupp. Directional Statistics, 1999. \n - Numerically stable sampling of the von Mises Fisher distribution on S2. Wenzel, 2012 - ''' + The VMF distribution is an isotropic distribution for + directional data. Its PDF is typically defined as + + .. math:: + + p_{vMF}(\mathbf{x}| \boldsymbol{\mu}, \kappa) = \frac{\kappa}{4\pi\cdot\text{sinh}(\kappa)}\exp(\kappa \boldsymbol{\mu}^T\mathbf{x}) + + Here, the numerically stable variant from (Wenzel, 2012) is used: + + .. math:: + + p_{vMF}(\mathbf{x}| \boldsymbol{\mu}, \kappa) = \frac{\kappa}{2\pi(1-\exp(-2\kappa))}\exp(\kappa( \boldsymbol{\mu}^T \mathbf{x}-1)) + + References: + + Mardia, Jupp. Directional Statistics, 1999. + + Wenzel. Numerically stable sampling of the von Mises Fisher distribution on S2. 2012 + """ def __init__(self, mu = None, kappa = None): self.mu = mu @@ -146,4 +159,40 @@ def fit(self, data): data : ndarray (n, 3) Vector data the distribution is fitted to ''' - self.mu, self.kappa = _fit(data) \ No newline at end of file + self.mu, self.kappa = _fit(data) + + def angle_within_probability_mass(self, alpha, deg = False): + r""" + Calculates the angle which contains probability mass alpha of the VMF density around the mean angle + + Reference: Fayat, 2021. Conversion of the von Mises-Fisher concentration parameter to an equivalent angle. + + https://github.com/rfayat/SphereProba/blob/main/ressources/vmf_integration.pdf + + Arguments + ---------- + alpha : float + Probability mass. Must be :math:`0<\alpha<1` + deg : optional, default False + If True, converts the result into degrees + + Returns + ---------- + angle : float + Resulting angle + + """ + if self.kappa is not None: + + nominator = np.log(1-alpha + alpha * np.exp(-2 * self.kappa)) + + angle = np.arccos(1+nominator/self.kappa) + + if deg == True: + angle=np.rad2deg(angle) + + return angle + + else: + + raise ValueError("Concentration parameter kappa unknown.") \ No newline at end of file