diff --git a/include/xtensor/core/xmath.hpp b/include/xtensor/core/xmath.hpp
index 6b86c4b74..cca2e2a91 100644
--- a/include/xtensor/core/xmath.hpp
+++ b/include/xtensor/core/xmath.hpp
@@ -1825,11 +1825,12 @@ namespace xt
     namespace detail
     {
         template <class T, class S, class ST>
-        inline auto mean_division(S&& s, ST e_size)
+        inline auto mean_division(S&& s, ST e_size, ST ddof)
         {
             using value_type = typename std::conditional_t<std::is_same<T, void>::value, double, T>;
             // Avoids floating point exception when s.size is 0
-            value_type div = s.size() != ST(0) ? static_cast<value_type>(e_size / s.size()) : value_type(0);
+            value_type div = s.size() != ST(0) ? static_cast<value_type>(e_size / s.size() - ddof)
+                                               : value_type(0);
             return std::move(s) / std::move(div);
         }
 
@@ -1848,7 +1849,7 @@ namespace xt
             const size_type size = e.size();
             XTENSOR_ASSERT(static_cast<size_type>(ddof) <= size);
             auto s = sum<T>(std::forward<E>(e), std::forward<X>(axes), es);
-            return mean_division<T>(std::move(s), size - static_cast<size_type>(ddof));
+            return mean_division<T>(std::move(s), size, static_cast<size_type>(ddof));
         }
 
         template <class T, class E, class I, std::size_t N, class D, class EVS>
@@ -1858,7 +1859,7 @@ namespace xt
             const size_type size = e.size();
             XTENSOR_ASSERT(static_cast<size_type>(ddof) <= size);
             auto s = sum<T>(std::forward<E>(e), axes, es);
-            return mean_division<T>(std::move(s), size - static_cast<size_type>(ddof));
+            return mean_division<T>(std::move(s), size, static_cast<size_type>(ddof));
         }
 
         template <class T, class E, class D, class EVS, XTL_REQUIRES(is_reducer_options<EVS>, xtl::is_integral<D>)>
diff --git a/test/test_xmath.cpp b/test/test_xmath.cpp
index 57771859b..4592a07cd 100644
--- a/test/test_xmath.cpp
+++ b/test/test_xmath.cpp
@@ -885,6 +885,18 @@ namespace xt
         EXPECT_EQ(xt::average(v, w, {0, 1})(), m);
     }
 
+    TEST(xmath, variance_ddof)
+    {
+        xt::xtensor<double, 2> v = {{1.0, 3.0, 5.0}, {2.0, 4.0, 6.0}};
+        xt::xtensor<double, 1> v1_ddof1 = {4.0, 4.0};
+        xt::xtensor<double, 1> v1_ddof2 = {8.0, 8.0};
+        std::vector<std::size_t> axes = {1};
+
+        EXPECT_TRUE(xt::all(xt::equal(xt::variance(v, {1}, 1), v1_ddof1)));
+        EXPECT_TRUE(xt::all(xt::equal(xt::variance(v, {1}, 2), v1_ddof2)));
+        EXPECT_TRUE(xt::all(xt::equal(xt::variance(v, axes, 2), v1_ddof2)));
+    }
+
     /************************
      * Linear interpolation *
      ************************/