diff --git a/include/xsf/stats.h b/include/xsf/stats.h
index 61bf7cfcd6..530df4e821 100644
--- a/include/xsf/stats.h
+++ b/include/xsf/stats.h
@@ -3,6 +3,7 @@
 #include "xsf/bessel.h"
 #include "xsf/cephes/bdtr.h"
 #include "xsf/cephes/chdtr.h"
+#include "xsf/cephes/const.h"
 #include "xsf/cephes/fdtr.h"
 #include "xsf/cephes/gdtr.h"
 #include "xsf/cephes/incbet.h"
@@ -300,4 +301,70 @@ inline float tukeylambdacdf(float x, double lmbda) {
     return tukeylambdacdf(static_cast<double>(x), static_cast<double>(lmbda));
 }
 
+inline double von_mises_cdf_series(double k, double x, unsigned int p) {
+    double s, c, sn, cn, R, V;
+    unsigned int n;
+    s = std::sin(x);
+    c = std::cos(x);
+    sn = std::sin(p * x);
+    cn = std::cos(p * x);
+    R = 0;
+    V = 0;
+    for (n = p - 1; n > 0; --n) {
+        double sn_new = sn * c - cn * s;
+        double cn_new = cn * c + sn * s;
+        sn = sn_new;
+        cn = cn_new;
+        R = k / (2 * n + k * R);
+        V = R * (sn / n + V);
+    }
+    return 0.5 + x / (2.0 * M_PI) + V / M_PI;
+}
+
+inline double von_mises_cdf_normalapprox(double k, double x) {
+    double b = xsf::cephes::detail::SQRT2OPI / cephes::i0e(k); // Check for negative k
+    double z = b * std::sin(x / 2.0);
+    return ndtr(z);
+}
+
+inline std::vector<double> von_mises_cdf(const std::vector<double> &k_obj, const std::vector<double> &x_obj) {
+    // Compute the CDF of the von Mises distribution with concentration parameter k at angle x
+
+    if (k_obj.size() != 1 && x_obj.size() != 1 && k_obj.size() != x_obj.size()) {
+        throw std::invalid_argument("Incompatible sizes for broadcasting.");
+    }
+
+    bool zerodim = (k_obj.size() == 1 && x_obj.size() == 1);
+    size_t n = std::max(k_obj.size(), x_obj.size());
+    std::vector<double> k(n), x(n);
+    for (size_t i = 0; i < n; ++i) {
+        k[i] = k_obj.size() == 1 ? k_obj[0] : k_obj[i];
+        x[i] = x_obj.size() == 1 ? x_obj[0] : x_obj[i];
+    }
+    std::vector<double> ix(n);
+    for (size_t i = 0; i < n; ++i) {
+        ix[i] = std::round(x[i] / (2 * M_PI));
+        x[i] -= ix[i] * 2.0 * M_PI;
+    }
+
+    // These values should give 12 decimal digits
+    const double CK = 50.0;
+    const double a1 = 28.0, a2 = 0.5, a3 = 100.0, a4 = 5.0;
+    std::vector<double> result(n);
+    for (size_t i = 0; i < n; ++i) {
+        if (k[i] < CK) {
+            unsigned int p = static_cast<unsigned int>(1 + a1 + a2 * k[i] - a3 / (k[i] + a4));
+            double val = von_mises_cdf_series(k[i], x[i], p);
+            val = std::min(std::max(val, 0.0), 1.0);
+            result[i] = val;
+        } else {
+            result[i] = von_mises_cdf_normalapprox(k[i], x[i]);
+        }
+    }
+    for (size_t i = 0; i < n; ++i) {
+        result[i] += ix[i];
+    }
+    return result;
+}
+
 } // namespace xsf
diff --git a/tests/xsf_tests/test_von_mises_cdf.cpp b/tests/xsf_tests/test_von_mises_cdf.cpp
new file mode 100644
index 0000000000..a05269a19d
--- /dev/null
+++ b/tests/xsf_tests/test_von_mises_cdf.cpp
@@ -0,0 +1,40 @@
+#include "../testing_utils.h"
+#include <xsf/stats.h>
+
+TEST_CASE("von_mises_cdf test", "[von_mises_cdf][xsf_tests]") {
+    // Reference values computed with scipy.stats.von_mises_cdf
+    // import numpy as np
+    // from scipy.stats._stats import von_mises_cdf
+
+    // np.set_printoptions(precision=20)
+    // k_obj = np.linspace(1e-3, 3.0, 10)
+    // x_obj = np.linspace(-10.0, 10.0, 10)
+    // von_mises_cdf(k_obj, x_obj)
+
+    SECTION("vector inputs (SciPy reference values)") {
+        const size_t n = 10;
+        const std::vector<double> k_obj = linspace(1e-3, 3.0, n);
+        const std::vector<double> x_obj = linspace(-10.0, 10.0, n);
+        const std::vector<double> expected = {-1.0914628654411138,   -0.7904352686647403, -0.3085050816322099,
+                                              -0.008913110513925071, 0.1450905974469251,  0.8857870521643215,
+                                              1.0018330872501384,    1.1579435355869516,  1.9789394168440955,
+                                              2.0011107464769506};
+
+        const double rtol = 1e-8;
+        const std::vector<double> result = xsf::von_mises_cdf(k_obj, x_obj);
+        for (size_t i = 0; i < n; ++i) {
+            const double ref = expected[i];
+            const auto rel_error = xsf::extended_relative_error(result[i], ref);
+            CAPTURE(i, k_obj[i], x_obj[i], result[i], ref, rtol, rel_error);
+            REQUIRE(rel_error <= rtol);
+        }
+    }
+
+    SECTION("broadcast scalar k") {
+        std::vector<double> k = {1.0};
+        std::vector<double> x = linspace(-5.0, 5.0, 10);
+
+        auto res = xsf::von_mises_cdf(k, x);
+        REQUIRE(res.size() == x.size());
+    }
+}