From 0852c16a2a23f25e284b796a9f69cf98ff86bea9 Mon Sep 17 00:00:00 2001 From: Rasmus Henningsson Date: Mon, 19 Jun 2023 15:18:15 +0800 Subject: [PATCH] median_absolute_deviation bug fix for non-zero center (#997) Fix bug causing median_absolute_deviation to return incorrect values when a non-zero center (such as median which is the default) was used. --- .../math/statistics/univariate_statistics.hpp | 8 +-- ...e_statistics_backwards_compatible_test.cpp | 65 +++++++++++++++++-- test/univariate_statistics_test.cpp | 65 +++++++++++++++++-- 3 files changed, 124 insertions(+), 14 deletions(-) diff --git a/include/boost/math/statistics/univariate_statistics.hpp b/include/boost/math/statistics/univariate_statistics.hpp index 082517d650..898ecb0c7b 100644 --- a/include/boost/math/statistics/univariate_statistics.hpp +++ b/include/boost/math/statistics/univariate_statistics.hpp @@ -505,14 +505,14 @@ auto median_absolute_deviation(ExecutionPolicy&& exec, RandomAccessIterator firs { auto middle = first + (num_elems - 1)/2; std::nth_element(exec, first, middle, last, comparator); - return abs(*middle); + return abs(*middle-center); } else { auto middle = first + num_elems/2 - 1; std::nth_element(exec, first, middle, last, comparator); std::nth_element(exec, middle, middle+1, last, comparator); - return (abs(*middle) + abs(*(middle+1)))/abs(static_cast(2)); + return (abs(*middle-center) + abs(*(middle+1)-center))/abs(static_cast(2)); } } @@ -1043,14 +1043,14 @@ Real median_absolute_deviation(RandomAccessIterator first, RandomAccessIterator { auto middle = first + (num_elems - 1)/2; std::nth_element(first, middle, last, comparator); - return abs(*middle); + return abs(*middle-center); } else { auto middle = first + num_elems/2 - 1; std::nth_element(first, middle, last, comparator); std::nth_element(middle, middle+1, last, comparator); - return (abs(*middle) + abs(*(middle+1)))/abs(static_cast(2)); + return (abs(*middle-center) + abs(*(middle+1)-center))/abs(static_cast(2)); } } diff --git a/test/univariate_statistics_backwards_compatible_test.cpp b/test/univariate_statistics_backwards_compatible_test.cpp index dfb8f30c26..9c50d6f2dc 100644 --- a/test/univariate_statistics_backwards_compatible_test.cpp +++ b/test/univariate_statistics_backwards_compatible_test.cpp @@ -612,11 +612,7 @@ void test_median_absolute_deviation() v = {-1, 1}; m = boost::math::statistics::median_absolute_deviation(v.begin(), v.end(), 0); BOOST_TEST_EQ(m, 1); - // The median is zero, so coincides with the default: - m = boost::math::statistics::median_absolute_deviation(v.begin(), v.end()); - BOOST_TEST_EQ(m, 1); - - m = boost::math::statistics::median_absolute_deviation(v); + m = boost::math::statistics::median_absolute_deviation(v, 0); BOOST_TEST_EQ(m, 1); @@ -649,6 +645,65 @@ void test_median_absolute_deviation() u[5] = -3; m = boost::math::statistics::median_absolute_deviation(u, 0); BOOST_TEST_EQ(m, 2); + + + v = {-1, 2, -3, 4, -5, 6, -7}; + m = boost::math::statistics::median_absolute_deviation(v.begin(), v.end()); + BOOST_TEST_EQ(m, 4); + + g = std::mt19937(12); + std::shuffle(v.begin(), v.end(), g); + m = boost::math::statistics::median_absolute_deviation(v); + BOOST_TEST_EQ(m, 4); + + v = {1, -2, -3, 3, -4, -5}; + m = boost::math::statistics::median_absolute_deviation(v.begin(), v.end()); + BOOST_TEST_EQ(m, 2); + std::shuffle(v.begin(), v.end(), g); + m = boost::math::statistics::median_absolute_deviation(v.begin(), v.end()); + BOOST_TEST_EQ(m, 2); + + v = {-1}; + m = boost::math::statistics::median_absolute_deviation(v.begin(), v.end()); + BOOST_TEST_EQ(m, 0); + + v = {-1, 1}; + m = boost::math::statistics::median_absolute_deviation(v.begin(), v.end()); + BOOST_TEST_EQ(m, 1); + + m = boost::math::statistics::median_absolute_deviation(v); + BOOST_TEST_EQ(m, 1); + + + v = {2, -4}; + m = boost::math::statistics::median_absolute_deviation(v.begin(), v.end()); + BOOST_TEST_EQ(m, 3); + + v = {1, -1, 1}; + m = boost::math::statistics::median_absolute_deviation(v.begin(), v.end()); + BOOST_TEST_EQ(m, 0); + + v = {1, 2, -3}; + m = boost::math::statistics::median_absolute_deviation(v.begin(), v.end()); + BOOST_TEST_EQ(m, 1); + std::shuffle(v.begin(), v.end(), g); + m = boost::math::statistics::median_absolute_deviation(v.begin(), v.end()); + BOOST_TEST_EQ(m, 1); + + w = {1, 2, -3}; + m = boost::math::statistics::median_absolute_deviation(w); + BOOST_TEST_EQ(m, 1); + + // boost.ublas vector? + boost::numeric::ublas::vector u2(6); + u2[0] = 1; + u2[1] = 2; + u2[2] = -3; + u2[3] = 1; + u2[4] = 2; + u2[5] = -3; + m = boost::math::statistics::median_absolute_deviation(u2); + BOOST_TEST_EQ(m, 1); } diff --git a/test/univariate_statistics_test.cpp b/test/univariate_statistics_test.cpp index 68cfc71621..77b8e8fc09 100644 --- a/test/univariate_statistics_test.cpp +++ b/test/univariate_statistics_test.cpp @@ -577,11 +577,7 @@ void test_median_absolute_deviation(ExecutionPolicy&& exec) v = {-1, 1}; m = boost::math::statistics::median_absolute_deviation(exec, v.begin(), v.end(), 0); BOOST_TEST_EQ(m, 1); - // The median is zero, so coincides with the default: - m = boost::math::statistics::median_absolute_deviation(exec, v.begin(), v.end()); - BOOST_TEST_EQ(m, 1); - - m = boost::math::statistics::median_absolute_deviation(exec, v); + m = boost::math::statistics::median_absolute_deviation(exec, v, 0); BOOST_TEST_EQ(m, 1); @@ -614,6 +610,65 @@ void test_median_absolute_deviation(ExecutionPolicy&& exec) u[5] = -3; m = boost::math::statistics::median_absolute_deviation(exec, u, 0); BOOST_TEST_EQ(m, 2); + + + v = {-1, 2, -3, 4, -5, 6, -7}; + m = boost::math::statistics::median_absolute_deviation(exec, v.begin(), v.end()); + BOOST_TEST_EQ(m, 4); + + g = std::mt19937(12); + std::shuffle(v.begin(), v.end(), g); + m = boost::math::statistics::median_absolute_deviation(exec, v); + BOOST_TEST_EQ(m, 4); + + v = {1, -2, -3, 3, -4, -5}; + m = boost::math::statistics::median_absolute_deviation(exec, v.begin(), v.end()); + BOOST_TEST_EQ(m, 2); + std::shuffle(v.begin(), v.end(), g); + m = boost::math::statistics::median_absolute_deviation(exec, v.begin(), v.end()); + BOOST_TEST_EQ(m, 2); + + v = {-1}; + m = boost::math::statistics::median_absolute_deviation(exec, v.begin(), v.end()); + BOOST_TEST_EQ(m, 0); + + v = {-1, 1}; + m = boost::math::statistics::median_absolute_deviation(exec, v.begin(), v.end()); + BOOST_TEST_EQ(m, 1); + + m = boost::math::statistics::median_absolute_deviation(exec, v); + BOOST_TEST_EQ(m, 1); + + + v = {2, -4}; + m = boost::math::statistics::median_absolute_deviation(exec, v.begin(), v.end()); + BOOST_TEST_EQ(m, 3); + + v = {1, -1, 1}; + m = boost::math::statistics::median_absolute_deviation(exec, v.begin(), v.end()); + BOOST_TEST_EQ(m, 0); + + v = {1, 2, -3}; + m = boost::math::statistics::median_absolute_deviation(exec, v.begin(), v.end()); + BOOST_TEST_EQ(m, 1); + std::shuffle(v.begin(), v.end(), g); + m = boost::math::statistics::median_absolute_deviation(exec, v.begin(), v.end()); + BOOST_TEST_EQ(m, 1); + + w = {1, 2, -3}; + m = boost::math::statistics::median_absolute_deviation(exec, w); + BOOST_TEST_EQ(m, 1); + + // boost.ublas vector? + boost::numeric::ublas::vector u2(6); + u2[0] = 1; + u2[1] = 2; + u2[2] = -3; + u2[3] = 1; + u2[4] = 2; + u2[5] = -3; + m = boost::math::statistics::median_absolute_deviation(exec, u2); + BOOST_TEST_EQ(m, 1); }