Central moment (#53)

* Added first version of central_moment, variance and standard_deviation

* Bumped visual studio version for tests

---------

Co-authored-by: Florian Weik <[email protected]>

README.md

@@ -54,6 +54,38 @@ The mean is the 1st (raw) [moment](#nrngmoment).
 auto const mean = nrng::mean(std::array{3., 1., 4., 1., 5., 9.});
 ```
 
+#### [nrng::central_moment](include/nrng/central_moment.hpp)
+
+The n-th (raw) central moment $M_n$ of a sample $S$, defined as
+$$M_n = \frac{1}{|S|} \sum_{ x_i \in S } (x_i - \mu)^n. $$
+
+where $\mu_S$ is the mean of $S$. 
+
+The order of summation is arbitraty, e.g. it behaves like `reduce` and not like `accumulate`. If the sample is empty,
+the moments are undefined.
+
+Example Usage
+
+```c++
+auto const snd = nrng::central_moment<2>(std::array{1, 2, 3, 4, 5});
+```
+
+#### [nrng::variance](include/nrng/central_moment.hpp)
+
+The variance is the 2nd (raw) [central_moment](#nrng---centralmoment).
+
+```c++
+auto const mean = nrng::variance(std::array{3., 1., 4., 1., 5., 9.});
+```
+
+#### [nrng::standard_deviation](include/nrng/central_moment.hpp)
+
+The standard deviation is the square root of the [variance](#nrng---variance).
+
+```c++
+auto const mean = nrng::standard_deviation(std::array{3., 1., 4., 1., 5., 9.});
+```
+
 ### Floating-point algorithms
 
 #### [nrng::kahan_plus](include/nrng/kahan.hpp)

appveyor.yml

@@ -1,5 +1,5 @@
 image:
-  - Visual Studio 2019
+  - Visual Studio 2022
 clone_folder: c:\projects\source
 
 build_script:
@@ -8,7 +8,7 @@ build_script:
 
     cd build
 
-    cmake c:\projects\source -G "Visual Studio 16 2019" -DCMAKE_BUILD_TYPE:STRING=Release
+    cmake c:\projects\source -G "Visual Studio 17 2022" -DCMAKE_BUILD_TYPE:STRING=Release
 
     cmake --build . --config "Release"
 

include/nrng/central_moment.hpp

@@ -0,0 +1,64 @@
+#ifndef NRNG_CENTRAL_MOMENT_HPP
+#define NRNG_CENTRAL_MOMENT_HPP
+
+#include <nrng/moment.hpp>
+#include <nrng/type_traits.hpp>
+
+#include <cmath>
+
+namespace nrng {
+/**
+ * @brief N-th central moment of range.
+ */
+template <unsigned Exponent, std::input_iterator I, std::sentinel_for<I> S>
+constexpr auto central_moment(I first, S last) {
+  auto const mean = ::nrng::mean(first, last);
+
+  return ::nrng::transform_reduce(
+             first, last, std::iter_value_t<I>{}, std::plus<>{},
+             [mean](auto v) { return power<Exponent>(v - mean); }) /
+         std::ranges::distance(first, last);
+}
+
+/**
+ * @brief N-th central moment of range.
+ */
+template <unsigned Exponent, std::ranges::input_range Rng>
+constexpr auto central_moment(Rng values) {
+  return central_moment<Exponent>(std::ranges::begin(values),
+                                  std::ranges::end(values));
+}
+
+template <std::ranges::input_range Rng> constexpr auto variance(Rng values) {
+  return central_moment<2>(values);
+}
+
+template <std::input_iterator I, std::sentinel_for<I> S>
+constexpr auto variance(I first, S last) {
+  return central_moment<2>(first, last);
+}
+
+struct sqrt {
+  template <class T>
+  requires std::is_arithmetic_v<T>
+  auto operator()(T const x) const {
+    return std::sqrt(x);
+  }
+};
+
+template <std::input_iterator I, std::sentinel_for<I> S>
+requires unary_closed_under<std::iter_value_t<I>, sqrt>
+constexpr auto standard_deviation(I first, S last) {
+  return sqrt{}(variance(first, last));
+}
+
+template <std::ranges::input_range Rng>
+requires unary_closed_under<std::ranges::range_value_t<Rng>,
+                            sqrt>
+constexpr auto standard_deviation(Rng values) {
+  return standard_deviation(std::ranges::begin(values),
+                            std::ranges::end(values));
+}
+} // namespace nrng
+
+#endif

tests/central_moment.test.cpp

@@ -0,0 +1,37 @@
+#include <doctest/doctest.h>
+
+#include <nrng/central_moment.hpp>
+
+#include <array>
+
+using nrng::moment;
+using nrng::central_moment;
+using nrng::mean;
+using nrng::variance;
+using nrng::standard_deviation;
+
+TEST_CASE("first central moment") {
+  auto const values = std::array{1, 2, 3, 4, 5};
+
+  CHECK(central_moment<1>(values) == 0);
+  CHECK(central_moment<1>(values.begin(), values.end()) == central_moment<1>(values));
+}
+
+TEST_CASE("second central moment") {
+  auto const values = std::array{1, 2, 3, 4, 5};
+  auto shift_by_mean = [mu = mean(values)](auto v) { return v - mu; };
+
+  CHECK(central_moment<2>(values) == moment<2>(values | std::views::transform(shift_by_mean)));
+}
+
+TEST_CASE("variance") {
+  auto const values = std::array{1, 2, 3, 4, 5};
+
+  CHECK(variance(values) == central_moment<2>(values));
+}
+
+TEST_CASE("standard deviation") {
+  auto const values = std::array{1, 2, 3, 4, 5};
+
+  CHECK(standard_deviation(values) == std::sqrt(variance(values)));
+}

tests/kahan.test.cpp

@@ -6,7 +6,7 @@
 
 TEST_CASE("kahan_plus result test") {
   /* Test numbers, chosen such that there are not
-   * enought significant places in float for direct
+   * enough significant places in float for direct
    * summation to yield the correct result. */
   auto const init = 1000000.0F;
   auto const t1 = 3.14159F;