Take the 2-minute tour ×

I am supposed to compute the standard deviation function in some monte carlo simulations. The formula is this one: enter image description here

I think my results are way off what they should be. My function uses tuples from the boost library and it looks like this:

double add_square(double prev_sum, double new_val)
  return prev_sum + new_val*new_val;

template <typename V>
double vec_add_squares(const V<double>& v)
  return std::accumulate(v.begin(), v.end(), 0.0, add_square);

    template <class T> 
    boost::tuple<double,double> get_std_dev_and_error(const vector<T>& input, double r, double N)
 double M = double(input.size());

 double sum = std::accumulate(input.begin(),input.end(),0.0);
 double Squared_sum = vec_add_squares(input);

 std::cout << "sum " << Squared_sum << endl;

 // Calls Sum
 double term1 = Squared_sum - (sum/M)*sum;

 double SD = (sqrt(term1) * exp(-2.0 * r *N))/(M-1) ;
 double SE = SD/sqrt(M);
 std::cout << "SD = " << SD << endl;
 std::cout << "SE = " << SE << endl;

 return boost::tuple<double,double>(SD, SE) ;
  1. Can anyone see any mistakes here?
  2. also, there is the "accumulate" funciton in the STL library - does there exist an accumulate squared (members of the container)?
share|improve this question
For accumulate you can write your own functor –  PSIAlt Nov 15 '12 at 19:14
Thanks for the answer PSIA. I did write my own as in the code above but I get incorrect results, so it must be wrong or wrongly implemented in the program. –  Mathias Nov 15 '12 at 19:57
Fyi, the error indicator's own error can be pretty huge. Getting confidence intervals for a monte-carlo simulation is harder than usually believed. –  Alexandre C. Nov 15 '12 at 22:02
Yes, it's quite big... the reason I know it's wrong is because I was told the results. –  Mathias Nov 15 '12 at 22:30
Maybe stats.stackexchange.com is a good place to ask about the computational formula? –  Riga Nov 16 '12 at 20:36

1 Answer 1

Just use Boost.Accumulators (as you already use boost):

#include <boost/accumulators/accumulators.hpp>
#include <boost/accumulators/statistics.hpp>
#include <boost/range/algorithm.hpp>
#include <iostream>
#include <ostream>

using namespace boost;
using namespace boost::accumulators;
using namespace std;

int main()
    accumulator_set<double, stats<tag::sum , tag::variance, tag::mean > > acc;
    double data[] = {1., 2., 3.};
    acc = for_each(data, acc);
    cout << "sum = " << sum(acc) << endl;
    cout << "variance = " << variance(acc) << endl;
    cout << "sqrt(variance()) = " << sqrt(variance(acc)) << endl;
    cout << "mean = " << mean(acc) << endl;

Output is:

sum = 6
variance = 0.666667
sqrt(variance()) = 0.816497
mean = 2
share|improve this answer
I really like Boost.Accumulators - I was a bit too fast about the standard deviation, it's not part of Accumulators (unfortunately)... I will add a picture of the formula –  Mathias Nov 15 '12 at 21:47
@Mathias: you have the variance available in boost::accumulators. Just take its square root (don't forget your extra factor exp (-2rT) / (M - 1)). Anyway, don't use the formula you posted -- it is terribly instable in many cases (boost uses a stable on-line fomula for the variance). –  Alexandre C. Nov 15 '12 at 22:01
When I try to re-write the equations from boost.org/doc/libs/1_52_0/doc/html/boost/accumulators/impl/… I still can't get the same formula as posted above (except for the exp(-2Tr) ) –  Mathias Nov 15 '12 at 22:29
@AlexandreC. sorry forgot to tag you –  Mathias Nov 15 '12 at 22:46
@Mathias: Your formula reads approx. SD = sqrt (variance) * exp (-2rT) / sqrt(M - 1). There should be M /(M - 1) in front of the first sum, and the denominator ought to be M instead of M - 1 to be exactly the standard stdev estimator that you would get out of boost::accumulators. This should not change much the result however; –  Alexandre C. Nov 15 '12 at 23:22

Your Answer


