add: wrapped stable distribution; bump: v0.1.6

circstat · Dec 6, 2024 · 11b1203 · 11b1203
1 parent c2d7d41
commit 11b1203
Show file tree

Hide file tree

Showing 11 changed files with 323 additions and 48 deletions.
diff --git a/docs/docs/images/circ-mod-wrapstable.png b/docs/docs/images/circ-mod-wrapstable.png
diff --git a/docs/mkdocs.yml b/docs/mkdocs.yml
@@ -2,7 +2,7 @@ site_name: pycircstat2
 theme:
     name: material
     features:
-        - navigation.instant
+        - navigation.sections
         - content.code.copy
         - content.tabs.link
         - content.mathjax

diff --git a/examples/B1-Fisher-1993.ipynb b/examples/B1-Fisher-1993.ipynb
diff --git a/examples/B2-Zar-2010.ipynb b/examples/B2-Zar-2010.ipynb
@@ -889,7 +889,7 @@
     {
      "data": {
       "text/plain": [
-       "[<matplotlib.lines.Line2D at 0x7f8b4561f400>]"
+       "[<matplotlib.lines.Line2D at 0x7f8c4bd51c30>]"
       ]
      },
      "execution_count": 18,
@@ -1447,19 +1447,17 @@
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "Last updated: 2024-11-21 16:44:05CET\n",
+      "Last updated: 2024-12-06 14:58:21CET\n",
       "\n",
       "Python implementation: CPython\n",
       "Python version       : 3.10.13\n",
       "IPython version      : 8.29.0\n",
       "\n",
-      "pycircstat2: 0.1.0\n",
-      "\n",
-      "matplotlib : 3.9.2\n",
-      "polars     : 1.14.0\n",
       "scipy      : 1.14.1\n",
-      "pycircstat2: 0.1.0\n",
+      "matplotlib : 3.9.2\n",
+      "pycircstat2: 0.1.6\n",
       "numpy      : 2.1.3\n",
+      "polars     : 1.14.0\n",
       "\n",
       "Watermark: 2.5.0\n",
       "\n"
@@ -1468,7 +1466,7 @@
    ],
    "source": [
     "%load_ext watermark\n",
-    "%watermark --time --date --timezone --updated --python --iversions --watermark -p pycircstat2"
+    "%watermark --time --date --timezone --updated --python --iversions --watermark"
    ]
   }
  ],

diff --git a/examples/B3-Pewsey-2014.ipynb b/examples/B3-Pewsey-2014.ipynb
@@ -836,7 +836,7 @@
     {
      "data": {
       "text/plain": [
-       "<matplotlib.legend.Legend at 0x7fbfab3739d0>"
+       "<matplotlib.legend.Legend at 0x7fd682262f80>"
       ]
      },
      "execution_count": 19,
@@ -1254,18 +1254,16 @@
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "Last updated: 2024-12-05 11:51:22CET\n",
+      "Last updated: 2024-12-06 14:58:35CET\n",
       "\n",
       "Python implementation: CPython\n",
       "Python version       : 3.10.13\n",
       "IPython version      : 8.29.0\n",
       "\n",
-      "pycircstat2: 0.1.4\n",
-      "\n",
-      "matplotlib : 3.9.2\n",
-      "pycircstat2: 0.1.4\n",
       "numpy      : 2.1.3\n",
       "polars     : 1.14.0\n",
+      "pycircstat2: 0.1.6\n",
+      "matplotlib : 3.9.2\n",
       "\n",
       "Watermark: 2.5.0\n",
       "\n"
@@ -1274,7 +1272,7 @@
    ],
    "source": [
     "%load_ext watermark\n",
-    "%watermark --time --date --timezone --updated --python --iversions --watermark -p pycircstat2"
+    "%watermark --time --date --timezone --updated --python --iversions --watermark"
    ]
   }
  ],

diff --git a/examples/T0-utils.ipynb b/examples/T0-utils.ipynb
@@ -305,16 +305,14 @@
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "Last updated: 2024-11-21 16:44:20CET\n",
+      "Last updated: 2024-12-06 14:57:40CET\n",
       "\n",
       "Python implementation: CPython\n",
       "Python version       : 3.10.13\n",
       "IPython version      : 8.29.0\n",
       "\n",
-      "pycircstat2: 0.1.0\n",
-      "\n",
       "numpy      : 2.1.3\n",
-      "pycircstat2: 0.1.0\n",
+      "pycircstat2: 0.1.6\n",
       "\n",
       "Watermark: 2.5.0\n",
       "\n"
@@ -323,7 +321,7 @@
    ],
    "source": [
     "%load_ext watermark\n",
-    "%watermark --time --date --timezone --updated --python --iversions --watermark -p pycircstat2"
+    "%watermark --time --date --timezone --updated --python --iversions --watermark"
    ]
   }
  ],

diff --git a/examples/T1-descriptive-statistics.ipynb b/examples/T1-descriptive-statistics.ipynb
@@ -420,7 +420,7 @@
       "median                    279.0                     247.5                     245.0               \n",
       "mean                      280.8                     248.7                     247.6               \n",
       "95% median CI             [245. 315.]               [229. 277.]               [229. 267.]         \n",
-      "95% bootstrap mean CI     [265.1 298.5]             [228.1 269.4]             [229.9 262. ]       \n",
+      "95% bootstrap mean CI     [263.1 298.7]             [218.  264.8]             [231.  260.2]       \n",
       "95% large-sample mean CI  -                         -                         [232.7 262.5]       \n"
      ]
     }
@@ -481,16 +481,14 @@
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "Last updated: 2024-11-28 14:31:06CET\n",
+      "Last updated: 2024-12-06 14:57:57CET\n",
       "\n",
       "Python implementation: CPython\n",
       "Python version       : 3.10.13\n",
       "IPython version      : 8.29.0\n",
       "\n",
-      "pycircstat2: 0.1.2\n",
-      "\n",
       "numpy      : 2.1.3\n",
-      "pycircstat2: 0.1.2\n",
+      "pycircstat2: 0.1.6\n",
       "\n",
       "Watermark: 2.5.0\n",
       "\n"
@@ -499,7 +497,7 @@
    ],
    "source": [
     "%load_ext watermark\n",
-    "%watermark --time --date --timezone --updated --python --iversions --watermark -p pycircstat2"
+    "%watermark --time --date --timezone --updated --python --iversions --watermark"
    ]
   }
  ],

diff --git a/examples/T2-hypothesis-testing.ipynb b/examples/T2-hypothesis-testing.ipynb
@@ -110,7 +110,7 @@
       "Modified Rayleigh's Test of Uniformity\n",
       "--------------------------------------\n",
       "H0: ρ = 0\n",
-      "HA: ρ ≠ 0 and μ = 1.57080 rad\n",
+      "HA: ρ ≠ 0 and μ = {angle:.5f} rad\n",
       "\n",
       "Test Statistics: 9.49761\n",
       "P-value: 0.00001 ***\n"
@@ -667,15 +667,13 @@
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "Last updated: 2024-11-21 16:44:43CET\n",
+      "Last updated: 2024-12-06 14:58:46CET\n",
       "\n",
       "Python implementation: CPython\n",
       "Python version       : 3.10.13\n",
       "IPython version      : 8.29.0\n",
       "\n",
-      "pycircstat2: 0.1.0\n",
-      "\n",
-      "pycircstat2: 0.1.0\n",
+      "pycircstat2: 0.1.6\n",
       "numpy      : 2.1.3\n",
       "\n",
       "Watermark: 2.5.0\n",
@@ -685,7 +683,7 @@
    ],
    "source": [
     "%load_ext watermark\n",
-    "%watermark --time --date --timezone --updated --python --iversions --watermark -p pycircstat2"
+    "%watermark --time --date --timezone --updated --python --iversions --watermark"
    ]
   }
  ],

diff --git a/examples/T3-circular-models.ipynb b/examples/T3-circular-models.ipynb
diff --git a/pycircstat2/distributions.py b/pycircstat2/distributions.py
@@ -18,6 +18,7 @@
     "jonespewsey_sineskewed",
     "jonespewsey_asym",
     "inverse_batschelet",
+    "wrapstable",
 ]
 
 
@@ -1699,3 +1700,113 @@ def kernel(x):
 
         c = 1 / (mult * (con1 - con2 * intval))
     return c
+
+
+class wrapstable_gen(rv_continuous):
+    r"""
+    Wrapped Stable Distribution
+
+    - is symmetric around $\delta$ when $\beta = 0$, and to be skewed to the right (left) if $\beta > 0$ ($\beta < 0$).
+    - can be reduced to
+        - the wrapped normal distribution when $\alpha = 2, \beta = 0$.
+        - the wrapped Cauchy distribution when $\alpha = 1, \beta = 0$.
+        - the wrappd Lévy distribution when $\alpha = 1/2, \beta = 1$
+
+    ![wrapstable](../images/circ-mod-wrapstable.png)
+
+    References
+    ----------
+    - Pewsey, A. (2008). The wrapped stable family of distributions as a flexible model for circular data. Computational Statistics & Data Analysis, 52(3), 1516-1523.
+    """
+
+    def _validate_params(self, delta, alpha, beta, gamma):
+        return (
+            (0 <= delta <= np.pi * 2)
+            and (0 < alpha <= 2)
+            and (-1 < beta < 1)
+            and (gamma > 0)
+        )
+
+    def _argcheck(self, delta, alpha, beta, gamma):
+        if self._validate_params(delta, alpha, beta, gamma):
+            return True
+        else:
+            return False
+
+    def _pdf(self, x, delta, alpha, beta, gamma):
+
+        def rho_p(p, alpha, gamma):
+            return np.exp(-((gamma * p) ** alpha))
+
+        def mu_p(p, alpha, beta, gamma, delta):
+            if np.all(alpha == 1):
+                mu = delta * p - 2 * beta * gamma * p * np.log(gamma * p) / np.pi
+            else:
+                mu = delta * p + beta * np.tan(np.pi * alpha / 2) * (
+                    (gamma * p) ** alpha - gamma * p
+                )
+            return mu
+
+        series_sum = 0
+        for p in np.arange(1, 150):
+            rho = rho_p(p, alpha, gamma)
+            mu = mu_p(p, alpha, beta, gamma, delta)
+            series_sum += rho * np.cos(p * x - mu)
+
+        pdf_values = 1 / (2 * np.pi) * (1 + 2 * series_sum)
+
+        return pdf_values
+
+    def pdf(self, x, delta, alpha, beta, gamma, *args, **kwargs):
+        r"""
+        Probability density function of the Wrapped Stable distribution.
+
+        $$
+        f(\theta) = \frac{1}{2\pi} \left[1 + 2 \sum_{p=1}^{\infty} \rho_p \cos\left(p(\theta - \mu_p)\right)\right]
+        $$
+
+        , where $\rho_p$ is the $p$th mean resultant length and $\mu_p$ is the $p$th mean direction:
+
+        $$
+        \rho_p = \exp\left(-(\gamma p)^\alpha\right)
+        $$
+
+        $$
+        \mu_p = 
+        \begin{cases}
+            \delta p + \beta \tan\left(\frac{\pi \alpha}{2}\right) \left((\gamma p)^\alpha - \gamma p\right), & \alpha \neq 1 \\
+            \delta p - \beta \frac{2}{\pi} \log(\gamma p), & \text{if } \alpha = 1
+        \end{cases}
+        $$
+
+        Parameters
+        ----------
+        x : array_like
+            Points at which to evaluate the PDF, defined on the interval $[0, 2\pi)$.
+        delta : float
+            Location parameter, $0 \leq \delta \leq 2\pi$. This is the mean direction of the distribution.
+        alpha : float
+            Stability parameter, $0 < \alpha \leq 2$. Higher values indicate heavier tails.
+        beta : float
+            Skewness parameter, $-1 < \beta < 1$. Controls the asymmetry of the distribution.
+        gamma : float
+            Scale parameter, $\gamma > 0$. Scales the distribution.
+
+        Returns
+        -------
+        pdf_values : array_like
+            Values of the probability density function at the specified points.
+        """
+        return super().pdf(x, delta, alpha, beta, gamma, *args, **kwargs)
+
+    def _cdf(self, x, delta, alpha, beta, gamma):
+
+        @np.vectorize
+        def _cdf_single(x, delta, alpha, beta, gamma):
+            integral, _ = quad(self._pdf, a=0, b=x, args=(delta, alpha, beta, gamma))
+            return integral
+
+        return _cdf_single(x, delta, alpha, beta, gamma)
+
+
+wrapstable = wrapstable_gen(name="wrapstable")
diff --git a/pyproject.toml b/pyproject.toml
@@ -4,7 +4,7 @@ build-backend = "setuptools.build_meta"
 
 [project]
 name = "pycircstat2"
-version = "0.1.5"  
+version = "0.1.6"  
 description = "Circular statistcs with Python."
 authors = [
 ]