Integration of Normal Distribution PDF using CAPD3 and CAPD4

Normal_Distribution_PDF_Integration_CAPD3_CAPD4.md

We want to integrate the probability density function for normal distribution:

To simplify the problem, let's pick mean = 0.0, standard deviation = 1.0. In theory, integration of this fuction from x = -oo to x = oo should give P = 1.0. Using CAPD, let's compute for the range X = [-10, 10] whose prob must be less than 1.0.

Output

CAPD3: [0.999999081962,1.00000075862]
CAPD4: [1.00000002444, 1.00000002444]

CAPD4's answer is not correct in a sense that its lower-bound is greater than 1.0.

Question

• Is it a bug or a consequence of not using multi-precision floating-point library?
• Is there any way to avoid this result (hopefully without using multi-precision libraries)?

Source Code

Here are the code that I used to produce the result. They are the modified version of CAPD's example code.

CAPD3

/////////////////////////////////////////////////////////////////////////////
//
/// @file encloseTrajectoryBetweenTimeSteps.cpp
///
/// @author Daniel Wilczak
//
/////////////////////////////////////////////////////////////////////////////
// Copyright (C) CAPD group
//
// This file constitutes a part of the CAPD library,
// distributed under the terms of the GNU General Public License.
// Consult  http://capd.ii.uj.edu.pl/ for details.

#include "capd/capdlib.h"
using namespace capd;
using namespace std;
int main() {
  try {
    cout.precision(12);
    //
    // probability density function = 1 / (sigma * sqrt(2pi)) * exp(-1/2 ((x-mu)/sigma)^2)
    // where mu = 0, sigma = 1.0
    //
    IMap vectorField("var:p, x;fun:(1.0/sqrt(6.283185))*exp(-1/2.0*(x^2.0)), 1.0;");
    ITaylor solver(vectorField,20,.1);
    ITimeMap timeMap(solver);
    IVector x(2);
    x[0]=0.0;   // p
    x[1]=-10.0; // x
    C0Rect2Set s(x);
    double T=20;
    timeMap.stopAfterStep(true);
    interval prevTime(0.);
    do {
      timeMap(T,s);
      interval stepMade = solver.getStep();
      cout << "\nstep made: " << stepMade;
      const ITaylor::CurveType& curve = solver.getCurve();
      interval domain = interval(0,1)*stepMade;
      int grid=5;
      for(int i=0;i<grid;++i)
      {
        interval subsetOfDomain = interval(i,i+1)*stepMade/grid;
        intersection(domain,subsetOfDomain,subsetOfDomain);
        IVector v = curve(subsetOfDomain);
        std::cout << "\nenclosure for t=" << prevTime + subsetOfDomain << ":  " << v;      
      }
      prevTime = timeMap.getCurrentTime();
      cout << "\ncurrent time: " << prevTime << endl << endl;
    } while(!timeMap.completed());
  } catch(exception& e) {
    cout << "\n\nException caught!\n" << e.what() << endl << endl;
  }
}

CAPD4

/////////////////////////////////////////////////////////////////////////////
//
/// @file encloseTrajectoryBetweenTimeSteps.cpp
///
/// @author Daniel Wilczak
//
/////////////////////////////////////////////////////////////////////////////
// Copyright (C) CAPD group
//
// This file constitutes a part of the CAPD library,
// distributed under the terms of the GNU General Public License.
// Consult  http://capd.ii.uj.edu.pl/ for details.

#include "capd/capdlib.h"
using namespace capd;
using namespace std;
int main() {
  try {
    cout.precision(12);
    //
    // probability density function = 1 / (sigma * sqrt(2pi)) * exp(-1/2 ((x-mu)/sigma)^2)
    // where mu = 0, sigma = 1.0
    //
    IMap vectorField("var:p, x;fun:(1.0/sqrt(6.283185))*exp(-1/2.0*(x^2.0)), 1.0;");
    IOdeSolver solver(vectorField,20);
    ITimeMap timeMap(solver);
    IVector x(2);
    x[0]=0.0;
    x[1]=-10.0;
    C0Rect2Set s(x);
    double T=20;
    timeMap.stopAfterStep(true);
    interval prevTime(0.);
    do {
      timeMap(T,s);
      interval stepMade = solver.getStep();
      cout << "\nstep made: " << stepMade;
      const IOdeSolver::SolutionCurve& curve = solver.getCurve();
      interval domain = interval(0,1)*stepMade;
      int grid=5;
      for(int i=0;i<grid;++i) {
        interval subsetOfDomain = interval(i,i+1)*stepMade/grid;
        intersection(domain,subsetOfDomain,subsetOfDomain);
        IVector v = curve(subsetOfDomain);
        std::cout << "\nenclosure for t=" << prevTime + subsetOfDomain << ":  " << v;      
      }
      prevTime = timeMap.getCurrentTime();
      cout << "\ncurrent time: " << prevTime << endl << endl;

    } while(!timeMap.completed());
  } catch(exception& e) {
    cout << "\n\nException caught!\n" << e.what() << endl << endl;
  }
}  // END

Hello,
your test contains a bug because 6.283185 is not an interval enclosing 2*PI. Use parameters instead:

IMap vectorField("par:pi;var:p, x;fun:(1.0/sqrt(pi))*exp(-1/2.0*(x^2.0)), 1.0;");
vectorField.setParameter("pi",interval::pi() * 2.0);

Daniel Wilczak