feat():initial version

install/boost_1_75_0/include/boost/math/statistics/anderson_darling.hpp

@@ -0,0 +1,112 @@
+/*
+ * Copyright Nick Thompson, 2019
+ * Use, modification and distribution are subject to the
+ * Boost Software License, Version 1.0. (See accompanying file
+ * LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+ */
+
+#ifndef BOOST_MATH_STATISTICS_ANDERSON_DARLING_HPP
+#define BOOST_MATH_STATISTICS_ANDERSON_DARLING_HPP
+
+#include <cmath>
+#include <algorithm>
+#include <boost/math/statistics/univariate_statistics.hpp>
+#include <boost/math/special_functions/erf.hpp>
+
+namespace boost { namespace math { namespace statistics {
+
+template<class RandomAccessContainer>
+auto anderson_darling_normality_statistic(RandomAccessContainer const & v,
+                                          typename RandomAccessContainer::value_type mu = std::numeric_limits<typename RandomAccessContainer::value_type>::quiet_NaN(),
+                                          typename RandomAccessContainer::value_type sd = std::numeric_limits<typename RandomAccessContainer::value_type>::quiet_NaN())
+{
+    using Real = typename RandomAccessContainer::value_type;
+    using std::log;
+    using std::sqrt;
+    using boost::math::erfc;
+
+    if (std::isnan(mu)) {
+        mu = boost::math::statistics::mean(v);
+    }
+    if (std::isnan(sd)) {
+        sd = sqrt(boost::math::statistics::sample_variance(v));
+    }
+
+    typedef boost::math::policies::policy<
+          boost::math::policies::promote_float<false>,
+          boost::math::policies::promote_double<false> >
+          no_promote_policy;
+
+    // This is where Knuth's literate programming could really come in handy!
+    // I need some LaTeX. The idea is that before any observation, the ecdf is identically zero.
+    // So we need to compute:
+    // \int_{-\infty}^{v_0} \frac{F(x)F'(x)}{1- F(x)} \, \mathrm{d}x, where F(x) := \frac{1}{2}[1+\erf(\frac{x-\mu}{\sigma \sqrt{2}})]
+    // Astonishingly, there is an analytic evaluation to this integral, as you can validate with the following Mathematica command:
+    // Integrate[(1/2 (1 + Erf[(x - mu)/Sqrt[2*sigma^2]])*Exp[-(x - mu)^2/(2*sigma^2)]*1/Sqrt[2*\[Pi]*sigma^2])/(1 - 1/2 (1 + Erf[(x - mu)/Sqrt[2*sigma^2]])),
+    // {x, -Infinity, x0}, Assumptions -> {x0 \[Element] Reals && mu \[Element] Reals && sigma > 0}]
+    // This gives (for s = x-mu/sqrt(2sigma^2))
+    // -1/2 + erf(s) + log(2/(1+erf(s)))
+
+
+    Real inv_var_scale = 1/(sd*sqrt(Real(2)));
+    Real s0 = (v[0] - mu)*inv_var_scale;
+    Real erfcs0 = erfc(s0, no_promote_policy());
+    // Note that if erfcs0 == 0, then left_tail = inf (numerically), and hence the entire integral is numerically infinite:
+    if (erfcs0 <= 0) {
+        return std::numeric_limits<Real>::infinity();
+    }
+
+    // Note that we're going to add erfcs0/2 when we compute the integral over [x_0, x_1], so drop it here:
+    Real left_tail = -1 + log(Real(2));
+
+
+    // For the right tail, the ecdf is identically 1.
+    // Hence we need the integral:
+    // \int_{v_{n-1}}^{\infty} \frac{(1-F(x))F'(x)}{F(x)} \, \mathrm{d}x
+    // This also has an analytic evaluation! It can be found via the following Mathematica command:
+    // Integrate[(E^(-(z^2/2)) *(1 - 1/2 (1 + Erf[z/Sqrt[2]])))/(Sqrt[2 \[Pi]] (1/2 (1 + Erf[z/Sqrt[2]]))),
+    // {z, zn, \[Infinity]}, Assumptions -> {zn \[Element] Reals && mu \[Element] Reals}]
+    // This gives (for sf = xf-mu/sqrt(2sigma^2))
+    // -1/2 + erf(sf)/2 + 2log(2/(1+erf(sf)))
+
+    Real sf = (v[v.size()-1] - mu)*inv_var_scale;
+    //Real erfcsf = erfc<Real>(sf, no_promote_policy());
+    // This is the actual value of the tail integral. However, the -erfcsf/2 cancels from the integral over [v_{n-2}, v_{n-1}]:
+    //Real right_tail = -erfcsf/2 + log(Real(2)) - log(2-erfcsf);
+
+    // Use erfc(-x) = 2 - erfc(x)
+    Real erfcmsf = erfc<Real>(-sf, no_promote_policy());
+    // Again if this is precisely zero then the integral is numerically infinite:
+    if (erfcmsf == 0) {
+        return std::numeric_limits<Real>::infinity();
+    }
+    Real right_tail = log(2/erfcmsf);
+
+    // Now we need each integral:
+    // \int_{v_i}^{v_{i+1}} \frac{(i+1/n - F(x))^2F'(x)}{F(x)(1-F(x))}  \, \mathrm{d}x
+    // Again we get an analytical evaluation via the following Mathematica command:
+    // Integrate[((E^(-(z^2/2))/Sqrt[2 \[Pi]])*(k1 - F[z])^2)/(F[z]*(1 - F[z])),
+    // {z, z1, z2}, Assumptions -> {z1 \[Element] Reals && z2 \[Element] Reals &&k1 \[Element] Reals}] // FullSimplify
+
+    Real integrals = 0;
+    int64_t N = v.size();
+    for (int64_t i = 0; i < N - 1; ++i) {
+        if (v[i] > v[i+1]) {
+            throw std::domain_error("Input data must be sorted in increasing order v[0] <= v[1] <= . . .  <= v[n-1]");
+        }
+
+        Real k = (i+1)/Real(N);
+        Real s1 = (v[i+1]-mu)*inv_var_scale;
+        Real erfcs1 = erfc<Real>(s1, no_promote_policy());
+        Real term = k*(k*log(erfcs0*(-2 + erfcs1)/(erfcs1*(-2 + erfcs0))) + 2*log(erfcs1/erfcs0));
+
+        integrals += term;
+        s0 = s1;
+        erfcs0 = erfcs1;
+    }
+    integrals -= log(erfcs0);
+    return v.size()*(left_tail + right_tail + integrals);
+}
+
+}}}
+#endif

install/boost_1_75_0/include/boost/math/statistics/bivariate_statistics.hpp

@@ -0,0 +1,96 @@
+//  (C) Copyright Nick Thompson 2018.
+//  Use, modification and distribution are subject to the
+//  Boost Software License, Version 1.0. (See accompanying file
+//  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+#ifndef BOOST_MATH_STATISTICS_BIVARIATE_STATISTICS_HPP
+#define BOOST_MATH_STATISTICS_BIVARIATE_STATISTICS_HPP
+
+#include <iterator>
+#include <tuple>
+#include <boost/assert.hpp>
+
+
+namespace boost{ namespace math{ namespace statistics {
+
+template<class Container>
+auto means_and_covariance(Container const & u, Container const & v)
+{
+    using Real = typename Container::value_type;
+    using std::size;
+    BOOST_ASSERT_MSG(size(u) == size(v), "The size of each vector must be the same to compute covariance.");
+    BOOST_ASSERT_MSG(size(u) > 0, "Computing covariance requires at least one sample.");
+
+    // See Equation III.9 of "Numerically Stable, Single-Pass, Parallel Statistics Algorithms", Bennet et al.
+    Real cov = 0;
+    Real mu_u = u[0];
+    Real mu_v = v[0];
+
+    for(size_t i = 1; i < size(u); ++i)
+    {
+        Real u_tmp = (u[i] - mu_u)/(i+1);
+        Real v_tmp = v[i] - mu_v;
+        cov += i*u_tmp*v_tmp;
+        mu_u = mu_u + u_tmp;
+        mu_v = mu_v + v_tmp/(i+1);
+    }
+
+    return std::make_tuple(mu_u, mu_v, cov/size(u));
+}
+
+template<class Container>
+auto covariance(Container const & u, Container const & v)
+{
+    auto [mu_u, mu_v, cov] = boost::math::statistics::means_and_covariance(u, v);
+    return cov;
+}
+
+template<class Container>
+auto correlation_coefficient(Container const & u, Container const & v)
+{
+    using Real = typename Container::value_type;
+    using std::size;
+    BOOST_ASSERT_MSG(size(u) == size(v), "The size of each vector must be the same to compute covariance.");
+    BOOST_ASSERT_MSG(size(u) > 0, "Computing covariance requires at least two samples.");
+
+    Real cov = 0;
+    Real mu_u = u[0];
+    Real mu_v = v[0];
+    Real Qu = 0;
+    Real Qv = 0;
+
+    for(size_t i = 1; i < size(u); ++i)
+    {
+        Real u_tmp = u[i] - mu_u;
+        Real v_tmp = v[i] - mu_v;
+        Qu = Qu + (i*u_tmp*u_tmp)/(i+1);
+        Qv = Qv + (i*v_tmp*v_tmp)/(i+1);
+        cov += i*u_tmp*v_tmp/(i+1);
+        mu_u = mu_u + u_tmp/(i+1);
+        mu_v = mu_v + v_tmp/(i+1);
+    }
+
+    // If both datasets are constant, then they are perfectly correlated.
+    if (Qu == 0 && Qv == 0)
+    {
+        return Real(1);
+    }
+    // If one dataset is constant and the other isn't, then they have no correlation:
+    if (Qu == 0 || Qv == 0)
+    {
+        return Real(0);
+    }
+
+    // Make sure rho in [-1, 1], even in the presence of numerical noise.
+    Real rho = cov/sqrt(Qu*Qv);
+    if (rho > 1) {
+        rho = 1;
+    }
+    if (rho < -1) {
+        rho = -1;
+    }
+    return rho;
+}
+
+}}}
+#endif

install/boost_1_75_0/include/boost/math/statistics/linear_regression.hpp

@@ -0,0 +1,95 @@
+/*
+ * Copyright Nick Thompson, 2019
+ * Use, modification and distribution are subject to the
+ * Boost Software License, Version 1.0. (See accompanying file
+ * LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+ */
+
+#ifndef BOOST_MATH_STATISTICS_LINEAR_REGRESSION_HPP
+#define BOOST_MATH_STATISTICS_LINEAR_REGRESSION_HPP
+
+#include <cmath>
+#include <algorithm>
+#include <utility>
+#include <boost/math/statistics/univariate_statistics.hpp>
+#include <boost/math/statistics/bivariate_statistics.hpp>
+
+namespace boost::math::statistics {
+
+
+template<class RandomAccessContainer>
+auto simple_ordinary_least_squares(RandomAccessContainer const & x,
+                                   RandomAccessContainer const & y)
+{
+    using Real = typename RandomAccessContainer::value_type;
+    if (x.size() <= 1)
+    {
+        throw std::domain_error("At least 2 samples are required to perform a linear regression.");
+    }
+
+    if (x.size() != y.size())
+    {
+        throw std::domain_error("The same number of samples must be in the independent and dependent variable.");
+    }
+    auto [mu_x, mu_y, cov_xy] = boost::math::statistics::means_and_covariance(x, y);
+
+    auto var_x = boost::math::statistics::variance(x);
+
+    if (var_x <= 0) {
+        throw std::domain_error("Independent variable has no variance; this breaks linear regression.");
+    }
+
+
+    Real c1 = cov_xy/var_x;
+    Real c0 = mu_y - c1*mu_x;
+
+    return std::make_pair(c0, c1);
+}
+
+template<class RandomAccessContainer>
+auto simple_ordinary_least_squares_with_R_squared(RandomAccessContainer const & x,
+                                   RandomAccessContainer const & y)
+{
+    using Real = typename RandomAccessContainer::value_type;
+    if (x.size() <= 1)
+    {
+        throw std::domain_error("At least 2 samples are required to perform a linear regression.");
+    }
+
+    if (x.size() != y.size())
+    {
+        throw std::domain_error("The same number of samples must be in the independent and dependent variable.");
+    }
+    auto [mu_x, mu_y, cov_xy] = boost::math::statistics::means_and_covariance(x, y);
+
+    auto var_x = boost::math::statistics::variance(x);
+
+    if (var_x <= 0) {
+        throw std::domain_error("Independent variable has no variance; this breaks linear regression.");
+    }
+
+
+    Real c1 = cov_xy/var_x;
+    Real c0 = mu_y - c1*mu_x;
+
+    Real squared_residuals = 0;
+    Real squared_mean_deviation = 0;
+    for(decltype(y.size()) i = 0; i < y.size(); ++i) {
+        squared_mean_deviation += (y[i] - mu_y)*(y[i]-mu_y);
+        Real ei = (c0 + c1*x[i]) - y[i];
+        squared_residuals += ei*ei;
+    }
+
+    Real Rsquared;
+    if (squared_mean_deviation == 0) {
+        // Then y = constant, so the linear regression is perfect.
+        Rsquared = 1;
+    } else {
+        Rsquared = 1 - squared_residuals/squared_mean_deviation;
+    }
+
+    return std::make_tuple(c0, c1, Rsquared);
+}
+
+}
+#endif

install/boost_1_75_0/include/boost/math/statistics/ljung_box.hpp

@@ -0,0 +1,70 @@
+//  (C) Copyright Nick Thompson 2019.
+//  Use, modification and distribution are subject to the
+//  Boost Software License, Version 1.0. (See accompanying file
+//  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+#ifndef BOOST_MATH_STATISTICS_LJUNG_BOX_HPP
+#define BOOST_MATH_STATISTICS_LJUNG_BOX_HPP
+
+#include <cmath>
+#include <iterator>
+#include <utility>
+#include <boost/math/distributions/chi_squared.hpp>
+#include <boost/math/statistics/univariate_statistics.hpp>
+
+namespace boost::math::statistics {
+
+template<class RandomAccessIterator>
+auto ljung_box(RandomAccessIterator begin, RandomAccessIterator end, int64_t lags = -1, int64_t fit_dof = 0) {
+    using Real = typename std::iterator_traits<RandomAccessIterator>::value_type;
+    int64_t n = std::distance(begin, end);
+    if (lags >= n) {
+      throw std::domain_error("Number of lags must be < number of elements in array.");
+    }
+
+    if (lags == -1) {
+      // This is the same default as Mathematica; it seems sensible enough . . .
+      lags = static_cast<int64_t>(std::ceil(std::log(Real(n))));
+    }
+
+    if (lags <= 0) {
+      throw std::domain_error("Must have at least one lag.");
+    }
+
+    auto mu = boost::math::statistics::mean(begin, end);
+
+    std::vector<Real> r(lags + 1, Real(0));
+    for (size_t i = 0; i < r.size(); ++i) {
+      for (auto it = begin + i; it != end; ++it) {
+        Real ak = *(it) - mu;
+        Real akml = *(it-i) - mu;
+        r[i] += ak*akml;
+      }
+    }
+
+    Real Q = 0;
+
+    for (size_t k = 1; k < r.size(); ++k) {
+      Q += r[k]*r[k]/(r[0]*r[0]*(n-k));
+    }
+    Q *= n*(n+2);
+
+    typedef boost::math::policies::policy<
+          boost::math::policies::promote_float<false>,
+          boost::math::policies::promote_double<false> >
+          no_promote_policy;
+
+    auto chi = boost::math::chi_squared_distribution<Real, no_promote_policy>(Real(lags - fit_dof));
+
+    Real pvalue = 1 - boost::math::cdf(chi, Q);
+    return std::make_pair(Q, pvalue);
+}
+
+
+template<class RandomAccessContainer>
+auto ljung_box(RandomAccessContainer const & v, int64_t lags = -1, int64_t fit_dof = 0) {
+    return ljung_box(v.begin(), v.end(), lags, fit_dof);
+}
+
+}
+#endif

install/boost_1_75_0/include/boost/math/statistics/runs_test.hpp

@@ -0,0 +1,121 @@
+/*
+ * Copyright Nick Thompson, 2019
+ * Use, modification and distribution are subject to the
+ * Boost Software License, Version 1.0. (See accompanying file
+ * LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+ */
+
+#ifndef BOOST_MATH_STATISTICS_RUNS_TEST_HPP
+#define BOOST_MATH_STATISTICS_RUNS_TEST_HPP
+
+#include <cmath>
+#include <algorithm>
+#include <utility>
+#include <boost/math/statistics/univariate_statistics.hpp>
+#include <boost/math/distributions/normal.hpp>
+
+namespace boost::math::statistics {
+
+template<class RandomAccessContainer>
+auto runs_above_and_below_threshold(RandomAccessContainer const & v,
+                          typename RandomAccessContainer::value_type threshold)
+{
+    using Real = typename RandomAccessContainer::value_type;
+    using std::sqrt;
+    using std::abs;
+    if (v.size() <= 1)
+    {
+        throw std::domain_error("At least 2 samples are required to get number of runs.");
+    }
+    typedef boost::math::policies::policy<
+          boost::math::policies::promote_float<false>,
+          boost::math::policies::promote_double<false> >
+          no_promote_policy;
+
+    decltype(v.size()) nabove = 0;
+    decltype(v.size()) nbelow = 0;
+
+    decltype(v.size()) imin = 0;
+
+    // Take care of the case that v[0] == threshold:
+    while (imin < v.size() && v[imin] == threshold) {
+        ++imin;
+    }
+
+    // Take care of the constant vector case:
+    if (imin == v.size()) {
+        return std::make_pair(std::numeric_limits<Real>::quiet_NaN(), Real(0));
+    }
+
+    bool run_up = (v[imin] > threshold);
+    if (run_up) {
+        ++nabove;
+    } else {
+        ++nbelow;
+    }
+    decltype(v.size()) runs = 1;
+    for (decltype(v.size()) i = imin + 1; i < v.size(); ++i) {
+      if (v[i] == threshold) {
+        // skip values precisely equal to threshold (following R's randtests package)
+        continue;
+      }
+      bool above = (v[i] > threshold);
+      if (above) {
+          ++nabove;
+      } else {
+          ++nbelow;
+      }
+      if (run_up == above) {
+        continue;
+      }
+      else {
+        run_up = above;
+        runs++;
+      }
+    }
+
+    // If you make n an int, the subtraction is gonna be bad in the variance:
+    Real n = nabove + nbelow;
+
+    Real expected_runs = Real(1) + Real(2*nabove*nbelow)/Real(n);
+    Real variance = 2*nabove*nbelow*(2*nabove*nbelow-n)/Real(n*n*(n-1));
+
+    // Bizarre, pathological limits:
+    if (variance == 0)
+    {
+        if (runs == expected_runs)
+        {
+            Real statistic = 0;
+            Real pvalue = 1;
+            return std::make_pair(statistic, pvalue);
+        }
+        else
+        {
+            return std::make_pair(std::numeric_limits<Real>::quiet_NaN(), Real(0));
+        }
+    }
+
+    Real sd = sqrt(variance);
+    Real statistic = (runs - expected_runs)/sd;
+
+    auto normal = boost::math::normal_distribution<Real, no_promote_policy>(0,1);
+    Real pvalue = 2*boost::math::cdf(normal, -abs(statistic));
+    return std::make_pair(statistic, pvalue);
+}
+
+template<class RandomAccessContainer>
+auto runs_above_and_below_median(RandomAccessContainer const & v)
+{
+    using Real = typename RandomAccessContainer::value_type;
+    using std::log;
+    using std::sqrt;
+
+    // We have to memcpy v because the median does a partial sort,
+    // and that would be catastrophic for the runs test.
+    auto w = v;
+    Real median = boost::math::statistics::median(w);
+    return runs_above_and_below_threshold(v, median);
+}
+
+}
+#endif

install/boost_1_75_0/include/boost/math/statistics/signal_statistics.hpp

@@ -0,0 +1,343 @@
+//  (C) Copyright Nick Thompson 2018.
+//  Use, modification and distribution are subject to the
+//  Boost Software License, Version 1.0. (See accompanying file
+//  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+#ifndef BOOST_MATH_TOOLS_SIGNAL_STATISTICS_HPP
+#define BOOST_MATH_TOOLS_SIGNAL_STATISTICS_HPP
+
+#include <algorithm>
+#include <iterator>
+#include <boost/assert.hpp>
+#include <boost/math/tools/complex.hpp>
+#include <boost/math/tools/roots.hpp>
+#include <boost/math/statistics/univariate_statistics.hpp>
+
+
+namespace boost::math::statistics {
+
+template<class ForwardIterator>
+auto absolute_gini_coefficient(ForwardIterator first, ForwardIterator last)
+{
+    using std::abs;
+    using RealOrComplex = typename std::iterator_traits<ForwardIterator>::value_type;
+    BOOST_ASSERT_MSG(first != last && std::next(first) != last, "Computation of the Gini coefficient requires at least two samples.");
+
+    std::sort(first, last,  [](RealOrComplex a, RealOrComplex b) { return abs(b) > abs(a); });
+
+
+    decltype(abs(*first)) i = 1;
+    decltype(abs(*first)) num = 0;
+    decltype(abs(*first)) denom = 0;
+    for (auto it = first; it != last; ++it)
+    {
+        decltype(abs(*first)) tmp = abs(*it);
+        num += tmp*i;
+        denom += tmp;
+        ++i;
+    }
+
+    // If the l1 norm is zero, all elements are zero, so every element is the same.
+    if (denom == 0)
+    {
+        decltype(abs(*first)) zero = 0;
+        return zero;
+    }
+    return ((2*num)/denom - i)/(i-1);
+}
+
+template<class RandomAccessContainer>
+inline auto absolute_gini_coefficient(RandomAccessContainer & v)
+{
+    return boost::math::statistics::absolute_gini_coefficient(v.begin(), v.end());
+}
+
+template<class ForwardIterator>
+auto sample_absolute_gini_coefficient(ForwardIterator first, ForwardIterator last)
+{
+    size_t n = std::distance(first, last);
+    return n*boost::math::statistics::absolute_gini_coefficient(first, last)/(n-1);
+}
+
+template<class RandomAccessContainer>
+inline auto sample_absolute_gini_coefficient(RandomAccessContainer & v)
+{
+    return boost::math::statistics::sample_absolute_gini_coefficient(v.begin(), v.end());
+}
+
+
+// The Hoyer sparsity measure is defined in:
+// https://arxiv.org/pdf/0811.4706.pdf
+template<class ForwardIterator>
+auto hoyer_sparsity(const ForwardIterator first, const ForwardIterator last)
+{
+    using T = typename std::iterator_traits<ForwardIterator>::value_type;
+    using std::abs;
+    using std::sqrt;
+    BOOST_ASSERT_MSG(first != last && std::next(first) != last, "Computation of the Hoyer sparsity requires at least two samples.");
+
+    if constexpr (std::is_unsigned<T>::value)
+    {
+        T l1 = 0;
+        T l2 = 0;
+        size_t n = 0;
+        for (auto it = first; it != last; ++it)
+        {
+            l1 += *it;
+            l2 += (*it)*(*it);
+            n += 1;
+        }
+
+        double rootn = sqrt(n);
+        return (rootn - l1/sqrt(l2) )/ (rootn - 1);
+    }
+    else {
+        decltype(abs(*first)) l1 = 0;
+        decltype(abs(*first)) l2 = 0;
+        // We wouldn't need to count the elements if it was a random access iterator,
+        // but our only constraint is that it's a forward iterator.
+        size_t n = 0;
+        for (auto it = first; it != last; ++it)
+        {
+            decltype(abs(*first)) tmp = abs(*it);
+            l1 += tmp;
+            l2 += tmp*tmp;
+            n += 1;
+        }
+        if constexpr (std::is_integral<T>::value)
+        {
+            double rootn = sqrt(n);
+            return (rootn - l1/sqrt(l2) )/ (rootn - 1);
+        }
+        else
+        {
+            decltype(abs(*first)) rootn = sqrt(static_cast<decltype(abs(*first))>(n));
+            return (rootn - l1/sqrt(l2) )/ (rootn - 1);
+        }
+    }
+}
+
+template<class Container>
+inline auto hoyer_sparsity(Container const & v)
+{
+    return boost::math::statistics::hoyer_sparsity(v.cbegin(), v.cend());
+}
+
+
+template<class Container>
+auto oracle_snr(Container const & signal, Container const & noisy_signal)
+{
+    using Real = typename Container::value_type;
+    BOOST_ASSERT_MSG(signal.size() == noisy_signal.size(),
+                     "Signal and noisy_signal must be have the same number of elements.");
+    if constexpr (std::is_integral<Real>::value)
+    {
+        double numerator = 0;
+        double denominator = 0;
+        for (size_t i = 0; i < signal.size(); ++i)
+        {
+            numerator += signal[i]*signal[i];
+            denominator += (noisy_signal[i] - signal[i])*(noisy_signal[i] - signal[i]);
+        }
+        if (numerator == 0 && denominator == 0)
+        {
+            return std::numeric_limits<double>::quiet_NaN();
+        }
+        if (denominator == 0)
+        {
+            return std::numeric_limits<double>::infinity();
+        }
+        return numerator/denominator;
+    }
+    else if constexpr (boost::math::tools::is_complex_type<Real>::value)
+
+    {
+        using std::norm;
+        typename Real::value_type numerator = 0;
+        typename Real::value_type denominator = 0;
+        for (size_t i = 0; i < signal.size(); ++i)
+        {
+            numerator += norm(signal[i]);
+            denominator += norm(noisy_signal[i] - signal[i]);
+        }
+        if (numerator == 0 && denominator == 0)
+        {
+            return std::numeric_limits<typename Real::value_type>::quiet_NaN();
+        }
+        if (denominator == 0)
+        {
+            return std::numeric_limits<typename Real::value_type>::infinity();
+        }
+
+        return numerator/denominator;
+    }
+    else
+    {
+        Real numerator = 0;
+        Real denominator = 0;
+        for (size_t i = 0; i < signal.size(); ++i)
+        {
+            numerator += signal[i]*signal[i];
+            denominator += (signal[i] - noisy_signal[i])*(signal[i] - noisy_signal[i]);
+        }
+        if (numerator == 0 && denominator == 0)
+        {
+            return std::numeric_limits<Real>::quiet_NaN();
+        }
+        if (denominator == 0)
+        {
+            return std::numeric_limits<Real>::infinity();
+        }
+
+        return numerator/denominator;
+    }
+}
+
+template<class Container>
+auto mean_invariant_oracle_snr(Container const & signal, Container const & noisy_signal)
+{
+    using Real = typename Container::value_type;
+    BOOST_ASSERT_MSG(signal.size() == noisy_signal.size(), "Signal and noisy signal must be have the same number of elements.");
+
+    Real mu = boost::math::statistics::mean(signal);
+    Real numerator = 0;
+    Real denominator = 0;
+    for (size_t i = 0; i < signal.size(); ++i)
+    {
+        Real tmp = signal[i] - mu;
+        numerator += tmp*tmp;
+        denominator += (signal[i] - noisy_signal[i])*(signal[i] - noisy_signal[i]);
+    }
+    if (numerator == 0 && denominator == 0)
+    {
+        return std::numeric_limits<Real>::quiet_NaN();
+    }
+    if (denominator == 0)
+    {
+        return std::numeric_limits<Real>::infinity();
+    }
+
+    return numerator/denominator;
+
+}
+
+template<class Container>
+auto mean_invariant_oracle_snr_db(Container const & signal, Container const & noisy_signal)
+{
+    using std::log10;
+    return 10*log10(boost::math::statistics::mean_invariant_oracle_snr(signal, noisy_signal));
+}
+
+
+// Follows the definition of SNR given in Mallat, A Wavelet Tour of Signal Processing, equation 11.16.
+template<class Container>
+auto oracle_snr_db(Container const & signal, Container const & noisy_signal)
+{
+    using std::log10;
+    return 10*log10(boost::math::statistics::oracle_snr(signal, noisy_signal));
+}
+
+// A good reference on the M2M4 estimator:
+// D. R. Pauluzzi and N. C. Beaulieu, "A comparison of SNR estimation techniques for the AWGN channel," IEEE Trans. Communications, Vol. 48, No. 10, pp. 1681-1691, 2000.
+// A nice python implementation:
+// https://github.com/gnuradio/gnuradio/blob/master/gr-digital/examples/snr_estimators.py
+template<class ForwardIterator>
+auto m2m4_snr_estimator(ForwardIterator first, ForwardIterator last, decltype(*first) estimated_signal_kurtosis=1, decltype(*first) estimated_noise_kurtosis=3)
+{
+    BOOST_ASSERT_MSG(estimated_signal_kurtosis > 0, "The estimated signal kurtosis must be positive");
+    BOOST_ASSERT_MSG(estimated_noise_kurtosis > 0, "The estimated noise kurtosis must be positive.");
+    using Real = typename std::iterator_traits<ForwardIterator>::value_type;
+    using std::sqrt;
+    if constexpr (std::is_floating_point<Real>::value || std::numeric_limits<Real>::max_exponent)
+    {
+        // If we first eliminate N, we obtain the quadratic equation:
+        // (ka+kw-6)S^2 + 2M2(3-kw)S + kw*M2^2 - M4 = 0 =: a*S^2 + bs*N + cs = 0
+        // If we first eliminate S, we obtain the quadratic equation:
+        // (ka+kw-6)N^2 + 2M2(3-ka)N + ka*M2^2 - M4 = 0 =: a*N^2 + bn*N + cn = 0
+        // I believe these equations are totally independent quadratics;
+        // if one has a complex solution it is not necessarily the case that the other must also.
+        // However, I can't prove that, so there is a chance that this does unnecessary work.
+        // Future improvements: There are algorithms which can solve quadratics much more effectively than the naive implementation found here.
+        // See: https://stackoverflow.com/questions/48979861/numerically-stable-method-for-solving-quadratic-equations/50065711#50065711
+        auto [M1, M2, M3, M4] = boost::math::statistics::first_four_moments(first, last);
+        if (M4 == 0)
+        {
+            // The signal is constant. There is no noise:
+            return std::numeric_limits<Real>::infinity();
+        }
+        // Change to notation in Pauluzzi, equation 41:
+        auto kw = estimated_noise_kurtosis;
+        auto ka = estimated_signal_kurtosis;
+        // A common case, since it's the default:
+        Real a = (ka+kw-6);
+        Real bs = 2*M2*(3-kw);
+        Real cs = kw*M2*M2 - M4;
+        Real bn = 2*M2*(3-ka);
+        Real cn = ka*M2*M2 - M4;
+        auto [S0, S1] = boost::math::tools::quadratic_roots(a, bs, cs);
+        if (S1 > 0)
+        {
+            auto N = M2 - S1;
+            if (N > 0)
+            {
+                return S1/N;
+            }
+            if (S0 > 0)
+            {
+                N = M2 - S0;
+                if (N > 0)
+                {
+                    return S0/N;
+                }
+            }
+        }
+        auto [N0, N1] = boost::math::tools::quadratic_roots(a, bn, cn);
+        if (N1 > 0)
+        {
+            auto S = M2 - N1;
+            if (S > 0)
+            {
+                return S/N1;
+            }
+            if (N0 > 0)
+            {
+                S = M2 - N0;
+                if (S > 0)
+                {
+                    return S/N0;
+                }
+            }
+        }
+        // This happens distressingly often. It's a limitation of the method.
+        return std::numeric_limits<Real>::quiet_NaN();
+    }
+    else
+    {
+        BOOST_ASSERT_MSG(false, "The M2M4 estimator has not been implemented for this type.");
+        return std::numeric_limits<Real>::quiet_NaN();
+    }
+}
+
+template<class Container>
+inline auto m2m4_snr_estimator(Container const & noisy_signal,  typename Container::value_type estimated_signal_kurtosis=1, typename Container::value_type estimated_noise_kurtosis=3)
+{
+    return m2m4_snr_estimator(noisy_signal.cbegin(), noisy_signal.cend(), estimated_signal_kurtosis, estimated_noise_kurtosis);
+}
+
+template<class ForwardIterator>
+inline auto m2m4_snr_estimator_db(ForwardIterator first, ForwardIterator last, decltype(*first) estimated_signal_kurtosis=1, decltype(*first) estimated_noise_kurtosis=3)
+{
+    using std::log10;
+    return 10*log10(m2m4_snr_estimator(first, last, estimated_signal_kurtosis, estimated_noise_kurtosis));
+}
+
+
+template<class Container>
+inline auto m2m4_snr_estimator_db(Container const & noisy_signal,  typename Container::value_type estimated_signal_kurtosis=1, typename Container::value_type estimated_noise_kurtosis=3)
+{
+    using std::log10;
+    return 10*log10(m2m4_snr_estimator(noisy_signal, estimated_signal_kurtosis, estimated_noise_kurtosis));
+}
+
+}
+#endif

install/boost_1_75_0/include/boost/math/statistics/t_test.hpp

@@ -0,0 +1,50 @@
+//  (C) Copyright Nick Thompson 2019.
+//  Use, modification and distribution are subject to the
+//  Boost Software License, Version 1.0. (See accompanying file
+//  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+#ifndef BOOST_MATH_STATISTICS_T_TEST_HPP
+#define BOOST_MATH_STATISTICS_T_TEST_HPP
+
+#include <cmath>
+#include <iterator>
+#include <utility>
+#include <boost/math/distributions/students_t.hpp>
+#include <boost/math/statistics/univariate_statistics.hpp>
+
+namespace boost::math::statistics {
+
+template<typename Real>
+std::pair<Real, Real> one_sample_t_test(Real sample_mean, Real sample_variance, Real num_samples, Real assumed_mean) {
+    using std::sqrt;
+    typedef boost::math::policies::policy<
+          boost::math::policies::promote_float<false>,
+          boost::math::policies::promote_double<false> >
+          no_promote_policy;
+
+    Real test_statistic = (sample_mean - assumed_mean)/sqrt(sample_variance/num_samples);
+    auto student = boost::math::students_t_distribution<Real, no_promote_policy>(num_samples - 1);
+    Real pvalue;
+    if (test_statistic > 0) {
+        pvalue = 2*boost::math::cdf<Real>(student, -test_statistic);;
+    }
+    else {
+        pvalue = 2*boost::math::cdf<Real>(student, test_statistic);
+    }
+    return std::make_pair(test_statistic, pvalue);
+}
+
+template<class ForwardIterator>
+auto one_sample_t_test(ForwardIterator begin, ForwardIterator end, typename std::iterator_traits<ForwardIterator>::value_type assumed_mean) {
+    using Real = typename std::iterator_traits<ForwardIterator>::value_type;
+    auto [mu, s_sq] = mean_and_sample_variance(begin, end);
+    return one_sample_t_test(mu, s_sq, Real(std::distance(begin, end)), assumed_mean);
+}
+
+template<class Container>
+auto one_sample_t_test(Container const & v, typename Container::value_type assumed_mean) {
+    return one_sample_t_test(v.begin(), v.end(), assumed_mean);
+}
+
+}
+#endif

install/boost_1_75_0/include/boost/math/statistics/univariate_statistics.hpp

@@ -0,0 +1,575 @@
+//  (C) Copyright Nick Thompson 2018.
+//  Use, modification and distribution are subject to the
+//  Boost Software License, Version 1.0. (See accompanying file
+//  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+#ifndef BOOST_MATH_STATISTICS_UNIVARIATE_STATISTICS_HPP
+#define BOOST_MATH_STATISTICS_UNIVARIATE_STATISTICS_HPP
+
+#include <algorithm>
+#include <iterator>
+#include <tuple>
+#include <cmath>
+#include <vector>
+#include <boost/assert.hpp>
+
+namespace boost::math::statistics {
+
+template<class ForwardIterator>
+auto mean(ForwardIterator first, ForwardIterator last)
+{
+    using Real = typename std::iterator_traits<ForwardIterator>::value_type;
+    BOOST_ASSERT_MSG(first != last, "At least one sample is required to compute the mean.");
+    if constexpr (std::is_integral<Real>::value)
+    {
+        double mu = 0;
+        double i = 1;
+        for(auto it = first; it != last; ++it) {
+            mu = mu + (*it - mu)/i;
+            i += 1;
+        }
+        return mu;
+    }
+    else if constexpr (std::is_same_v<typename std::iterator_traits<ForwardIterator>::iterator_category, std::random_access_iterator_tag>)
+    {
+        size_t elements = std::distance(first, last);
+        Real mu0 = 0;
+        Real mu1 = 0;
+        Real mu2 = 0;
+        Real mu3 = 0;
+        Real i = 1;
+        auto end = last - (elements % 4);
+        for(auto it = first; it != end;  it += 4) {
+            Real inv = Real(1)/i;
+            Real tmp0 = (*it - mu0);
+            Real tmp1 = (*(it+1) - mu1);
+            Real tmp2 = (*(it+2) - mu2);
+            Real tmp3 = (*(it+3) - mu3);
+            // please generate a vectorized fma here
+            mu0 += tmp0*inv;
+            mu1 += tmp1*inv;
+            mu2 += tmp2*inv;
+            mu3 += tmp3*inv;
+            i += 1;
+        }
+        Real num1 = Real(elements  - (elements %4))/Real(4);
+        Real num2 = num1 + Real(elements % 4);
+
+        for (auto it = end; it != last; ++it)
+        {
+            mu3 += (*it-mu3)/i;
+            i += 1;
+        }
+
+        return (num1*(mu0+mu1+mu2) + num2*mu3)/Real(elements);
+    }
+    else
+    {
+        auto it = first;
+        Real mu = *it;
+        Real i = 2;
+        while(++it != last)
+        {
+            mu += (*it - mu)/i;
+            i += 1;
+        }
+        return mu;
+    }
+}
+
+template<class Container>
+inline auto mean(Container const & v)
+{
+    return mean(v.cbegin(), v.cend());
+}
+
+template<class ForwardIterator>
+auto variance(ForwardIterator first, ForwardIterator last)
+{
+    using Real = typename std::iterator_traits<ForwardIterator>::value_type;
+    BOOST_ASSERT_MSG(first != last, "At least one sample is required to compute mean and variance.");
+    // Higham, Accuracy and Stability, equation 1.6a and 1.6b:
+    if constexpr (std::is_integral<Real>::value)
+    {
+        double M = *first;
+        double Q = 0;
+        double k = 2;
+        for (auto it = std::next(first); it != last; ++it)
+        {
+            double tmp = *it - M;
+            Q = Q + ((k-1)*tmp*tmp)/k;
+            M = M + tmp/k;
+            k += 1;
+        }
+        return Q/(k-1);
+    }
+    else
+    {
+        Real M = *first;
+        Real Q = 0;
+        Real k = 2;
+        for (auto it = std::next(first); it != last; ++it)
+        {
+            Real tmp = (*it - M)/k;
+            Q += k*(k-1)*tmp*tmp;
+            M += tmp;
+            k += 1;
+        }
+        return Q/(k-1);
+    }
+}
+
+template<class Container>
+inline auto variance(Container const & v)
+{
+    return variance(v.cbegin(), v.cend());
+}
+
+template<class ForwardIterator>
+auto sample_variance(ForwardIterator first, ForwardIterator last)
+{
+    size_t n = std::distance(first, last);
+    BOOST_ASSERT_MSG(n > 1, "At least two samples are required to compute the sample variance.");
+    return n*variance(first, last)/(n-1);
+}
+
+template<class Container>
+inline auto sample_variance(Container const & v)
+{
+    return sample_variance(v.cbegin(), v.cend());
+}
+
+template<class ForwardIterator>
+auto mean_and_sample_variance(ForwardIterator first, ForwardIterator last)
+{
+    using Real = typename std::iterator_traits<ForwardIterator>::value_type;
+    BOOST_ASSERT_MSG(first != last, "At least one sample is required to compute mean and variance.");
+    // Higham, Accuracy and Stability, equation 1.6a and 1.6b:
+    if constexpr (std::is_integral<Real>::value)
+    {
+        double M = *first;
+        double Q = 0;
+        double k = 2;
+        for (auto it = std::next(first); it != last; ++it)
+        {
+            double tmp = *it - M;
+            Q = Q + ((k-1)*tmp*tmp)/k;
+            M = M + tmp/k;
+            k += 1;
+        }
+        return std::pair<double, double>{M, Q/(k-2)};
+    }
+    else
+    {
+        Real M = *first;
+        Real Q = 0;
+        Real k = 2;
+        for (auto it = std::next(first); it != last; ++it)
+        {
+            Real tmp = (*it - M)/k;
+            Q += k*(k-1)*tmp*tmp;
+            M += tmp;
+            k += 1;
+        }
+        return std::pair<Real, Real>{M, Q/(k-2)};
+    }
+}
+
+template<class Container>
+auto mean_and_sample_variance(Container const & v)
+{
+    return mean_and_sample_variance(v.begin(), v.end());
+}
+
+// Follows equation 1.5 of:
+// https://prod.sandia.gov/techlib-noauth/access-control.cgi/2008/086212.pdf
+template<class ForwardIterator>
+auto skewness(ForwardIterator first, ForwardIterator last)
+{
+    using Real = typename std::iterator_traits<ForwardIterator>::value_type;
+    using std::sqrt;
+    BOOST_ASSERT_MSG(first != last, "At least one sample is required to compute skewness.");
+    if constexpr (std::is_integral<Real>::value)
+    {
+        double M1 = *first;
+        double M2 = 0;
+        double M3 = 0;
+        double n = 2;
+        for (auto it = std::next(first); it != last; ++it)
+        {
+            double delta21 = *it - M1;
+            double tmp = delta21/n;
+            M3 = M3 + tmp*((n-1)*(n-2)*delta21*tmp - 3*M2);
+            M2 = M2 + tmp*(n-1)*delta21;
+            M1 = M1 + tmp;
+            n += 1;
+        }
+
+        double var = M2/(n-1);
+        if (var == 0)
+        {
+            // The limit is technically undefined, but the interpretation here is clear:
+            // A constant dataset has no skewness.
+            return double(0);
+        }
+        double skew = M3/(M2*sqrt(var));
+        return skew;
+    }
+    else
+    {
+        Real M1 = *first;
+        Real M2 = 0;
+        Real M3 = 0;
+        Real n = 2;
+        for (auto it = std::next(first); it != last; ++it)
+        {
+            Real delta21 = *it - M1;
+            Real tmp = delta21/n;
+            M3 += tmp*((n-1)*(n-2)*delta21*tmp - 3*M2);
+            M2 += tmp*(n-1)*delta21;
+            M1 += tmp;
+            n += 1;
+        }
+
+        Real var = M2/(n-1);
+        if (var == 0)
+        {
+            // The limit is technically undefined, but the interpretation here is clear:
+            // A constant dataset has no skewness.
+            return Real(0);
+        }
+        Real skew = M3/(M2*sqrt(var));
+        return skew;
+    }
+}
+
+template<class Container>
+inline auto skewness(Container const & v)
+{
+    return skewness(v.cbegin(), v.cend());
+}
+
+// Follows equation 1.5/1.6 of:
+// https://prod.sandia.gov/techlib-noauth/access-control.cgi/2008/086212.pdf
+template<class ForwardIterator>
+auto first_four_moments(ForwardIterator first, ForwardIterator last)
+{
+    using Real = typename std::iterator_traits<ForwardIterator>::value_type;
+    BOOST_ASSERT_MSG(first != last, "At least one sample is required to compute the first four moments.");
+    if constexpr (std::is_integral<Real>::value)
+    {
+        double M1 = *first;
+        double M2 = 0;
+        double M3 = 0;
+        double M4 = 0;
+        double n = 2;
+        for (auto it = std::next(first); it != last; ++it)
+        {
+            double delta21 = *it - M1;
+            double tmp = delta21/n;
+            M4 = M4 + tmp*(tmp*tmp*delta21*((n-1)*(n*n-3*n+3)) + 6*tmp*M2 - 4*M3);
+            M3 = M3 + tmp*((n-1)*(n-2)*delta21*tmp - 3*M2);
+            M2 = M2 + tmp*(n-1)*delta21;
+            M1 = M1 + tmp;
+            n += 1;
+        }
+
+        return std::make_tuple(M1, M2/(n-1), M3/(n-1), M4/(n-1));
+    }
+    else
+    {
+        Real M1 = *first;
+        Real M2 = 0;
+        Real M3 = 0;
+        Real M4 = 0;
+        Real n = 2;
+        for (auto it = std::next(first); it != last; ++it)
+        {
+            Real delta21 = *it - M1;
+            Real tmp = delta21/n;
+            M4 = M4 + tmp*(tmp*tmp*delta21*((n-1)*(n*n-3*n+3)) + 6*tmp*M2 - 4*M3);
+            M3 = M3 + tmp*((n-1)*(n-2)*delta21*tmp - 3*M2);
+            M2 = M2 + tmp*(n-1)*delta21;
+            M1 = M1 + tmp;
+            n += 1;
+        }
+
+        return std::make_tuple(M1, M2/(n-1), M3/(n-1), M4/(n-1));
+    }
+}
+
+template<class Container>
+inline auto first_four_moments(Container const & v)
+{
+    return first_four_moments(v.cbegin(), v.cend());
+}
+
+
+// Follows equation 1.6 of:
+// https://prod.sandia.gov/techlib-noauth/access-control.cgi/2008/086212.pdf
+template<class ForwardIterator>
+auto kurtosis(ForwardIterator first, ForwardIterator last)
+{
+    auto [M1, M2, M3, M4] = first_four_moments(first, last);
+    if (M2 == 0)
+    {
+        return M2;
+    }
+    return M4/(M2*M2);
+}
+
+template<class Container>
+inline auto kurtosis(Container const & v)
+{
+    return kurtosis(v.cbegin(), v.cend());
+}
+
+template<class ForwardIterator>
+auto excess_kurtosis(ForwardIterator first, ForwardIterator last)
+{
+    return kurtosis(first, last) - 3;
+}
+
+template<class Container>
+inline auto excess_kurtosis(Container const & v)
+{
+    return excess_kurtosis(v.cbegin(), v.cend());
+}
+
+
+template<class RandomAccessIterator>
+auto median(RandomAccessIterator first, RandomAccessIterator last)
+{
+    size_t num_elems = std::distance(first, last);
+    BOOST_ASSERT_MSG(num_elems > 0, "The median of a zero length vector is undefined.");
+    if (num_elems & 1)
+    {
+        auto middle = first + (num_elems - 1)/2;
+        std::nth_element(first, middle, last);
+        return *middle;
+    }
+    else
+    {
+        auto middle = first + num_elems/2 - 1;
+        std::nth_element(first, middle, last);
+        std::nth_element(middle, middle+1, last);
+        return (*middle + *(middle+1))/2;
+    }
+}
+
+
+template<class RandomAccessContainer>
+inline auto median(RandomAccessContainer & v)
+{
+    return median(v.begin(), v.end());
+}
+
+template<class RandomAccessIterator>
+auto gini_coefficient(RandomAccessIterator first, RandomAccessIterator last)
+{
+    using Real = typename std::iterator_traits<RandomAccessIterator>::value_type;
+    BOOST_ASSERT_MSG(first != last && std::next(first) != last, "Computation of the Gini coefficient requires at least two samples.");
+
+    std::sort(first, last);
+    if constexpr (std::is_integral<Real>::value)
+    {
+        double i = 1;
+        double num = 0;
+        double denom = 0;
+        for (auto it = first; it != last; ++it)
+        {
+            num += *it*i;
+            denom += *it;
+            ++i;
+        }
+
+        // If the l1 norm is zero, all elements are zero, so every element is the same.
+        if (denom == 0)
+        {
+            return double(0);
+        }
+
+        return ((2*num)/denom - i)/(i-1);
+    }
+    else
+    {
+        Real i = 1;
+        Real num = 0;
+        Real denom = 0;
+        for (auto it = first; it != last; ++it)
+        {
+            num += *it*i;
+            denom += *it;
+            ++i;
+        }
+
+        // If the l1 norm is zero, all elements are zero, so every element is the same.
+        if (denom == 0)
+        {
+            return Real(0);
+        }
+
+        return ((2*num)/denom - i)/(i-1);
+    }
+}
+
+template<class RandomAccessContainer>
+inline auto gini_coefficient(RandomAccessContainer & v)
+{
+    return gini_coefficient(v.begin(), v.end());
+}
+
+template<class RandomAccessIterator>
+inline auto sample_gini_coefficient(RandomAccessIterator first, RandomAccessIterator last)
+{
+    size_t n = std::distance(first, last);
+    return n*gini_coefficient(first, last)/(n-1);
+}
+
+template<class RandomAccessContainer>
+inline auto sample_gini_coefficient(RandomAccessContainer & v)
+{
+    return sample_gini_coefficient(v.begin(), v.end());
+}
+
+template<class RandomAccessIterator>
+auto median_absolute_deviation(RandomAccessIterator first, RandomAccessIterator last, typename std::iterator_traits<RandomAccessIterator>::value_type center=std::numeric_limits<typename std::iterator_traits<RandomAccessIterator>::value_type>::quiet_NaN())
+{
+    using std::abs;
+    using Real = typename std::iterator_traits<RandomAccessIterator>::value_type;
+    using std::isnan;
+    if (isnan(center))
+    {
+        center = boost::math::statistics::median(first, last);
+    }
+    size_t num_elems = std::distance(first, last);
+    BOOST_ASSERT_MSG(num_elems > 0, "The median of a zero-length vector is undefined.");
+    auto comparator = [&center](Real a, Real b) { return abs(a-center) < abs(b-center);};
+    if (num_elems & 1)
+    {
+        auto middle = first + (num_elems - 1)/2;
+        std::nth_element(first, middle, last, comparator);
+        return abs(*middle);
+    }
+    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<Real>(2));
+    }
+}
+
+template<class RandomAccessContainer>
+inline auto median_absolute_deviation(RandomAccessContainer & v, typename RandomAccessContainer::value_type center=std::numeric_limits<typename RandomAccessContainer::value_type>::quiet_NaN())
+{
+    return median_absolute_deviation(v.begin(), v.end(), center);
+}
+
+template<class ForwardIterator>
+auto interquartile_range(ForwardIterator first, ForwardIterator last)
+{
+    using Real = typename std::iterator_traits<ForwardIterator>::value_type;
+    static_assert(!std::is_integral<Real>::value, "Integer values have not yet been implemented.");
+    auto m = std::distance(first,last);
+    BOOST_ASSERT_MSG(m >= 3, "At least 3 samples are required to compute the interquartile range.");
+    auto k = m/4;
+    auto j = m - (4*k);
+    // m = 4k+j.
+    // If j = 0 or j = 1, then there are an even number of samples below the median, and an even number above the median.
+    //    Then we must average adjacent elements to get the quartiles.
+    // If j = 2 or j = 3, there are an odd number of samples above and below the median, these elements may be directly extracted to get the quartiles.
+
+    if (j==2 || j==3)
+    {
+        auto q1 = first + k;
+        auto q3 = first + 3*k + j - 1;
+        std::nth_element(first, q1, last);
+        Real Q1 = *q1;
+        std::nth_element(q1, q3, last);
+        Real Q3 = *q3;
+        return Q3 - Q1;
+    } else {
+        // j == 0 or j==1:
+        auto q1 = first + k - 1;
+        auto q3 = first + 3*k - 1 + j;
+        std::nth_element(first, q1, last);
+        Real a = *q1;
+        std::nth_element(q1, q1 + 1, last);
+        Real b = *(q1 + 1);
+        Real Q1 = (a+b)/2;
+        std::nth_element(q1, q3, last);
+        a = *q3;
+        std::nth_element(q3, q3 + 1, last);
+        b = *(q3 + 1);
+        Real Q3 = (a+b)/2;
+        return Q3 - Q1;
+    }
+}
+
+template<class RandomAccessContainer>
+inline auto interquartile_range(RandomAccessContainer & v)
+{
+    return interquartile_range(v.begin(), v.end());
+}
+
+template<class ForwardIterator, class OutputIterator>
+auto sorted_mode(ForwardIterator first, ForwardIterator last, OutputIterator output) -> decltype(output)
+{
+    using Z = typename std::iterator_traits<ForwardIterator>::value_type;
+    static_assert(std::is_integral<Z>::value, "Floating point values have not yet been implemented.");
+    using Size = typename std::iterator_traits<ForwardIterator>::difference_type;
+
+    std::vector<Z> modes {};
+    modes.reserve(16);
+    Size max_counter {0};
+
+    while(first != last)
+    {
+        Size current_count {0};
+        auto end_it {first};
+        while(end_it != last && *end_it == *first)
+        {
+            ++current_count;
+            ++end_it;
+        }
+
+        if(current_count > max_counter)
+        {
+            modes.resize(1);
+            modes[0] = *first;
+            max_counter = current_count;
+        }
+
+        else if(current_count == max_counter)
+        {
+            modes.emplace_back(*first);
+        }
+
+        first = end_it;
+    }
+
+    return std::move(modes.begin(), modes.end(), output);
+}
+
+template<class Container, class OutputIterator>
+inline auto sorted_mode(Container & v, OutputIterator output) -> decltype(output)
+{
+    return sorted_mode(v.begin(), v.end(), output);
+}
+
+template<class RandomAccessIterator, class OutputIterator>
+auto mode(RandomAccessIterator first, RandomAccessIterator last, OutputIterator output) -> decltype(output)
+{
+    std::sort(first, last);
+    return sorted_mode(first, last, output);
+}
+
+template<class RandomAccessContainer, class OutputIterator>
+inline auto mode(RandomAccessContainer & v, OutputIterator output) -> decltype(output)
+{
+    return mode(v.begin(), v.end(), output);
+}
+
+}
+#endif