// Computes average and standard deviation 
// for data stored in std::vector 
#include <vector>
#include <algorithm>
#include <numeric>
#include <functional>
#include <cmath>

void get_stats(const std::vector<double> & x, double & average, double & sigma )
{
  // average = sum_i=1^N x_i /N
  //
  //              sum_i=1^N (x_i-<x>)^2          (x-<x>).(x-<x>)
  // sigma = sqrt( -------------------- ) = sqrt(---------------)
  //                    N-1                         N-1
  using namespace std;
  int N = x.size();  
  average  = accumulate(x.begin(),x.end(),0.0)/N;   // sum up values
  vector<double> xx(x); // xx = x-<x>; subtract average from all elements
  transform(x.begin(),x.end(),xx.begin(),bind2nd(minus<double>(),average));
  sigma = sqrt(inner_product(xx.begin(),xx.end(),xx.begin(),0.0)/(N-1));
}