mcsusefulmath.cpp

蒙特卡洛仿真源代码很有参考价值要按照步骤进行操作

/*
  These functions provide useful mathematical functions for other
  classes to use.  
*/

/*
    MCS-libre, a library to facilite Monte-Carlo simulation.
    Copyright (C) 2002  D. E. Williams

    This library is free software; you can redistribute it and/or
    modify it under the terms of the GNU Lesser General Public
    License as published by the Free Software Foundation; either
    version 2.1 of the License, or (at your option) any later version.

    This library is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
    Lesser General Public License for more details.

    You should have received a copy of the GNU Lesser General Public
    License along with this library; if not, write to the Free Software
    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

    The full text of the GNU Lesser General Public License is
    contained in the files COPYING and LICENSE.

    See the 'MCS-libre-HOWTO' in the 'doc' directory for author's
    contact information.  

    Have a nice day!
*/

#include <cmath>

#include "mcsUsefulMath.h"

#define PI 3.14159
#define TOL .00001

double sign(double x)
{
  if (x == 0)
    {
      return 0;
    }
  if (x > 0)
    {
      return 1;
    }
  if (x < 0)
    {
      return -1;
    }
}

double inverse_error_function(double x)
{
  double oldval = 1;
  double newval = sqrt(-1*log(1-fabs(x)))*sign(x);
  while (fabs(erf(newval)-x) > (TOL*fabs(x)))
    {
      oldval = newval;
      newval = oldval - (erf(oldval)-x)*exp(pow(oldval,2))*((sqrt(PI))/2);
    }
  return newval;
}

double inverse_cdf_stdnormal(double U)
{
  return sqrt(2) * inverse_error_function(2*U - 1);
}