105

Is there a way to calculate mean and standard deviation for a vector containing samples using Boost?

Or do I have to create an accumulator and feed the vector into it?

Improve this question
0

10 Answers 10

Reset to default
251

I don't know if Boost has more specific functions, but you can do it with the standard library.

Given std::vector<double> v, this is the naive way:

#include <numeric> double sum = std::accumulate(v.begin(), v.end(), 0.0); double mean = sum / v.size(); double sq_sum = std::inner_product(v.begin(), v.end(), v.begin(), 0.0); double stdev = std::sqrt(sq_sum / v.size() - mean * mean);

This is susceptible to overflow or underflow for huge or tiny values. A slightly better way to calculate the standard deviation is:

double sum = std::accumulate(v.begin(), v.end(), 0.0); double mean = sum / v.size(); std::vector<double> diff(v.size()); std::transform(v.begin(), v.end(), diff.begin(), std::bind2nd(std::minus<double>(), mean)); double sq_sum = std::inner_product(diff.begin(), diff.end(), diff.begin(), 0.0); double stdev = std::sqrt(sq_sum / v.size());

UPDATE for C++11:

The call to std::transform can be written using a lambda function instead of std::minus and std::bind2nd(now deprecated):

std::transform(v.begin(), v.end(), diff.begin(), [mean](double x) { return x - mean; });
Improve this answer
Sign up to request clarification or add additional context in comments.

13 Comments

The first set of equations does not work. I put int 10 & 2, and got an output of 4. At a glance I think it's b/c it assumes that (a-b)^2 = a^2-b^2
@CharlesL.: It should work, and 4 is the correct answer.
@StudentT: No, but you can substitute (v.size() - 1) for v.size() in the last line above: std::sqrt(sq_sum / (v.size() - 1)). (For the first method, it's a little complicated: std::sqrt(sq_sum / (v.size() - 1) - mean * mean * v.size() / (v.size() - 1)).
Using std::inner_product for sum of squares is very neat.
i can confirm firsthand that the first implementation does overflow/underflow for tiny numbers. i had to change to the second implementation and then i did not get NAN value for the standard deviation. The two extra lines of code are worth it to avoid overflow/underflow!
|
85

If performance is important to you, and your compiler supports lambdas, the stdev calculation can be made faster and simpler: In tests with VS 2012 I've found that the following code is over 10 X quicker than the Boost code given in the chosen answer; it's also 5 X quicker than the safer version of the answer using standard libraries given by musiphil.

Note I'm using sample standard deviation, so the below code gives slightly different results (Why there is a Minus One in Standard Deviations)

double sum = std::accumulate(std::begin(v), std::end(v), 0.0); double m = sum / v.size(); double accum = 0.0; std::for_each (std::begin(v), std::end(v), [&](const double d) { accum += (d - m) * (d - m); }); double stdev = sqrt(accum / (v.size()-1));
Improve this answer

8 Comments

Thanks for sharing this answer even a year later. Now I come another year later and made this one generic for both the value type and the container type. See here (Note: I guess that my range-based for loop is as fast as your lambda code.)
what is the difference between using std::end(v) instead of v.end()?
The std::end() function was added by C++11 standard for cases when there is nothing like v.end(). The std::end can be overloaded for the less standard container -- see en.cppreference.com/w/cpp/iterator/end
Can you explain why is this faster?
Well for one thing, the "safe" answer (which is like my answer) makes 3 passes through the array: Once for the sum, once for the diff-mean, and once for the squaring. In my code there's only 2 passes -- It's conflating the second two passes into one. And (when I last looked, quite a while ago now!) the inner_product calls were not optimized away. In addition the "safe" code copies v into an entirely new array of diffs, which adds more delay. In my opinion my code is more readable too - and is easily ported to JavaScript and other languages :)
|
58

Using accumulators is the way to compute means and standard deviations in Boost.

accumulator_set<double, stats<tag::variance> > acc; for_each(a_vec.begin(), a_vec.end(), bind<void>(ref(acc), _1)); cout << mean(acc) << endl; cout << sqrt(variance(acc)) << endl;

 

Improve this answer

2 Comments

Note, that tag::variance calculates variance by an approximate formula. tag::variance(lazy) calculates by an exact formula, specifically: second moment - squared mean which will produce incorrect result if variance is very small because of rounding errors. It can actually produce negative variance.
Use the recursive (online) algorithm if you know you are going to have a lots of numbers. This will take care of both under and overflow problems.
8

Improving on the answer by musiphil, you can write a standard deviation function without the temporary vector diff, just using a single inner_product call with the C++11 lambda capabilities:

double stddev(std::vector<double> const & func) { double mean = std::accumulate(func.begin(), func.end(), 0.0) / func.size(); double sq_sum = std::inner_product(func.begin(), func.end(), func.begin(), 0.0, [](double const & x, double const & y) { return x + y; }, [mean](double const & x, double const & y) { return (x - mean)*(y - mean); }); return std::sqrt(sq_sum / func.size() - 1); }

I suspect doing the subtraction multiple times is cheaper than using up additional intermediate storage, and I think it is more readable, but I haven't tested the performance yet.

As for explanation why to use N-1 (as in func.size() - 1), see these questions - note how the question states we have a "vector containing samples".

Improve this answer

5 Comments

I think this is computing the variance, not the standard deviation.
The std deviation is calculed dividing by N and not by N-1. Why do you divide the sq_sum by func.size()-1?
Why N-1? You divide by N when computing a population std dev, and divide by N-1 when computing a sample std dev. stats.stackexchange.com/q/3931 stats.stackexchange.com/q/485326
@JonathonReinhart you mean to say initial computation using N-1 was correct, right?
I'm saying that they are both correct, depending on what you are calculating.
7

It seems the following elegant recursive solution has not been mentioned, although it has been around for a long time. Referring to Knuth's Art of Computer Programming, 

mean_1 = x_1, variance_1 = 0; //initial conditions; edge case; //for k >= 2, mean_k = mean_k-1 + (x_k - mean_k-1) / k; variance_k = variance_k-1 + (x_k - mean_k-1) * (x_k - mean_k);

then for a list of n>=2 values, the estimate of the standard deviation is: 

stddev = std::sqrt(variance_n / (n-1)).

Hope this helps!

Improve this answer

2 Comments

This is pretty cool. I implemented it with an index loop ( pastebin.com/aRd1ChjD ), but it runs three times slower than the stl based solution.
1. Check this out: "Welford's online algorithm". It also mentions "numerical instability" which must be considered especially if a huge amount of elements is calculated. 2. the variance_k in this answer does not mean the real variance. Rather, it's the sum of squares of deviations, i.e. sum((x_k - mean)^2). 3. One can loop from k=1(i=0 if one uses index instead) i.e. the 1st element. 4. @DarioP tiny mistake(it should be k+1 instead of k).
1

My answer is similar as Josh Greifer but generalised to sample covariance. Sample variance is just sample covariance but with the two inputs identical. This includes Bessel's correlation.

 template <class Iter> typename Iter::value_type cov(const Iter &x, const Iter &y) { double sum_x = std::accumulate(std::begin(x), std::end(x), 0.0); double sum_y = std::accumulate(std::begin(y), std::end(y), 0.0); double mx = sum_x / x.size(); double my = sum_y / y.size(); double accum = 0.0; for (auto i = 0; i < x.size(); i++) { accum += (x.at(i) - mx) * (y.at(i) - my); } return accum / (x.size() - 1); }
Improve this answer

Comments

0

2x faster than the versions before mentioned - mostly because transform() and inner_product() loops are joined. Sorry about my shortcut/typedefs/macro: Flo = float. CR const ref. VFlo - vector. Tested in VS2010

#define fe(EL, CONTAINER) for each (auto EL in CONTAINER) //VS2010 Flo stdDev(VFlo CR crVec) { SZ n = crVec.size(); if (n < 2) return 0.0f; Flo fSqSum = 0.0f, fSum = 0.0f; fe(f, crVec) fSqSum += f * f; // EDIT: was Cit(VFlo, crVec) { fe(f, crVec) fSum += f; Flo fSumSq = fSum * fSum; Flo fSumSqDivN = fSumSq / n; Flo fSubSqSum = fSqSum - fSumSqDivN; Flo fPreSqrt = fSubSqSum / (n - 1); return sqrt(fPreSqrt); }
Improve this answer

4 Comments

Can the Cit() loop be written as for( float f : crVec ) { fSqSum += f * f; fSum += f; } ?
Yes in C++11. Trying to use macros that make it version independent. Updated the code. PS. For readability I usually prefer 1 action per LOC. Compiler should see that those are constant iterations and join them if it "thinks" it's faster to iterate once. Doing it in small short steps (without using std::inner_product() e.g.), kind of assembly-style, explains to new reader what it means. Binary will be smaller by side-effect (in some cases).
"Trying to use macros that make it version independent" - yet you limit yourself to the non-standard Visual C++ "for each" construct (stackoverflow.com/questions/197375/…)
@codeling It's just 1 macro for an illustration for 1 version of C++ for that post only. The was the algorithm - not coding std. Back then I used even uglier Cit(CFlo, crVec), which had default const-iter "cit", but re-indicates the container type. List of all compiler/OS-specific macros, is good when portability is in questions. In examples with boost it's also not easy to port it to std C++. I didn't explain the ugly short Flo, VFlo, CR, SZ neither -> float, vector<float>, const&, size - for shorten iteration lines of std C++. Same style Crit(MSZPFlo, crMap) foo(*crit.second); //rev-iter
0

In order to calculate the sample mean with a better presicion the following r-step recursion can be used:

mean_k=1/k*[(k-r)*mean_(k-r) + sum_over_i_from_(n-r+1)_to_n(x_i)],

where r is chosen to make summation components closer to each other.

Improve this answer

Comments

-6

Create your own container:

template <class T> class statList : public std::list<T> { public: statList() : std::list<T>::list() {} ~statList() {} T mean() { return accumulate(begin(),end(),0.0)/size(); } T stddev() { T diff_sum = 0; T m = mean(); for(iterator it= begin(); it != end(); ++it) diff_sum += ((*it - m)*(*it -m)); return diff_sum/size(); } };

It does have some limitations, but it works beautifully when you know what you are doing.

Improve this answer

2 Comments

To answer the question: because there’s absolutely no need. Creating your own container has absolutely no benefits compared to writing a free function.
I don't even know where to start with this. You're using a list as the underlying data structure, you don't even cache the values, which would be one of the few reasons I can think of to use a container-like structure. Especially if the values chance infrequently and the mean/stddev are needed often.
-8

//means deviation in c++

/A deviation that is a difference between an observed value and the true value of a quantity of interest (such as a population mean) is an error and a deviation that is the difference between the observed value and an estimate of the true value (such an estimate may be a sample mean) is a residual. These concepts are applicable for data at the interval and ratio levels of measurement./

#include <iostream> #include <conio.h> using namespace std; /* run this program using the console pauser or add your own getch, system("pause") or input loop */ int main(int argc, char** argv) { int i,cnt; cout<<"please inter count:\t"; cin>>cnt; float *num=new float [cnt]; float *s=new float [cnt]; float sum=0,ave,M,M_D; for(i=0;i<cnt;i++) { cin>>num[i]; sum+=num[i]; } ave=sum/cnt; for(i=0;i<cnt;i++) { s[i]=ave-num[i]; if(s[i]<0) { s[i]=s[i]*(-1); } cout<<"\n|ave - number| = "<<s[i]; M+=s[i]; } M_D=M/cnt; cout<<"\n\n Average: "<<ave; cout<<"\n M.D(Mean Deviation): "<<M_D; getch(); return 0;

}

Improve this answer

Comments

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.