Mean and Variance for Math Library

Question

Taking forward the code written for the math library I'm developing. I wrote templatized functions for mean and variance.

1. std::accumulate() (or summation in general) is likely to cause overflow. Should I go for an incremental approach?

2. Another point of concern is for templatized argument, I am forced to use the same type as the return type, which means the mean of integers shall also result in an integer, clearly not a good choice. At the same time if I enforce double as the return type, then extending into multi-dimensional types like vectors or complex numbers will be problematic.

3. Is the comparison for double appropriate? Is there a way to write a generalized comparison function for all float types?

CODE

#include <vector>
#include <numeric>
#include <string>
#include <functional>

namespace Statistics
{
    template <typename T>
    T average(std::vector<T> distributionVector)
    {
        if (distributionVector.size() == 0)
        {
            throw std::invalid_argument("Statistics::average - The distribution provided is empty");
        }
        return std::accumulate(distributionVector.begin(), distributionVector.end(), T()) 
            / (distributionVector.size());
    }

    template <typename T>
    T variance(std::vector<T> distributionVector)
    {
        if (distributionVector.size() == 0)
        {
            throw std::invalid_argument("Statistics::expectation - The distribution provided is empty");
        }
        T sumOfSquare = std::accumulate(distributionVector.begin(), distributionVector.end(), T(), [](T a,T b) { return a + b*b; });
        T meanOfSquare = sumOfSquare / distributionVector.size();
        T squareOfMean = (average(distributionVector)) * (average(distributionVector));
        return (meanOfSquare - squareOfMean);
    }
}

Test code

#include "pch.h"
#include <vector>
#include "../MathLibrary/Statistics.h"

void compareDoubles(double a, double b)
{
    const double THRESHOLD = 0.01;
    ASSERT_TRUE(abs(a - b) < THRESHOLD);
}
TEST(Statistics_mean, small_distributions)
{
    std::vector<int> testVector = { -2,-1,0,1,2 };
    EXPECT_EQ(Statistics::average(testVector), 0);
    std::vector<double> testVectorDouble = {5,5,6,6};
    compareDoubles(Statistics::average(testVectorDouble), 5.5);
}

TEST(Statistics_mean, empty_distribution)
{
    std::vector<int> testVector;
    EXPECT_THROW(Statistics::average(testVector), std::invalid_argument);
}

TEST(Statistics_variance, small_distribution)
{
    std::vector<double> testVector = { 0,0 };
    compareDoubles(Statistics::variance(testVector), 0);
    std::vector<double> testVector2 = {1,2,3,4};
    compareDoubles(Statistics::variance(testVector2), 1.25);
    std::vector<double> testVectorRandom = { 1,2,3,4,6,8,9,34,45,78,89 };
    compareDoubles(Statistics::variance(testVectorRandom), 938.2314);
}

Answer 1

The calculation of variance as Σ(x²)/n - (Σx/n)² doesn't work well for small variances, as it can be a difference of two large numbers, resulting in very poor relative accuracy.

You'll get more accuracy using the basic definition Σ(x-̅x)²/n. Here's a version I wrote recently for a different project:

#include <cmath>
#include <numeric>
#include <ranges>
#include <utility>

template<typename Container>
auto mean_and_stddev(const Container& values)
{
    auto len = static_cast<double>(values.size());
    auto mean = std::accumulate(values.begin(), values.end(), 0.0) / len;
    auto sdv =                   // squared differences view
        values
        | std::views::transform([mean](auto d){ auto x = d - mean; return x*x; })
        | std::views::common;
    auto stddev = std::sqrt(std::accumulate(sdv.begin(), sdv.end(), 0.0) / len);
    return std::pair{ mean, stddev };
}

That's more numerically stable.

---

Other issues:

• We're passing potentially large containers by value. There's no need for that - pass as const ref instead.

• We only support std::vector as container, so we can't pass a std::array for example, or a sub-range of a vector. Consider passing start and end iterators instead.

• Throwing exceptions might not be the best interface for the mean and variance of an empty collection - consider returning a NaN value instead.

• Still missing necessary headers (<cmath> for std::abs() (which was misspelt); <stdexcept> for the exception types, <gtest/gtest.h> for the test library).

• Using integer types as return for integer inputs - consider std::common_type_t<T, double> for the computation and return value.

• Users may want both sample and population variances; it's not much work to provide both.

• compareDoubles() loses a lot of debugging information from failures - we don't get the source location, or the actual and expected values. Consider using EXPECT_DOUBLE_EQ(), or examining its implementation if you really need to provide your own.

• It's hard to see how we obtained the expected variance for testVectorRandom. Prefer to write tests where the results can be calculated with mental arithmetic.

  Example, from my implementation above:

    #include <gtest/gtest.h>
    #include <array>
    TEST(mean_and_stddev, int_values)
    {
        auto [mean, stddev] = mean_and_stddev(std::array{50, 150});
        EXPECT_DOUBLE_EQ(mean, 100.);
        EXPECT_DOUBLE_EQ(stddev, 50.);
    }
  

Answer 2

Review

Welcome to Code Review. There are some suggestions as follows.

Answer your questions

std::accumulate()(or summation in general) is likely to cause overflow. Should I go for an incremental approach?

Why you think that std::accumulate is likely to cause overflow? In your test cases it works well. Whether overflow occurs or not depends on the type you choose. Choose type smart and no overflow happen.

I am forced to use the same type as the return type, which means the mean of integers shall also result in an integer, clearly not a good choice.

You noticed that the result Statistics::average(testVector2) in the following code is 3 not 3.3333.

std::vector<int> testVector2 = { 0, 3, 7 };
std::cout << Statistics::average(testVector2) << "\n";

To solve this issue, one of solutions is to specify the return type. We can use double (as default return type) here because the average of numbers is always a floating number.

template <typename OutputT = double, typename T>
OutputT average(std::vector<T> distributionVector)
{
    if (distributionVector.size() == 0)
    {
        throw std::invalid_argument("Statistics::average - The distribution provided is empty");
    }
    return std::accumulate(distributionVector.begin(), distributionVector.end(), OutputT())
        / (distributionVector.size());
}

Is the comparison for double appropriate? Is there a way to write a generalized comparison function for all float types?

Maybe you can check the post on SO.

Input container for average and variance template function

You specify the input type in std::vector. How about std::array or std::list?

Actually, you don't need to specify std::vector<T> and can use OutputT average(T distributionVector) only. Moreover, with C++20 you can do this:

template <typename OutputT = double, std::ranges::input_range T>
constexpr OutputT average(const T& distributionVector)

Answer 3

I'll just address the concerns about the summation.

std::accumulate() (or summation in general) is likely to cause overflow. Should I go for an incremental approach?

You can query the maximum value for a double using std::numeric_limits<double>::max(). This will most likely tell you that it's \$1.79769 \cdot 10^{308}\$. Since you are probably not summing more than \$2^{64} \approx 1.8 \cdot 10^{19}\$ values, and you are squaring numbers which means doubling their exponent, if the maximum value of the elements in the input is no more than roughly \$10^{\frac{308-19}{2}} \approx 10^{144}\$, you'll be fine.

However, what is more of a concern is the accuracy of the results. If you add a very small number to a very large one, the limited precision of a double will mean that the contribution of the small number gets lost. If you have a few big values and lots of small ones, this can be a problem. The solution to this problem is to use pairwise summation. Also heed Toby Speight's advice and calculate the variance using the formula \$\frac{\sum(x-\bar x)^2}{n}\$.

Answer 4

overflow checking

Generally, all the code in the standard library and elsewhere is written without inserting any overflow checking. Unless you have a special need, don't worry about it in your functions.

Now that it's a template, the caller could specify an extended-precision integer class, or a "safe" integer class, to get such checking. The checking is part of the type used, and need not be explicitly addressed in your code. It will be built into the operator+ etc. for that type. For example, see the videos from the 2021 C++now conference — I recall a presentation on Simplest Safe Integers, and there have been others in the past.

parameter, return, and other types

Another point of concern is for templatized argument, I am forced to use the same type as the return type, which means the mean of integers shall also result in an integer, clearly not a good choice. At the same time if I enforce double as the return type, then extending into multi-dimensional types like vectors or complex numbers will be problematic.

No, you are not forced to use the same type for the argument as the return.

You can take the input parameter as another template argument. Using classic syntax:

template <typename R, typename P>
R factorial (P x)
{ ⋯ }

With C++20, you can use Concepts to specify that the template parameters must be integral types.

Note that you put R first, since you will deduce P from the arguments but must give R. Example use:

const auto y = factorial<bignum_t>(432);

Likewise, you can take additional template arguments, perhaps with defaults, to use for internal computations if that becomes necessary.

It's probably much more efficient to use built-in integers for inputs and where you can inside the body of the function, and the extended-precision (or "safe") class for the accumulation of the result.

comparison for floating-point

That should be built-into the unit testing framework. (example)

the new functions

template <typename T>
T average(std::vector<T> distributionVector)

First of all, why is it limited to/specific to vector?
Second, why are you passing it by value? Do you understand that this copies the entire vector?

Classically, such functions should take a pair of iterators to specify the input. As of C++20, you could use a Range. By taking any parameter that satisfies the Range concept, that includes std::vector and any other sequential collection as well.

If you change the template declaration to accept any container, that doesn't affect the code! Except... it would work for std::deque, boost::small_vector, etc. but would have trouble when you pass it a raw array, or certain types of range views. For this reason, use the non-member begin, end, size, etc. to better abstract the collection/view.

Classically, this would be done with the "std:: two-step". With C++20, just use the newer std::ranges::begin etc.

if (distributionVector.size() == 0)

Use empty instead of comparing the size with 0. And again, use the non-member.

When testing, include a plain C array as one of the tests; e.g.

constexpr int test_val_1[] = { 1,2,3,4,6,8,9,34,45,78,89 };
   ⋮
const auto result = variance<double>(test_val_1);

Notice that in this example I also specified that I want a floating-point calculation of the variance, even though the input is integers.

the API

In many calculators, the stats are computed incrementally. You submit each value (or value pair) as you compute them or enter them, often with a key labeled "Σ+". You can model this with a class that keeps intermediate accumulators but does not need the entire list of input at once. It can be used incrementally.

The class would use template arguments for the type to use for the accumulators, which IIRC are things like the count of items (n), the sum, the sum of the squares, and more for two-variable input.

The class would have functions to submit a value, or submit a range of values at once. These member functions can themselves be templates, to be flexible as to the parameter type.

It would have functions like average and variance that operate on the stored data. This way you can call for multiple stats without having to re-compute things or make multiple passes over the same data.

Answer 5

A common way to avoid overflow (when using integers) or accuracy loss (when using floating point type is to use a guess. If you do not know the exact value of the mean but know it to an order of magnitude, you can use the formula: \$\overline {x - guess} = \bar x - guess\$. If you choose your guess wisely, the accuracy gain can be important.