Standard deviation in one pass in C++

Question

(See the next iteration.)

I have this template function for computing standard deviation in one pass over the data. The underlying math is here.

coderodde_sd.h:

#ifndef CODERODDE_SD_H
#define CODERODDE_SD_H

#include <cmath>
#include <iterator>
#include <sstream>
#include <stdexcept>

namespace net {

    namespace coderodde {

        namespace stat {

            template<class Iter>
            double sd(Iter begin, Iter end)
            {
                typedef typename std::iterator_traits<Iter>::difference_type 
                                                             difference_type;

                difference_type distance = std::distance(begin, end);

                if (distance < 2)
                {
                    std::stringstream ss;
                    ss << "The standard deviation cannot be computed for "
                          "less than two elements. The input sequence has "
                       << distance 
                       << " elements.";

                    throw std::runtime_error(ss.str());
                }

                typedef typename std::iterator_traits<Iter>::value_type 
                                                             value_type;

                double x = 0.0;
                double x_squared = 0.0;

                for (Iter it = begin; it != end; ++it) 
                {
                    x += *it;
                    x_squared += (*it) * (*it);
                }

                return std::sqrt((x_squared - (x * x) / distance) / 
                                 (distance - 1)
                                );
            }

        } /* net::coderodde::stat */

    } /* net::coderodde */

} /* net */

#endif  /* CODERODDE_SD_H */

main.cpp:

#include <iostream>
#include <list>
#include "coderodde_sd.h"

using net::coderodde::stat::sd;
using std::cout;
using std::list;

int main(int argc, char** argv) {
    double bad_array[]{1.0};

    try 
    {
        sd(bad_array, bad_array);
    }
    catch (std::runtime_error& error)
    {
        cout << "ERROR: " << error.what() << "\n";
    }

    try 
    {
        sd(bad_array, bad_array + 1);
    }
    catch (std::runtime_error& error)
    {
        cout << "ERROR: " << error.what() << "\n";
    }

    std::list<int> my_list = { 1, 5, 2, 4, 3 };

    cout << "Standard deviation: " 
         << sd(my_list.begin(), my_list.end())
         << "\n";

    return 0;
}

I tried hard to write idiomatic C++; please tell me where I failed to reach that ideal + any other critique.

Answer 1

In your example, you have two passes over the list - std::distance works in linear time unless it operates on random access iterators.

A naive way to fix this is to calculate the number of elements inside your for loop and check it afterwards. But this way you'll pessimize performance for random access iterators.

To avoid this pessimization you can examine the iterator_tag of you Iter at compile time and decide whether you're going to call std::distance or count your elements in a loop.

The simplest way to do this is in C++14 is by tag dispatching. Here is an example:

template<typename RandomAccessIter>
double sd(RandomAccessIter begin, RandomAccessIter end, std::random_access_iterator_tag)
{
    // calling std::distance
}

template<typename ForwardIter>
double sd(ForwardIter begin, ForwardIter end, std::forward_iterator_tag)
{
    // counting elements in a loop
}

template<typename Iter>
double sd(Iter begin, Iter end)
{
    // dispatching
    return sd(begin, end, typename std::iterator_traits<Iter>::iterator_category());
}

This code will corretly handle bidirectional iterators as well since std::bidirectional_iterator_tag inherits from std::forward_iterator_tag. And code duplication can be avoided by incaplusating common logic in separate functions.

Viola! Now you have only one pass.

Answer 2

I'd rename your template parameter name - instead of Iter use a name that indicates what type of iterator is expected. In this case InputIterator should be fine.

Assuming you're using at least C++11, it would be more idiomatic to use

using difference_type = std::iterator_traits<Iter>::difference_type;
using value_type = std::iterator_traits<Iter>::value_type;

You may also consider using some more functional features

struct StdDev {
    double value = 0.0;
    double squared = 0.0;
};

template <class InputIterator>
StdDev adder(const StdDev& left, const InputIterator& right) {
    double rightVal = *right;
    double rightSquared = rightVal * rightVal;

    return StdDev{left.value + rightVal, left.squared + rightSquared };
}

...

auto val = std::accumulate(begin, end, StdDev(), adder<Iter>);

return std::sqrt((val.squared - (val.value * val.value) / distance) / (distance - 1));