Monte-Carlo simulation is a sometimes elegant (and sometimes crude) method for simulating complex systems. Parameters that affect the system are selected from random distributions and the system response to these values is then calculated. Repeating this process many times produces often useful information about the system. The method is especially useful for examining non-linear systems that are difficult to write equations for.

Monte-Carlo simulation proves useful in the following areas:

• Discrete-Event Simulation
• Systems Science
• Bayesian Inference

This list is by no means exhaustive. Please feel free to comment if your favorite application of Monte-Carlo simulation is missing.

MCS-libre is a C++ library to facilitate Monte-Carlo simulation that I wrote about ten years ago. I’m writing about its use now so that I can reference it when I demonstrate SWIG in a future post. The library is freely available at http://sourceforge.net/projects/mcs-libre. Usage information, with a full example, is given below:

Usage

Sampling

A class called mcsLibre facilites sampling during simulation. Begin by including the necessary header file:

#include "mcsLibre.h"

Then declare an object of type mcsLibre:

mcsLibre samplesource;

This initializes MCS-libre using the default 100 random number seeds. Users can provide an alternate set of seeds if they wish.

One can then sample from normal and uniform distributions across multiple random number streams.

double i = samplesource.normal(mean, stdev, stream);
double j = samplesource.uniform(min, max, stream);

The stream argument used above is an integer that specifies which random number stream to generate values from. It is imperative that statistically independent variables sample from different streams!

Collecting Simulation Statistics

A class called mcsLibreCollectStats assists with collection of statistical information during simulation. Begin by including the necessary header file:

#include “mcsLibreCollectStats.h”

Then declare an object of type mcsLibreCollectStats as follows:

mcsLibreCollectStats stats_var1;

To add random samples of a variable to the collection object, use the add_value function:

stats_var1.add_value(30.3);

Retrieving Statistical Information

The mcsLibreCollectStats class provides functions for extracting statistical information out of its objects:

double m1 = stats_var1.get_mean_value(); // mean
double s1 = stats_var1.get_stdev_value(); // standard deviation
double mn1 = stats_var1.get_min_value(); // minimum value
double mx1 = stats_var1.get_max_value(); // maximum value
double n1 = stats_var1.get_number_samples(); // number of samples
double sum1 = stats_var1.get_sum_of_samples(); // sum of samples

Example

The rest of this post demonstrates the use of MCS-libre with an example simulation.

Simulating Reynold’s Number

Reynold’s number is a non-dimensional value used to describe flow around an object in a fluid dynamics system. It is the ratio of inertial forces to viscous forces, and is described by the equation

where V is the velocity of the object relative to the fluid, x is the distance the fluid travels across the object, and nu is the kinematic viscosity of the fluid. In this example we sample values for V, x, and nu from different distributions and then calculate the Reynold’s number from the sampled values. We do this thousands of times, storing each result, to gain information about how the Reynold’s number value responds to the sampled variable space.

Full Implementation of Example

/* Demonstration of MCS-libre toolkit for computing propagation of error in the Reynolds number. */

#include <iostream>

// include mcsLibre simulation library
#include "mcsLibre.h"

// include mcsLibre statistical probe class
#include "mcsLibreCollectStats.h"

main()
{
  std::cout << "Performing Monte-Carlo simulation..." << std::endl;

  // create object for simulation
  mcsLibre reynolds;

  // create probe for data collection
  mcsLibreCollectStats reynolds_stats;

  // define parameter types
  double Re, Re_mean, Re_stdev;    // Reynolds number of flow
  double nu, nu_mean, nu_stdev;    // kinematic viscosity of fluid
  double v, v_mean, v_stdev;       // velocity of flow
  double x, x_mean, x_stdev;       // distance considered

  // set mean and standard deviations from measurements
  nu_mean = 0.0000157; nu_stdev = 0.000002;  //  m^2 / s
  v_mean = 2; v_stdev = 0.1;                 //  m / s
  x_mean = 12; x_stdev = 1;                  //  m

  // simulation loop
  int samples = 15000;      // number of simulations
  for (int loopvar = 0; loopvar < samples; loopvar++)
  {
    // simulate nu from stream 1
    nu = reynolds.normal(nu_mean, nu_stdev, 1);

    // simulate v from stream 2
    v = reynolds.normal(v_mean, v_stdev, 2);

    // simulate x from stream 3
    x = reynolds.normal(x_mean, x_stdev, 3);

    // calculate a sample Reynolds number
    Re = (v * x) / nu;

    // send sample value to data collection probe
    reynolds_stats.add_value(Re);
  }

  // retrieve mean value from simulation
  Re_mean = reynolds_stats.get_mean_value();

  // retrieve standard deviation value from simulation
  Re_stdev = reynolds_stats.get_stdev_value();

  std::cout << std::endl << "Simulation of Error Propagation " << "in Reynold's Number" << std::endl;
  std::cout << "Re = (v * x) / nu" << std::endl << std::endl;
  std::cout << "Let v = normal(" << v_mean << ", " << v_stdev << "), x = normal(" << x_mean << ", "  << x_stdev << ")," << std::endl;
  std::cout << "and nu = normal(" << nu_mean << ", " << nu_stdev << ")." << std::endl << std::endl;
  std::cout << "Simulating " << samples << " Reynold's numbers yields:" << std::endl << std::endl;
  std::cout << "Sample Mean = " << Re_mean << "   Sample Standard Deviation = " << Re_stdev << std::endl;
  std::cout << "Sample Max = " << reynolds_stats.get_max_value();
  std::cout << "    Sample Min = " << reynolds_stats.get_min_value() << std::endl;
  std::cout << "Sum of Samples = " << reynolds_stats.get_sum_of_values() << std::endl;
  std::cout << "Number of Samples = " << reynolds_stats.get_number_samples() << std::endl;
  std::cout << std::endl;
}

Example Results