/*
    
   NB requires C++ V11 standard.
   
   http://www.cplusplus.com/reference/random/uniform_real_distribution/
   http://www.cplusplus.com/reference/random/normal_distribution/
*/
 
#include <iostream>
#include <random>

using namespace std;

int main()
{
   default_random_engine re;
//   mt19937 re;                // Mersenne Twister

// Uniform distribution P(x)=1/(b-a) for a<=x<=b and 0 elsewhere.  
   double a = -1.0;
   double b =  1.0;
   uniform_real_distribution<double> distribution(a,b);


// Code for a Normal AKA Gaussian distribution P(X) \propto Exp( -(X-a)^2/(2 \sigma^2) )
/*
   double average=5.0;
   double standarddev=2.0;
   normal_distribution<double> distribution( average, standarddev);
*/

   const int N=10000;
   double M1=0;
   double M2=0;
   double M3=0;
   double M4=0;

   for (int i=0 ; i<N ; i++)
     {
        double sample = distribution(re);
        
        M1+=sample;
        M2+=sample*sample;
        M3+=pow(sample,3);
        M4+=pow(sample,4);
     }

   M1/=N;  M2/=N;   M3/=N;   M4/=N;

   double avg=M1;                 // Average
   double var=M2-M1*M1;           // Variance
   double std=sqrt(var);          // standard deviation
   
   cout << "Avg=<X>="           << avg << endl;
   cout << "Var=<X^2>-<X>^2="   << var << endl;
   cout << "Std.dev=sqrt(Var)=" << std << endl;
   cout << "<X^2>="             << M2 << endl;
   cout << "<X^3>="             << M3 << endl;
   cout << "<X^4>="             << M4 << endl;

   return 0;
}