suprising behavior of boost::math::statistics::mean()

[lewandowski-a-j-gmail-com]

I do know this is not a bug, but please consider rejecting, at compile time, the input iterators for boost::math::statistics::mean()
It cost me some time to figure out why it doesn't give excactly the same results
as boost::math::statistics::first_four_moments() (I was pasing std::istream_iterator).
The point is, for the data I was initially playing with, the results were very close, I was even thinking that mean() uses some improved algorithm to reduce rounding errors, but I found an example (below) that shows that mean() yields random results for input iterators
(variance, skewness, kurtosis and first_four_moments work just fine in that case -- and this is even more confusing)

BTW, this is a great library, and I'm going to recommend it to my students, that is why I'm filling 3rd issue today ;)

I do realize that docs clearly states ForwardIterator, but the docs also says
experienced programmers will recognize that even trivial algorithms are easy to screw up
and I would add: and so are the requirements not enforced by the compiler.

Here is the code:

#include <boost/math/statistics/univariate_statistics.hpp>

#include <sstream>
#include <algorithm>
#include <cassert>
#include <iostream>

int main()
{
        auto ifs =std::istringstream("1 2 3 4 5 6 7 8 9");
        
        auto begin = std::istream_iterator<double>(ifs);
        const auto end   = decltype(begin)();
        
        namespace bms = boost::math::statistics;
        auto mean = bms::mean(begin,end);
        
        
        ifs.clear();
        ifs.seekg(0,std::istringstream::beg);
        
        begin = std::istream_iterator<double>(ifs);
        
        auto [mean2,x,y,z] = bms::first_four_moments(begin,end);
        
        std::cout << "mean:  "<<mean<<std::endl; //prints 1 (wrong)
        
        std::cout << "mean2: "<<mean2<<std::endl;//prints 5 (correct)       

        ifs.clear();
        ifs.seekg(0,std::istringstream::beg);
        
        mean = bms::mean(begin,end);
        std::cout << "mean:  "<<mean<<std::endl; //prints 1 (again wrong)
        
}