Monte Carlo Simulation Example

The following example code is written in C#, but it's almost identical in C++, and very similar in F#, Visual Basic, VB.NET, Java, and MATLAB.  Complete example source codein all these languages is included with the Solver Platform SDK.

// Example program calling the Monte Carlo Simulation engine.
// Given two uncertain variables X and Y with distributions: 
// X = [automatically fitted to sample data] and Y = PsiNormal(1,1) 
// Perform 1000 trials and calculate the functions: X + Y and X * Y
private Engine_Action SimEvaluator_OnEvaluate(Evaluator evaluator)
{
   // First, we ask for the current uncertain variable values
   double [] Var = evaluator.Problem.VarUncertain.Value.Array;
   m_Log.Dump("Evaluation x1 = " + Var[0] + ", x2 = " + Var[1] + ", Trial = " +
      evaluator.Problem.Engine.Stat.FunctionEvals);

   // We compute the function values given the current variable values
   evaluator.Problem.FcnUncertain.Value[0] = Var[0] + Var[1];  // X+Y
   evaluator.Problem.FcnUncertain.Value[1] = Var[0] * Var[1];  // X*Y

   return Engine_Action.Continue;
}private void Example18()
{
   m_Log.Clear();
   try
   {
      int nvars = 2, nfcns = 2;
      double [] correl = new double [] { 1.0, 0.5, 0.5, 1.0 };
      double [] sample = new double [] // for distribution fitting
         {-3.44,0.69,-0.11,1.25,0.73,-1.66,1.57,-1.64,-0.51,-0.19,
           1.24,-0.96,1.40,-0.01,-0.77,0.91,-0.72,0.85,-0.20,0.59,
           1.65,1.40,0.19,-1.38,-0.45,-0.71,1.25,-1.18,-0.33,2.33,
           0.73,-0.75,-0.63,-1.06,-0.71,0.33,-0.82,0.43,0.10,0.88,
           0.11,-0.88,0.02,-0.94,0.05,-0.06,-1.87,-2.03,-0.89,1.14};

      using (Problem problem = new Problem(Solver_Type.Simulate, nvars, nfcns))
      {
         problem.Evaluators[Eval_Type.Sample].OnEvaluate += new
            EvaluateEventHandler(SimEvaluator_OnEvaluate);
         problem.Solver.NumTrials = 1000;

         Distribution dist = new Distribution();
         dist.InitSample(sample); // set up the sample
         dist.AutoFit(true, 0.95); // continuous data, 95% confidence
         if (dist.FitSuccess)
         {
            m_Log.Dump("Distribution Type = " + dist.Name);
            for (int i=0; i< dist.NumParams;i++)
               m_Log.Dump("Param[" + i + "] = " + dist.Params[i]);
         }
         else return; // fitting unsuccessful

         problem.Model.VarDistribution[0] = dist; // assign the fitted distribution
         problem.Model.VarDistribution[1] = // define an analytic distribution
            new Distribution("PsiNormal(1,1)");

         problem.Model.DistCorrelations = // define the correlations
            new DoubleMatrix(Array_Order.ByRow, nvars, nvars, correl);

         problem.Solver.Simulate(); // perform the simulation
         m_Log.Dump("Status = " + problem.Solver.SimulateStatus);

         // display simulation results
         m_Log.Dump("\r\nStatistics on Uncertain Variables");
         for (int i = 0; i < nvars; i++)
               m_Log.Dump("\r\nVariable " + i); 
               m_Log.Dump("Minimum = " + problem.VarUncertain.Statistics.Minimum[i]);
               m_Log.Dump("Maximum = " + problem.VarUncertain.Statistics.Maximum[i]);
               m_Log.Dump("Mean = " + problem.VarUncertain.Statistics.Mean[i]);
               m_Log.Dump("StdDev = " + problem.VarUncertain.Statistics.StdDev[i]);
               m_Log.Dump("Variance = " + problem.VarUncertain.Statistics.Variance[i]);
               m_Log.Dump("Skewness = " + problem.VarUncertain.Statistics.Skewness[i]);
               m_Log.Dump("Kurtosis = " + problem.VarUncertain.Statistics.Kurtosis[i]);
               m_Log.Dump("10th Percentile = " + problem.VarUncertain.Percentiles[i, 9]);
               m_Log.Dump("90th Percentile = " + problem.VarUncertain.Percentiles[i, 89]);

         m_Log.Dump("\r\nStatistics on Uncertain Functions");
         for (int i = 0; i < nfcns; i++)
               m_Log.Dump("\r\nFunction " + i);
               m_Log.Dump("Minimum = " + problem.FcnUncertain.Statistics.Minimum[i]);
               m_Log.Dump("Maximum = " + problem.FcnUncertain.Statistics.Maximum[i]);
               m_Log.Dump("Mean = " + problem.FcnUncertain.Statistics.Mean[i]);
               m_Log.Dump("StdDev = " + problem.FcnUncertain.Statistics.StdDev[i]);
               m_Log.Dump("Variance = " + problem.FcnUncertain.Statistics.Variance[i]);
               m_Log.Dump("Skewness = " + problem.FcnUncertain.Statistics.Skewness[i]);
               m_Log.Dump("Kurtosis = " + problem.FcnUncertain.Statistics.Kurtosis[i]);
               m_Log.Dump("10th Percentile = " + problem.FcnUncertain.Percentiles[i,9]);
               m_Log.Dump("90th Percentile = " + problem.FcnUncertain.Percentiles[i,89]);
      }
   }
   catch (SolverException ex)
   {
      m_Log.Dump("Exception " + ex.ExceptionType);
      m_Log.Dump("Exception " + ex.Message);
   }
}

< Back to Monte Carlo Simulation Summary

< Back to Solver Platform SDK Product Overview