Logo Search packages:      
Sourcecode: openturns version File versions  Download package

DistributionImplementation.cxx

Go to the documentation of this file.
//                                               -*- C++ -*-
/**
 *  @file  DistributionImplementation.cxx
 *  @brief Abstract top-level class for all distributions
 *
 *  (C) Copyright 2005-2007 EDF-EADS-Phimeca
 *
 *  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.
 *
 *  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
 *
 *  @author: $LastChangedBy: dutka $
 *  @date:   $LastChangedDate: 2008-10-31 11:52:04 +0100 (ven 31 oct 2008) $
 *  Id:      $Id: DistributionImplementation.cxx 995 2008-10-31 10:52:04Z dutka $
 */
#include <cmath>
#include "PersistentObjectFactory.hxx"
#include "DistributionImplementation.hxx"
#include "Distribution.hxx"
#include "Exception.hxx"
#include "Curve.hxx"
#include "Contour.hxx"
#include "Log.hxx"
#include "Lapack.hxx"
#include "SquareMatrix.hxx"
#include "Collection.hxx"
#include "RandomGenerator.hxx"
#include "Box.hxx"
#include "Normal.hxx"
#include "IndependentCopula.hxx"
#include "MarginalTransformationEvaluation.hxx"
#include "MarginalTransformationGradient.hxx"
#include "MarginalTransformationHessian.hxx"
#include "InverseMarginalTransformationEvaluation.hxx"
#include "InverseMarginalTransformationGradient.hxx"
#include "InverseMarginalTransformationHessian.hxx"
#include "NatafIndependentCopulaEvaluation.hxx"
#include "NatafIndependentCopulaGradient.hxx"
#include "NatafIndependentCopulaHessian.hxx"
#include "InverseNatafIndependentCopulaEvaluation.hxx"
#include "InverseNatafIndependentCopulaGradient.hxx"
#include "InverseNatafIndependentCopulaHessian.hxx"
#include "RosenblattEvaluation.hxx"
#include "InverseRosenblattEvaluation.hxx"
#include "NumericalMathFunctionImplementation.hxx"
#include "SklarCopula.hxx"

namespace OpenTURNS {

  namespace Uncertainty {

    namespace Model {

      typedef Base::Common::NotYetImplementedException                 NotYetImplementedException;
      typedef Base::Common::InvalidArgumentException                   InvalidArgumentException;
      typedef Base::Common::Log                                        Log;
      typedef Base::Graph::Curve                                       Curve;
      typedef Base::Graph::Contour                                     Contour;
      typedef Base::Type::SquareMatrix                                 SquareMatrix;
      typedef Base::Type::Collection<Distribution>                     DistributionCollection;
      typedef Base::Stat::RandomGenerator                              RandomGenerator;
      typedef Algorithm::Box                                           Box;
      typedef Uncertainty::Distribution::Normal                        Normal;
      typedef Uncertainty::Distribution::IndependentCopula             IndependentCopula;
      typedef Uncertainty::Distribution::SklarCopula                   SklarCopula;
      typedef Algorithm::NatafIndependentCopulaEvaluation              NatafIndependentCopulaEvaluation;
      typedef Algorithm::NatafIndependentCopulaGradient                NatafIndependentCopulaGradient;
      typedef Algorithm::NatafIndependentCopulaHessian                 NatafIndependentCopulaHessian;
      typedef Algorithm::InverseNatafIndependentCopulaEvaluation       InverseNatafIndependentCopulaEvaluation;
      typedef Algorithm::InverseNatafIndependentCopulaGradient         InverseNatafIndependentCopulaGradient;
      typedef Algorithm::InverseNatafIndependentCopulaHessian          InverseNatafIndependentCopulaHessian;
      typedef Algorithm::MarginalTransformationEvaluation              MarginalTransformationEvaluation;
      typedef Algorithm::MarginalTransformationGradient                MarginalTransformationGradient;
      typedef Algorithm::MarginalTransformationHessian                 MarginalTransformationHessian;
      typedef Algorithm::InverseMarginalTransformationEvaluation       InverseMarginalTransformationEvaluation;
      typedef Algorithm::InverseMarginalTransformationGradient         InverseMarginalTransformationGradient;
      typedef Algorithm::InverseMarginalTransformationHessian          InverseMarginalTransformationHessian;
      typedef Algorithm::RosenblattEvaluation                          RosenblattEvaluation;
      typedef Algorithm::InverseRosenblattEvaluation                   InverseRosenblattEvaluation;
      typedef Base::Func::NumericalMathFunctionImplementation          NumericalMathFunctionImplementation;
      typedef MarginalTransformationEvaluation::DistributionCollection DistributionCollection;

      const UnsignedLong    DistributionImplementation::DefaultLevelNumber;
      const UnsignedLong    DistributionImplementation::DefaultIntegrationNodesNumber;
      const UnsignedLong    DistributionImplementation::DefaultPointNumber;
      const NumericalScalar DistributionImplementation::DefaultQuantileEpsilon;
      const UnsignedLong    DistributionImplementation::DefaultQuantileIteration;
      const NumericalScalar DistributionImplementation::DefaultPDFEpsilon;
      const NumericalScalar DistributionImplementation::DefaultCDFEpsilon;
      const NumericalScalar DistributionImplementation::QMin;
      const NumericalScalar DistributionImplementation::QMax;

      CLASSNAMEINIT(DistributionImplementation);

      static Base::Common::Factory<DistributionImplementation> RegisteredFactory("DistributionImplementation");

      /* Default constructor */
00108       DistributionImplementation::DistributionImplementation(const String & name)
      : PersistentObject(name),
        mean_(NumericalPoint(0)),
        covariance_(CovarianceMatrix(0)),
        gaussNodesAndWeights_(),
        integrationNodesNumber_(DefaultIntegrationNodesNumber),
        isAlreadyComputedMean_(false),
        isAlreadyComputedCovariance_(false),
        isAlreadyComputedGaussNodesAndWeights_(false),
        pdfEpsilon_(DefaultPDFEpsilon),
        cdfEpsilon_(DefaultCDFEpsilon),
        quantileEpsilon_(),
        dimension_(1),
        weight_(1.),
        range_(),
        description_(1)
      {
      // Initialize any other class members here
      // At last, allocate memory space if needed, but go to destructor to free it
      description_[0] = "marginal 1";
      }

      /* Virtual constructor */
00131       DistributionImplementation * DistributionImplementation::clone() const
      {
      return new DistributionImplementation(*this);
      }

      /* Comparison operator */
00137       Bool DistributionImplementation::operator ==(const DistributionImplementation & other) const {
      Bool sameObject = false;

      if (this != &other) { // Other is NOT me, so I have to realize the comparison
        // sameObject = ...
      } else sameObject = true;

      return sameObject;
      }
  
      /* String converter */
00148       String DistributionImplementation::str() const {
      OSS oss;
      oss << "class=" << DistributionImplementation::GetClassName()
            << " description=" << description_;
      return oss;
      }
  

      /* Weight accessor */
00157       void DistributionImplementation::setWeight(const NumericalScalar w)
      throw(InvalidArgumentException)
      {
      weight_ = w;
      }

      /* Weight accessor */
      NumericalScalar DistributionImplementation::getWeight() const
      {
      return weight_;
      }
      

      /* Dimension accessor */
00171       UnsignedLong DistributionImplementation::getDimension() const
      {
      return dimension_;
      }

      /* Get the roughness, i.e. the L2-norm of the PDF */
00177       NumericalScalar DistributionImplementation::getRoughness() const
      {
      throw NotYetImplementedException(HERE);
      }
      
      /* Dimension accessor */
00183       void DistributionImplementation::setDimension(const UnsignedLong dim)
      throw(InvalidArgumentException)
      {
      if (dim == 0) throw InvalidArgumentException(HERE) << "Dimension argument must be an integer >= 1, here dim = " << dim;
      if (dim != dimension_)
        {
          dimension_ = dim;
          isAlreadyComputedMean_ = false;
          isAlreadyComputedCovariance_ = false;
          isAlreadyComputedGaussNodesAndWeights_ = false;
          // Check if the current description is compatible with the new dimension
          if (description_.getSize() != dim)
            {
            description_ = Description(dim);
            for (UnsignedLong i = 0; i < dim; ++i)
              {
                description_[i] = OSS() << "marginal " << i+1;
              }
            }
        }
      }

      /* Get one realization of the distributionImplementation */
00206       DistributionImplementation::NumericalPoint DistributionImplementation::getRealization() const
      {
      // Use CDF inversion
      return computeQuantile(RandomGenerator::Generate());
      }
     
      /* Get a numerical sample whose elements follow the distributionImplementation */
00213       DistributionImplementation::NumericalSample DistributionImplementation::getNumericalSample(const UnsignedLong size) const
      {
      NumericalSample returnSample(size, dimension_);
      for (UnsignedLong i = 0; i < size; ++i)
        {
          NumericalPoint realization(getRealization());
          for (UnsignedLong j = 0; j < dimension_; ++j)
            {
            returnSample[i][j] = realization[j];
            }
        }
      returnSample.setName(getName());
      returnSample.setDescription(getDescription());
      return returnSample;
      }

      /* Get the DDF of the distributionImplementation */
      DistributionImplementation::NumericalPoint DistributionImplementation::computeDDF(const NumericalPoint & point) const
      {
      NumericalPoint ddf(dimension_);
      NumericalScalar cdfPoint(computeCDF(point));
      NumericalScalar idenom(1.0 / sqrt(cdfEpsilon_));
      for (UnsignedLong i = 0; i < dimension_; ++i)
        {
          NumericalPoint epsilon(dimension_);
          epsilon[i] = pow(cdfEpsilon_, 0.25);
          ddf[i] = (computeCDF(point + epsilon) - 2.0 * cdfPoint + computeCDF(point - epsilon)) * idenom;
        }
      return ddf;
      }

      /* Get the PDF of the distributionImplementation */
      NumericalScalar DistributionImplementation::computePDF(const NumericalPoint & point) const
      {
      NumericalPoint epsilon(dimension_, pow(cdfEpsilon_, 1.0 / 3.0));
      pdfEpsilon_ = pow(cdfEpsilon_, 2.0 / 3.0);
      // Centered finite differences of CDF
      return (computeCDF(point + epsilon) - computeCDF(point - epsilon)) / DefaultPDFEpsilon;
      }

      /* Get the CDF of the distributionImplementation */
      NumericalScalar DistributionImplementation::computeCDF(const NumericalPoint & point, const Bool tail) const
      {
      throw NotYetImplementedException(HERE) << "in DistributionImplementation::computeCDF(const NumericalPoint & point) const";
      }

      /* Compute the probability content of an interval */
00260       NumericalScalar DistributionImplementation::computeProbability(const Interval & interval) const
      {
      if (dimension_ == 1)
        {
          // Empty interval
          if (interval.isEmpty()) return 0.0;
          NumericalScalar upperCDF(1.0);
          NumericalScalar lowerCDF(0.0);
          if (interval.getFiniteLowerBound()[0]) lowerCDF = computeCDF(interval.getLowerBound());
          if (interval.getFiniteUpperBound()[0]) upperCDF = computeCDF(interval.getUpperBound());
          return upperCDF - lowerCDF;
        }
      throw NotYetImplementedException(HERE) << "in DistributionImplementation::computeProbability(const Interval & interval) const";
      }

      /* Get the characteristic function of the distribution, i.e. phi(u) = E(exp(I*u*X)) */
00276       NumericalComplex DistributionImplementation::computeCharacteristicFunction(const NumericalPoint & point) const
      {
      throw NotYetImplementedException(HERE) << "in DistributionImplementation::computeCharacteristicFunction(const NumericalPoint & point) const";
      }

      NumericalComplex DistributionImplementation::computeCharacteristicFunction(const NumericalScalar x) const
      {
      if (dimension_ != 1) throw InvalidDimensionException(HERE) << "Error:  cannot use the simplified interface of computeCharacteristicFunction with distributions of dimension > 1";
      return computeCharacteristicFunction(NumericalPoint(1, x));
      }

00287       NumericalComplex DistributionImplementation::computeCharacteristicFunction(const UnsignedLong index, const NumericalScalar step) const
      {
      if (dimension_ != 1) throw InvalidDimensionException(HERE) << "Error:  cannot use the simplified interface of computeCharacteristicFunction with distributions of dimension > 1";
      return computeCharacteristicFunction(NumericalPoint(1, index * step));
      }

      /* Get the DDF of the distributionImplementation */
      DistributionImplementation::NumericalSample DistributionImplementation::computeDDF(const NumericalSample & inSample) const
      {
      UnsignedLong size = inSample.getSize();
      NumericalSample outSample(size, 1);
      for(UnsignedLong i = 0; i < size; ++i) outSample[i] = computeDDF(inSample[i]);
      return outSample;
      }

      /* Get the PDF of the distributionImplementation */
      DistributionImplementation::NumericalSample DistributionImplementation::computePDF(const NumericalSample & inSample) const
      {
      UnsignedLong size = inSample.getSize();
      NumericalSample outSample(size, 1);
      for(UnsignedLong i = 0; i < size; ++i) outSample[i] = NumericalPoint(1, computePDF(inSample[i]));
      return outSample;
      }

      /* Get the CDF of the distributionImplementation */
      DistributionImplementation::NumericalSample DistributionImplementation::computeCDF(const NumericalSample & inSample, const Bool tail) const
      {
      UnsignedLong size = inSample.getSize();
      NumericalSample outSample(size, getDimension());
      for(UnsignedLong i = 0; i < size; ++i) outSample[i] = NumericalPoint(1, computeCDF(inSample[i], tail));
      return outSample;
      }

      /* Get the DDF of the distributionImplementation */
00321       NumericalScalar DistributionImplementation::computeDDF(const NumericalScalar scalar) const
      {
      if (dimension_ != 1) throw InvalidDimensionException(HERE) << "ERROR: cannot use the simplified interface of computeDDF with distributions of dimension > 1";
      return computeDDF(NumericalPoint(1, scalar))[0];
      }

      /* Get the PDF of the distributionImplementation */
00328       NumericalScalar DistributionImplementation::computePDF(const NumericalScalar scalar) const
      {
      if (dimension_ != 1) throw InvalidDimensionException(HERE) << "ERROR: cannot use the simplified interface of computePDF with distributions of dimension > 1";
      return computePDF(NumericalPoint(1, scalar));
      }

      /* Get the CDF of the distributionImplementation */
00335       NumericalScalar DistributionImplementation::computeCDF(const NumericalScalar scalar, const Bool tail) const
      {
      if (dimension_ != 1) throw InvalidDimensionException(HERE) << "ERROR: cannot use the simplified interface of computeCDF with distributions of dimension > 1";
      return computeCDF(NumericalPoint(1, scalar), tail);
      }

      /*  Compute the PDF of 1D distributions over a regular grid */
00342       DistributionImplementation::NumericalSample DistributionImplementation::computePDF(const NumericalScalar xMin, const NumericalScalar xMax, const UnsignedLong pointNumber) const
      {
      if (dimension_ != 1) throw InvalidArgumentException(HERE) << "Error: cannot compute the PDF over a regular 1D grid if the dimension is > 1";
      NumericalSample result(pointNumber, 2);
      NumericalScalar x(xMin);
      NumericalScalar step((xMax - xMin) / NumericalScalar(pointNumber - 1.0));
      for (UnsignedLong i = 0; i < pointNumber; ++i)
        {
          result[i][0] = x;
          result[i][1] = computePDF(x);
          x += step;
        }
      return result;
      }

      /*  Compute the CDF of 1D distributions over a regular grid */
00358       DistributionImplementation::NumericalSample DistributionImplementation::computeCDF(const NumericalScalar xMin, const NumericalScalar xMax, const UnsignedLong pointNumber, const Bool tail) const
      {
      if (dimension_ != 1) throw InvalidArgumentException(HERE) << "Error: cannot compute the CDF over a regular 1D grid if the dimension is > 1";
      NumericalSample result(pointNumber, 2);
      NumericalScalar x(xMin);
      NumericalScalar step((xMax - xMin) / NumericalScalar(pointNumber - 1.0));
      for (UnsignedLong i = 0; i < pointNumber; ++i)
        {
          result[i][0] = x;
          result[i][1] = computeCDF(x, tail);
          x += step;
        }
      return result;
      }

      /* Get the PDF gradient of the distribution */
00374       DistributionImplementation::NumericalPoint DistributionImplementation::computePDFGradient(const NumericalPoint & point) const
      {
      throw NotYetImplementedException(HERE) << "in DistributionImplementation::computeDerivedPDF(const NumericalPoint & point) const";
      }

      /* Get the CDF gradient of the distribution */
00380       DistributionImplementation::NumericalPoint DistributionImplementation::computeCDFGradient(const NumericalPoint & point) const
      {
      throw NotYetImplementedException(HERE) << "in DistributionImplementation::computeDerivedCDF(const NumericalPoint & point) const";
      }


      /* Compute the DDF of Xi | X1, ..., Xi-1. x = Xi, y = (X1,...,Xi-1) */
00387       NumericalScalar DistributionImplementation::computeConditionalDDF(const NumericalScalar x, const NumericalPoint & y) const
      {
      throw NotYetImplementedException(HERE) << "in DistributionImplementation::computeDerivedCDF(const NumericalPoint & point) const";
      }
      
      /* Compute the PDF of Xi | X1, ..., Xi-1. x = Xi, y = (X1,...,Xi-1) */
00393       NumericalScalar DistributionImplementation::computeConditionalPDF(const NumericalScalar x, const NumericalPoint & y) const
      {
      const UnsignedLong conditioningDimension(y.getDimension());
      if (conditioningDimension >= getDimension()) throw InvalidArgumentException(HERE) << "Error: cannot compute a conditional PDF with a conditioning point of dimension greater or equal to the distribution dimension.";
      // Special case for no conditioning or independent copula
      if ((conditioningDimension == 0) || (hasIndependentCopula())) return getMarginal(conditioningDimension)->computePDF(x);
      // General case
      Indices conditioning(conditioningDimension);
      for (UnsignedLong i = 0; i < conditioningDimension; ++i)
        {
          conditioning[i] = i;
        }
      Indices conditioned(conditioning);
      conditioned.add(conditioningDimension);
      Implementation conditioningDistribution(getMarginal(conditioning));
      const NumericalScalar pdfConditioning(conditioningDistribution->computePDF(y));
      if (pdfConditioning <= 0.0) return 0.0;
      NumericalPoint z(y);
      z.add(x);
      Implementation conditionedDistribution(getMarginal(conditioned));
      const NumericalScalar pdfConditioned(conditionedDistribution->computePDF(z));
      pdfEpsilon_ = conditionedDistribution->getPDFEpsilon() + conditioningDistribution->getPDFEpsilon();
      return pdfConditioned / pdfConditioning;
      }
      
      /* Compute the CDF of Xi | X1, ..., Xi-1. x = Xi, y = (X1,...,Xi-1) */
00419       NumericalScalar DistributionImplementation::computeConditionalCDF(const NumericalScalar x, const NumericalPoint & y) const
      {
      throw NotYetImplementedException(HERE) << "in DistributionImplementation::computeDerivedCDF(const NumericalPoint & point) const";
      }
      
      /* Compute the quantile of Xi | X1, ..., Xi-1, i.e. x such that CDF(x|y) = q with x = Xi, y = (X1,...,Xi-1) */
00425       NumericalScalar DistributionImplementation::computeConditionalQuantile(const NumericalScalar q, const NumericalPoint & y) const
      {
      const UnsignedLong conditioningDimension(y.getDimension());
      if (conditioningDimension >= dimension_) throw InvalidArgumentException(HERE) << "Error: cannot compute a conditional quantile with a conditioning point of dimension greater or equal to the distribution dimension.";
      if ((q < 0.0) || (q > 1.0)) throw InvalidArgumentException(HERE) << "Error: cannot compute a conditional quantile for a probability level outside of [0, 1]";
      // Initialize the conditional quantile with the quantile of the i-th marginal distribution
      Implementation marginalDistribution(getMarginal(conditioningDimension));
      NumericalScalar quantile(marginalDistribution->computeQuantile(q)[0]);
      // Special case for bording values
      if ((q == 0.0) || (q == 1.0)) return quantile;
      // Special case when no contitioning or independent copula
      if ((conditioningDimension == 0) || hasIndependentCopula()) return quantile;
      NumericalScalar step(marginalDistribution->getDispersionIndicator());
      // Search a bracketing interval [a, b] for the quantile
      NumericalScalar a(quantile);
      NumericalScalar aCDF(computeConditionalCDF(a, y));
      NumericalScalar b(a);
      NumericalScalar bCDF(aCDF);
      // While b is a lower bound for the quantile
      while (bCDF <= q)
        {
          b += step;
          bCDF = computeConditionalCDF(b, y);
          step *= 2.0;
        }
      // If b > a, then [b - 0.5 * step, b] is a bracketing interval
      if (b > a)
        {
          a = b - 0.5 * step;
        }
      else
        {
          // Here, we know that bCDF = aCDF > q, so the following loop is visited at least once
          while (aCDF > q)
            {
            a -= step;
            aCDF = computeConditionalCDF(a, y);
            step *= 2.0;
            }
          // The bracketing interval is [a, a + 0.5 * step]
          b = a + 0.5 * step;
        }
      // First, we try the Newton method initialized at 0.5 * (a + b)
      quantile = 0.5 * (a + b);
      Bool convergence(false);
      NumericalScalar residual(0.0);
      UnsignedLong iteration(0);
      while (!convergence && (iteration < DefaultQuantileIteration))
        {
          NumericalScalar pdf(computeConditionalPDF(quantile, y));
          // f pdf is zero, stop Newton iteration
          if (pdf == 0.0) break;
          NumericalScalar cdf(computeConditionalCDF(quantile, y));
          residual = (q - cdf) / pdf;
          quantile += residual;
          // Safety stop if the iteration seems to diverge
          if ((quantile > b) || (quantile < a)) break;
          convergence = fabs(residual) < DefaultQuantileEpsilon * (1.0 + fabs(quantile)) || (fabs(cdf - q) < 2.0 * cdfEpsilon_);
          ++iteration;
        }
      // If no convergence after Newton steps, try bisection method
      while (! convergence)
        {
          quantile = 0.5 * (a + b);
          NumericalScalar cdf(computeConditionalCDF(quantile, y));
          if (cdf > q)
            {
            b = quantile;
            }
          else
            {
            a = quantile;
            }
          convergence = (b - a < DefaultQuantileEpsilon * (1.0 + fabs(quantile))) || (fabs(cdf - q) < 2.0 * cdfEpsilon_);
        }
      return quantile;
      }


      /* Quantile computation for dimension=1 */
00505       NumericalScalar DistributionImplementation::computeScalarQuantile(const NumericalScalar prob, const NumericalScalar initialGuess, const NumericalScalar initialStep) const
      {
      if (getDimension() != 1) throw InvalidDimensionException(HERE) << "Error: the method computeScalarQuantile is only defined for 1D distributions";
      if ((prob < 0.0) || (prob > 1.0)) throw InvalidArgumentException(HERE) << "Error: cannot compute a quantile for a probability level outside of [0, 1]";
      // Special case for boarding values
      if (prob == 0.0) return range_.getLowerBound()[0];
      if (prob == 1.0) return range_.getUpperBound()[0];
      NumericalScalar quantile(initialGuess);
      NumericalScalar step(initialStep);
      // Search a bracketing interval [a, b] for the quantile
      NumericalScalar a(quantile);
      NumericalScalar aCDF(computeCDF(NumericalPoint(1, a)));
      NumericalScalar b(a);
      NumericalScalar bCDF(aCDF);
      // While b is a lower bound for the quantile
      while (bCDF <= prob)
        {
          b += step;
          bCDF = computeCDF(b);
          step *= 2.0;
        }
      // If b > a, then [b - 0.5 * step, b] is a bracketing interval
      if (b > a)
        {
          a = b - 0.5 * step;
        }
      else
        {
          // Here, we know that bCDF = aCDF > prob, so the following loop is visited at least once
          while (aCDF > prob)
            {
            a -= step;
            aCDF = computeCDF(a);
            step *= 2.0;
            }
          // The bracketing interval is [a, a + 0.5 * step]
          b = a + 0.5 * step;
        }
      // First, we try the Newton method initialized at 0.5 * (a + b)
      quantile = 0.5 * (a + b);
      Bool convergence(false);
      NumericalScalar residual(0.0);
      UnsignedLong iteration(0);
      while (!convergence && (iteration < DefaultQuantileIteration))
        {
          NumericalScalar pdf(computePDF(quantile));
          // f pdf is zero, stop Newton iteration
          if (pdf == 0.0) break;
          NumericalScalar cdf(computeCDF(quantile));
          residual = (prob - cdf) / pdf;
          quantile += residual;
          // Safety stop if the iteration seems to diverge
          if ((quantile > b) || (quantile < a)) break;
          convergence = fabs(residual) < DefaultQuantileEpsilon * (1.0 + fabs(quantile)) || (fabs(cdf - prob) < 2.0 * cdfEpsilon_);
          ++iteration;
        }
      // If no convergence after Newton steps, try bisection method
      while (! convergence)
        {
          quantile = 0.5 * (a + b);
          NumericalScalar cdf(computeCDF(quantile));
          if (cdf > prob)
            {
            b = quantile;
            }
          else
            {
            a = quantile;
            }
          convergence = (b - a < DefaultQuantileEpsilon * (1.0 + fabs(quantile))) || (fabs(cdf - prob) < 2.0 * cdfEpsilon_);
        }
      return quantile;
      } // computeScalarQuantile
      
      /* Generic implementation of the quantile computation */
00580       DistributionImplementation::NumericalPoint DistributionImplementation::computeQuantile(const NumericalScalar prob) const
      {
      if (prob < 0.0 || prob > 1.0) throw InvalidArgumentException(HERE) << "Error: cannot compute a quantile for a probability level outside of [0, 1]";
      // Special case for bording values
      if (prob == 0.0) return range_.getLowerBound();
      if (prob == 1.0) return range_.getUpperBound();
      // Special case for dimension 1
      if (dimension_ == 1) return NumericalPoint(1, computeScalarQuantile(prob));
      // Extract the marginal distributions
      Base::Type::Collection<Distribution> marginals(dimension_);
      for (UnsignedLong i = 0; i < dimension_; i++)
        {
          marginals[i] = getMarginal(i);
        }
      // The n-D quantile is defined as X(tau) = (F1^{-1}(tau), ..., Fn^{-1}(tau)),
      // with tau such as C(X(tau)) = prob.
      // We solve this last equation with respect to tau using the secant method
      // Lower end of the bracketing interval
      NumericalScalar leftTau(prob);
      NumericalPoint left(dimension_);
      NumericalScalar leftCDF(0.0);
      // Upper end of the bracketing interval
      NumericalScalar rightTau(1.0 - (1.0 - prob) / dimension_);
      NumericalPoint right(dimension_);
      NumericalScalar rightCDF(1.0);
      // Main loop of the bisection method
      while ((rightCDF - leftCDF > 2.0 * cdfEpsilon_) && (rightTau - leftTau > 2.0 * cdfEpsilon_))
        {
          NumericalScalar middleTau(0.5 * (rightTau + leftTau));
          NumericalPoint middle(dimension_);
          for (UnsignedLong i = 0; i < dimension_; i++)
            {
            middle[i] = marginals[i].computeQuantile(middleTau)[0];
            }
          NumericalScalar middleCDF(computeCDF(middle));
          if (middleCDF > prob)
            {
            right = middle;
            rightTau = middleTau;
            rightCDF = middleCDF;
            }
          else
            {
            left = middle;
            leftTau = middleTau;
            leftCDF = middleCDF;
            }
        }
      return 0.5 * (left + right);
      }

      /* Get the mathematical and numerical range of the distribution.
       Its mathematical range is the smallest closed interval outside
       of which the PDF is zero, and the numerical range is the interval
       outside of which the PDF is rounded to zero in double precision */
00635       DistributionImplementation::Interval DistributionImplementation::getRange() const
      {
      return range_;
      }

      void DistributionImplementation::setRange(const Interval & range)
      {
      if (range.getDimension() != dimension_) throw InvalidArgumentException(HERE) << "Error: the given range has a dimension incompatible with the dimension of the distribution.";
      range_ = range;
      }

      /* Compute the numerical range of the distribution given the parameters values */
00647       void DistributionImplementation::computeRange()
      {
      throw NotYetImplementedException(HERE);
      }

      /* Compute the mean of the distribution
         based on the relation:
       mean = Int_0^1 quantile(u) du
       Integration using the sinh tanh method, see:
       David H. Bailey, Jonathan M. Borwein, "Effective Error Bounds in Euler-Maclaurin-Based Quadrature Schemes", working paper, 21 June 2005.
      */
00658       void DistributionImplementation::computeMean() const throw(NotDefinedException)
      {
      UnsignedLong dimension(getDimension());
      mean_ = NumericalPoint(dimension);
      NumericalScalar epsilon(sqrt(DefaultQuantileEpsilon));
      UnsignedLong MaximumLevel(DefaultLevelNumber);
      // For each component
      for(UnsignedLong component = 0; component < dimension; ++component)
        {
          NumericalScalar h(0.5);
          UnsignedLong N(6);
          Distribution marginalDistribution(getMarginal(component));
          // Central term
          mean_[component] = h * 0.5 * marginalDistribution.computeQuantile(0.5)[0];
          // First block
          for (UnsignedLong j = 1; j <= N; ++j)
            {
            NumericalScalar hj(h * j);
            NumericalScalar expHj(exp(hj));
            NumericalScalar iexpHj(1.0 / expHj);
            NumericalScalar sinhHj(0.5 * (expHj - iexpHj));
            NumericalScalar expSinhHj(exp(sinhHj));
            NumericalScalar iexpSinhHj(1.0 / expSinhHj);
            NumericalScalar iTwoCoshSinhHj(1.0 / (expSinhHj + iexpSinhHj));
            NumericalScalar xjm(iexpSinhHj * iTwoCoshSinhHj);
            NumericalScalar xjp(expSinhHj * iTwoCoshSinhHj);
            NumericalScalar wj((expHj + iexpHj) * iTwoCoshSinhHj * iTwoCoshSinhHj);
            mean_[component] += h * wj * (marginalDistribution.computeQuantile(xjm)[0] + marginalDistribution.computeQuantile(xjp)[0]);
            } // End of first block
          // Sequential addition of half-blocks
          NumericalScalar error(1.0);
          UnsignedLong level(0);
          while( (error > epsilon) && (level < MaximumLevel))
            {
            ++level;
            h *= 0.5;
            mean_[component] *= 0.5;
            NumericalScalar delta(0.0);
            for (UnsignedLong j = 0; j <= N; ++j)
              {
                NumericalScalar hj(h * (2 * j + 1));
                NumericalScalar expHj(exp(hj));
                NumericalScalar iexpHj(1.0 / expHj);
                NumericalScalar sinhHj(0.5 * (expHj - iexpHj));
                NumericalScalar expSinhHj(exp(sinhHj));
                NumericalScalar iexpSinhHj(1.0 / expSinhHj);
                NumericalScalar iTwoCoshSinhHj(1.0 / (expSinhHj + iexpSinhHj));
                NumericalScalar xjm(iexpSinhHj * iTwoCoshSinhHj);
                NumericalScalar xjp(expSinhHj * iTwoCoshSinhHj);
                NumericalScalar wj((expHj + iexpHj) * iTwoCoshSinhHj * iTwoCoshSinhHj);
                delta += h * wj * (marginalDistribution.computeQuantile(xjm)[0] + marginalDistribution.computeQuantile(xjp)[0]);
              }
            error = fabs((delta - mean_[component]) / (1.0 + fabs(delta)));
            mean_[component] += delta;
            N *= 2;
            } // End of half-block
        } // End of each component
      isAlreadyComputedMean_ = true;
      }

      /* Get the mean of the distribution */
00719       DistributionImplementation::NumericalPoint DistributionImplementation::getMean() const throw(NotDefinedException)
      {
      if (!isAlreadyComputedMean_)
        {
          computeMean();
        }
      return mean_;
      }

      /* Compute the nodes and weights for a 1D gauss quadrature over [-1, 1] */
00729       void DistributionImplementation::computeGaussNodesAndWeights() const
      {
      int integrationNodesNumber(integrationNodesNumber_);
      gaussNodesAndWeights_ = NumericalSample(2, integrationNodesNumber);
      // First, build a symmetric tridiagonal matrix whose eigenvalues are the nodes of the
        // gauss integration rule
      char jobz('V');
      int ljobz(1);
      NumericalPoint d(integrationNodesNumber);
      NumericalPoint e(integrationNodesNumber);
      for (UnsignedLong i = 1; i < (UnsignedLong)integrationNodesNumber; ++i)
        {
          e[i - 1] = 0.5 / sqrt(1.0 - pow(2.0 * i, -2));
        }
      int ldz(integrationNodesNumber);
      SquareMatrix z(integrationNodesNumber);
      NumericalPoint work(2 * integrationNodesNumber - 2);
      int info;
      DSTEV_F77(&jobz, &integrationNodesNumber, &d[0], &e[0], &z(0, 0), &ldz, &work[0], &info, &ljobz);
      for (UnsignedLong i = 0; i < (UnsignedLong)integrationNodesNumber; ++i)
        {
          // Nodes
          gaussNodesAndWeights_[0][i] = d[i];
          // Weights
          gaussNodesAndWeights_[1][i] = 2.0 * pow(z(0, i), 2);
        }
      }

      /* Compute the covariance of the distribution */
00758       void DistributionImplementation::computeCovariance() const throw(NotDefinedException)
      {
      UnsignedLong dimension(getDimension());
      // We need this to initialize the covariance matrix in two cases:
      // + this is the first call to this routine (which could be checked by testing the dimension of the distribution and the dimension of the matrix
      // + the copula has changed from a non-independent one to the independent copula
      covariance_ = CovarianceMatrix(dimension);
      // First the diagonal terms, which are the marginal covariances
      // To ensure that the mean is up to date
      mean_ = getMean();
      // Get the standard deviation
      NumericalPoint standardDeviation(getStandardDeviation());
      for(UnsignedLong component = 0; component < dimension; ++component)
        {
          covariance_(component, component) = standardDeviation[component] * standardDeviation[component];
        }
      // Off-diagonal terms if the copula is not the independent copula
      if (!hasIndependentCopula())
        {
          NumericalScalar delta(2.0);
            Indices indices(2);
          int N(8 * 2 * 2 * 2 * 2 * 2);
          NumericalScalar h(0.5 / 2 / 2 / 2 / 2 / 2);
          for(UnsignedLong rowIndex = 0; rowIndex < dimension; ++rowIndex)
            {
            indices[0] = rowIndex;
            Distribution marginalI(getMarginal(rowIndex));
            NumericalScalar mi(marginalI.computeQuantile(0.5)[0]);
            NumericalScalar di(marginalI.computeQuantile(0.75)[0]-marginalI.computeQuantile(0.25)[0]);
            for(UnsignedLong columnIndex = rowIndex + 1; columnIndex < dimension; ++columnIndex)
              {
                indices[1] = columnIndex;
                Distribution marginalDistribution(getMarginal(indices));
                if (!marginalDistribution.getImplementation()->hasIndependentCopula())
                  {
                  Distribution marginalJ(getMarginal(columnIndex));
                  NumericalScalar mj(marginalJ.computeQuantile(0.5)[0]);
                  NumericalScalar dj(marginalJ.computeQuantile(0.75)[0]-marginalJ.computeQuantile(0.25)[0]);
                  NumericalPoint xij(2);
                  xij[0] = mi;
                  xij[1] = mj;
                  NumericalScalar covarianceIJ(0.0);
                  // Then we loop over the integration points
                  for(int rowNodeIndex = -N; rowNodeIndex < N + 1; ++rowNodeIndex)
                    {
                      NumericalScalar hi(h * rowNodeIndex);
                      NumericalScalar expHi(exp(hi));
                      NumericalScalar iexpHi(1.0 / expHi);
                      NumericalScalar sinhHi(0.5 * (expHi - iexpHi));
                      NumericalScalar expSinhHi(exp(sinhHi));
                      NumericalScalar iexpSinhHi(1.0 / expSinhHi);
                      NumericalScalar iTwoCoshSinhHi(1.0 / (expSinhHi + iexpSinhHi));
                      NumericalScalar xip(mi + expSinhHi * iTwoCoshSinhHi * di * delta);
                      NumericalScalar wi((expHi + iexpHi) * iTwoCoshSinhHi * iTwoCoshSinhHi);
                      NumericalScalar cdfip(marginalI.computeCDF(xip));
                      for(int columnNodeIndex = -N; columnNodeIndex < N + 1; ++columnNodeIndex)
                        {
                        NumericalScalar hj(h * columnNodeIndex);
                        NumericalScalar expHj(exp(hj));
                        NumericalScalar iexpHj(1.0 / expHj);
                        NumericalScalar sinhHj(0.5 * (expHj - iexpHj));
                        NumericalScalar expSinhHj(exp(sinhHj));
                        NumericalScalar iexpSinhHj(1.0 / expSinhHj);
                        NumericalScalar iTwoCoshSinhHj(1.0 / (expSinhHj + iexpSinhHj));
                        NumericalScalar xjp(mj + expSinhHj * iTwoCoshSinhHj * dj * delta);
                        NumericalScalar wj((expHj + iexpHj) * iTwoCoshSinhHj * iTwoCoshSinhHj);
                        NumericalScalar cdfjp(marginalJ.computeCDF(xjp));
                        NumericalPoint inpp(2);
                        inpp[0] = xip;
                        inpp[1] = xjp;
                        covarianceIJ += delta * delta * di * dj * h * h * wi * wj * (marginalDistribution.computeCDF(inpp) - cdfip * cdfjp);
                        } // loop over J integration nodes
                    } // loop over I integration nodes
                  covariance_(rowIndex, columnIndex) = covarianceIJ;
                  }
              } // loop over column indices
            } // loop over row indices
        } // if !hasIndependentCopula
      isAlreadyComputedCovariance_ = true;
      } // computeCovariance

      /* Get the covariance of the distribution */
00840       DistributionImplementation::CovarianceMatrix DistributionImplementation::getCovariance() const throw(NotDefinedException)
      {
      if (!isAlreadyComputedCovariance_)
        {
          computeCovariance();
        }
      return covariance_;
      }

      /* integrationNodesNumber accessors */
00850       UnsignedLong DistributionImplementation::getIntegrationNodesNumber() const
      {
      return integrationNodesNumber_;
      }

      void DistributionImplementation::setIntegrationNodesNumber(const UnsignedLong integrationNodesNumber)
      {
      if (integrationNodesNumber != integrationNodesNumber_)
        {
          isAlreadyComputedMean_ = false;
          isAlreadyComputedCovariance_ = false;
          isAlreadyComputedGaussNodesAndWeights_ = false;
          integrationNodesNumber_ = integrationNodesNumber;
        }
      }

      /* Gauss nodes and weights accessor */
00867       DistributionImplementation::NumericalSample DistributionImplementation::getGaussNodesAndWeights() const
      {
      if (!isAlreadyComputedGaussNodesAndWeights_)
        {
          computeGaussNodesAndWeights();
        }
      return gaussNodesAndWeights_;
      }

      /* Get the standard deviation of the distribution */
00877       DistributionImplementation::NumericalPoint DistributionImplementation::getStandardDeviation() const throw(NotDefinedException)
      {
      UnsignedLong dimension(getDimension());
      NumericalScalar epsilon(sqrt(DefaultQuantileEpsilon));
      UnsignedLong MaximumLevel(DefaultLevelNumber + 1);
      mean_ = getMean();
      NumericalPoint standardDeviation(dimension);
      // For each component
      for(UnsignedLong component = 0; component < dimension; ++component)
        {
          NumericalScalar h(0.5);
          UnsignedLong N(6);
          Distribution marginalDistribution(getMarginal(component));
          // Central term
          standardDeviation[component] = h * 0.5 * pow(marginalDistribution.computeQuantile(0.5)[0] - mean_[component], 2.0);
          // First block
          for (UnsignedLong j = 1; j <= N; ++j)
            {
            NumericalScalar hj(h * j);
            NumericalScalar expHj(exp(hj));
            NumericalScalar iexpHj(1.0 / expHj);
            NumericalScalar sinhHj(0.5 * (expHj - iexpHj));
            NumericalScalar expSinhHj(exp(sinhHj));
            NumericalScalar iexpSinhHj(1.0 / expSinhHj);
            NumericalScalar iTwoCoshSinhHj(1.0 / (expSinhHj + iexpSinhHj));
            NumericalScalar xjm(iexpSinhHj * iTwoCoshSinhHj);
            NumericalScalar xjp(expSinhHj * iTwoCoshSinhHj);
            NumericalScalar wj((expHj + iexpHj) * iTwoCoshSinhHj * iTwoCoshSinhHj);
            standardDeviation[component] += h * wj * (pow(marginalDistribution.computeQuantile(xjm)[0] - mean_[component], 2.0) + pow(marginalDistribution.computeQuantile(xjp)[0] - mean_[component], 2.0));
            } // End of first block
          // Sequential addition of half-blocks
          NumericalScalar error(1.0);
          UnsignedLong level(0);
          while( (error > epsilon) && (level < MaximumLevel))
            {
            ++level;
            h *= 0.5;
            standardDeviation[component] *= 0.5;
            NumericalScalar delta(0.0);
            for (UnsignedLong j = 0; j <= N; ++j)
              {
                NumericalScalar hj(h * (2 * j + 1));
                NumericalScalar expHj(exp(hj));
                NumericalScalar iexpHj(1.0 / expHj);
                NumericalScalar sinhHj(0.5 * (expHj - iexpHj));
                NumericalScalar expSinhHj(exp(sinhHj));
                NumericalScalar iexpSinhHj(1.0 / expSinhHj);
                NumericalScalar iTwoCoshSinhHj(1.0 / (expSinhHj + iexpSinhHj));
                NumericalScalar xjm(iexpSinhHj * iTwoCoshSinhHj);
                NumericalScalar xjp(expSinhHj * iTwoCoshSinhHj);
                NumericalScalar wj((expHj + iexpHj) * iTwoCoshSinhHj * iTwoCoshSinhHj);
                delta += h * wj * (pow(marginalDistribution.computeQuantile(xjm)[0] - mean_[component], 2.0) + pow(marginalDistribution.computeQuantile(xjp)[0] - mean_[component], 2.0));
              }
            error = fabs((delta - standardDeviation[component]) / (1.0 + fabs(delta)));
            standardDeviation[component] += delta;
            N *= 2;
            } // End of half-block
          standardDeviation[component] = sqrt(standardDeviation[component]);
        } // End of each component
      return standardDeviation;
      }

      /* Get the skewness of the distribution */
00940       DistributionImplementation::NumericalPoint DistributionImplementation::getSkewness() const throw(NotDefinedException)
      {
      UnsignedLong dimension(getDimension());
      NumericalScalar epsilon(sqrt(DefaultQuantileEpsilon));
      UnsignedLong MaximumLevel(DefaultLevelNumber + 2);
      mean_ = getMean();
      NumericalPoint standardDeviation(getStandardDeviation());
      NumericalPoint skewness(dimension);
      // For each component
      for(UnsignedLong component = 0; component < dimension; ++component)
        {
          NumericalScalar h(0.5);
          UnsignedLong N(6);
          Distribution marginalDistribution(getMarginal(component));
          // Central term
          skewness[component] = h * 0.5 * pow(marginalDistribution.computeQuantile(0.5)[0] - mean_[component], 3.0);
          // First block
          for (UnsignedLong j = 1; j <= N; ++j)
            {
            NumericalScalar hj(h * j);
            NumericalScalar expHj(exp(hj));
            NumericalScalar iexpHj(1.0 / expHj);
            NumericalScalar sinhHj(0.5 * (expHj - iexpHj));
            NumericalScalar expSinhHj(exp(sinhHj));
            NumericalScalar iexpSinhHj(1.0 / expSinhHj);
            NumericalScalar iTwoCoshSinhHj(1.0 / (expSinhHj + iexpSinhHj));
            NumericalScalar xjm(iexpSinhHj * iTwoCoshSinhHj);
            NumericalScalar xjp(expSinhHj * iTwoCoshSinhHj);
            NumericalScalar wj((expHj + iexpHj) * iTwoCoshSinhHj * iTwoCoshSinhHj);
            skewness[component] += h * wj * (pow(marginalDistribution.computeQuantile(xjm)[0] - mean_[component], 3.0) + pow(marginalDistribution.computeQuantile(xjp)[0] - mean_[component], 3.0));
            } // End of first block
          // Sequential addition of half-blocks
          NumericalScalar error(1.0);
          UnsignedLong level(0);
          while( (error > epsilon) && (level < MaximumLevel))
            {
            ++level;
            h *= 0.5;
            skewness[component] *= 0.5;
            NumericalScalar delta(0.0);
            for (UnsignedLong j = 0; j <= N; ++j)
              {
                NumericalScalar hj(h * (2 * j + 1));
                NumericalScalar expHj(exp(hj));
                NumericalScalar iexpHj(1.0 / expHj);
                NumericalScalar sinhHj(0.5 * (expHj - iexpHj));
                NumericalScalar expSinhHj(exp(sinhHj));
                NumericalScalar iexpSinhHj(1.0 / expSinhHj);
                NumericalScalar iTwoCoshSinhHj(1.0 / (expSinhHj + iexpSinhHj));
                NumericalScalar xjm(iexpSinhHj * iTwoCoshSinhHj);
                NumericalScalar xjp(expSinhHj * iTwoCoshSinhHj);
                NumericalScalar wj((expHj + iexpHj) * iTwoCoshSinhHj * iTwoCoshSinhHj);
                delta += h * wj * (pow(marginalDistribution.computeQuantile(xjm)[0] - mean_[component], 3.0) + pow(marginalDistribution.computeQuantile(xjp)[0] - mean_[component], 3.0));
              }
            error = fabs((delta - skewness[component]) / (1.0 + fabs(delta)));
            skewness[component] += delta;
            N *= 2;
            } // End of half-block
          skewness[component] = skewness[component] / pow(standardDeviation[component], 3.0);
        } // End of each component
      return skewness;
      }

      /* Get the kurtosis of the distribution */
01004       DistributionImplementation::NumericalPoint DistributionImplementation::getKurtosis() const throw(NotDefinedException)
      {
      UnsignedLong dimension(getDimension());
      NumericalScalar epsilon(sqrt(DefaultQuantileEpsilon));
      UnsignedLong MaximumLevel(DefaultLevelNumber + 3);
      mean_ = getMean();
      NumericalPoint standardDeviation(getStandardDeviation());
      NumericalPoint kurtosis(dimension);
      //    NumericalPoint values(MaximumLevel + 1);
      // For each component
      for(UnsignedLong component = 0; component < dimension; ++component)
        {
          NumericalScalar h(0.5);
          UnsignedLong N(6);
          Distribution marginalDistribution(getMarginal(component));
          // Central term
          kurtosis[component] = h * 0.5 * pow(marginalDistribution.computeQuantile(0.5)[0] - mean_[component], 4.0);
          // First block
          for (UnsignedLong j = 1; j <= N; ++j)
            {
            NumericalScalar hj(h * j);
            NumericalScalar expHj(exp(hj));
            NumericalScalar iexpHj(1.0 / expHj);
            NumericalScalar sinhHj(0.5 * (expHj - iexpHj));
            NumericalScalar expSinhHj(exp(sinhHj));
            NumericalScalar iexpSinhHj(1.0 / expSinhHj);
            NumericalScalar iTwoCoshSinhHj(1.0 / (expSinhHj + iexpSinhHj));
            NumericalScalar xjm(iexpSinhHj * iTwoCoshSinhHj);
            NumericalScalar xjp(expSinhHj * iTwoCoshSinhHj);
            NumericalScalar wj((expHj + iexpHj) * iTwoCoshSinhHj * iTwoCoshSinhHj);
            kurtosis[component] += h * wj * (pow(marginalDistribution.computeQuantile(xjm)[0] - mean_[component], 4.0) + pow(marginalDistribution.computeQuantile(xjp)[0] - mean_[component], 4.0));
            } // End of first block
          //values[0] = kurtosis[component];
          // Sequential addition of half-blocks
          NumericalScalar error(1.0);
          UnsignedLong level(0);
          while( (error > epsilon) && (level < MaximumLevel))
            {
            ++level;
            h *= 0.5;
            kurtosis[component] *= 0.5;
            NumericalScalar delta(0.0);
            for (UnsignedLong j = 0; j <= N; ++j)
              {
                NumericalScalar hj(h * (2 * j + 1));
                NumericalScalar expHj(exp(hj));
                NumericalScalar iexpHj(1.0 / expHj);
                NumericalScalar sinhHj(0.5 * (expHj - iexpHj));
                NumericalScalar expSinhHj(exp(sinhHj));
                NumericalScalar iexpSinhHj(1.0 / expSinhHj);
                NumericalScalar iTwoCoshSinhHj(1.0 / (expSinhHj + iexpSinhHj));
                NumericalScalar xjm(iexpSinhHj * iTwoCoshSinhHj);
                NumericalScalar xjp(expSinhHj * iTwoCoshSinhHj);
                NumericalScalar wj((expHj + iexpHj) * iTwoCoshSinhHj * iTwoCoshSinhHj);
                delta += h * wj * (pow(marginalDistribution.computeQuantile(xjm)[0] - mean_[component], 4.0) + pow(marginalDistribution.computeQuantile(xjp)[0] - mean_[component], 4.0));
              }
            error = fabs((delta - kurtosis[component]) / (1.0 + fabs(delta)));
            kurtosis[component] += delta;
            N *= 2;
            } // End of half-block
          kurtosis[component] = kurtosis[component] / pow(standardDeviation[component], 4.0);
        } // End of each component
      return kurtosis;
      }

      /* Check if the distribution is elliptical */
01070       Bool DistributionImplementation::isElliptical() const
      {
      return false;
      }

      /* Check if the distribution is continuos */
01076       Bool DistributionImplementation::isContinuous() const
      {
      return true;
      }

      /* Tell if the distribution has elliptical copula */
01082       Bool DistributionImplementation::hasEllipticalCopula() const
      {
      return true;
      }

      /* Tell if the distribution has independent copula */
01088       Bool DistributionImplementation::hasIndependentCopula() const
      {
      return true;
      }

      /* Compute the density generator of the elliptical generator, i.e.
       *  the function phi such that the density of the distribution can
       *  be written as p(x) = phi(t(x-mu)R(x-mu))
       */
01097       NumericalScalar DistributionImplementation::computeDensityGenerator(const NumericalScalar betaSquare) const
      {
      throw NotYetImplementedException(HERE) << "in DistributionImplementation::computeDensityGenerator(const NumericalScalar betaSquare) const";
      }

      /* Compute the derivative of the density generator */
01103       NumericalScalar DistributionImplementation::computeDensityGeneratorDerivative(const NumericalScalar betaSquare) const
      {
      throw NotYetImplementedException(HERE) << "in DistributionImplementation::computeDensityGeneratorDerivative(const NumericalScalar betaSquare) const";
      }

      /* Compute the seconde derivative of the density generator */
01109       NumericalScalar DistributionImplementation::computeDensityGeneratorSecondDerivative(const NumericalScalar betaSquare) const
      {
      throw NotYetImplementedException(HERE) << "in DistributionImplementation::computeDensityGeneratorSecondDerivative(const NumericalScalar betaSquare) const";
      }

      /* Get the i-th marginal distribution */
01115       DistributionImplementation::Implementation DistributionImplementation::getMarginal(const UnsignedLong i) const throw(InvalidArgumentException)
      {
      if ((i != 0) || (dimension_ != 1))
        {
          throw NotYetImplementedException(HERE) << "in DistributionImplementation::getMarginal(const UnsignedLong i) const throw(InvalidArgumentException)";
        }
      return clone();
      }

      /* Get the distribution of the marginal distribution corresponding to indices dimensions */
01125       DistributionImplementation::Implementation DistributionImplementation::getMarginal(const Indices & indices) const throw(InvalidArgumentException)
      {
      if ((indices.getSize() != 1) || (indices[0] != 0) || (dimension_ != 1))
        {
          throw NotYetImplementedException(HERE) << "in DistributionImplementation::getMarginal(const Indices & indices) const throw(InvalidArgumentException)";
        }
      return clone();
      }

      /* Get the copula of a distribution */
01135       DistributionImplementation::Implementation DistributionImplementation::getCopula() const
      {
      if (dimension_ == 1) return new IndependentCopula(1);
      return new SklarCopula(*this);
      }

      /* Get the isoprobabilist transformation */
01142       DistributionImplementation::IsoProbabilisticTransformation DistributionImplementation::getIsoProbabilisticTransformation() const
      {
      // Special case for dimension 1
      if (dimension_ == 1)
        {
          DistributionCollection collection(1);
          collection[0] = *this;
          // First, build the isoprobabilistic transformation of the independent copula of dimension 1
          IsoProbabilisticTransformation copulaTransformation;
          copulaTransformation.setEvaluationImplementation(new NatafIndependentCopulaEvaluation(1));
          copulaTransformation.setGradientImplementation(new NatafIndependentCopulaGradient(1));
          copulaTransformation.setHessianImplementation(new NatafIndependentCopulaHessian(1));
            // Second, build the marginal transformation
          IsoProbabilisticTransformation marginalTransformation;
          marginalTransformation.setEvaluationImplementation(new MarginalTransformationEvaluation(collection));
          marginalTransformation.setGradientImplementation(new MarginalTransformationGradient(collection));
          marginalTransformation.setHessianImplementation(new MarginalTransformationHessian(collection));
          marginalTransformation.setParameters(getParametersCollection()[0]);
          IsoProbabilisticTransformation transformation(copulaTransformation, marginalTransformation);
            NumericalPointWithDescription parameters(getParametersCollection()[0]);
          UnsignedLong parametersDimension(parameters.getDimension());
          Description parametersDescription(parameters.getDescription());
          String name(parameters.getName());
          for (UnsignedLong i = 0; i < parametersDimension; i++)
            {
            parametersDescription[i] = OSS() << name << "_" << parametersDescription[i];
            }
          parameters.setDescription(parametersDescription);
          transformation.setParameters(parameters);
          return transformation;
        }
      // General case, Rosenblatt transformation
      return NumericalMathFunctionImplementation(new RosenblattEvaluation(clone()));
      }

      /* Get the inverse isoprobabilist transformation */
01178       DistributionImplementation::InverseIsoProbabilisticTransformation DistributionImplementation::getInverseIsoProbabilisticTransformation() const
      {
      // Special case for dimension 1
      if (dimension_ == 1)
        {
          DistributionCollection collection(1);
          collection[0] = *this;
          // First, build the inverse isoprobabilistic transformation of the independent copula of dimension 1
          InverseIsoProbabilisticTransformation inverseCopulaTransformation;
          inverseCopulaTransformation.setEvaluationImplementation(new InverseNatafIndependentCopulaEvaluation(1));
          inverseCopulaTransformation.setGradientImplementation(new InverseNatafIndependentCopulaGradient(1));
          inverseCopulaTransformation.setHessianImplementation(new InverseNatafIndependentCopulaHessian(1));
          // Second, build the inverse marginal transformation
          InverseIsoProbabilisticTransformation inverseMarginalTransformation;
          inverseMarginalTransformation.setEvaluationImplementation(new InverseMarginalTransformationEvaluation(collection));
          inverseMarginalTransformation.setGradientImplementation(new InverseMarginalTransformationGradient(collection));
          inverseMarginalTransformation.setHessianImplementation(new InverseMarginalTransformationHessian(collection));
          inverseMarginalTransformation.setParameters(getParametersCollection()[0]);
          InverseIsoProbabilisticTransformation inverseTransformation(inverseMarginalTransformation, inverseCopulaTransformation);
            NumericalPointWithDescription parameters(getParametersCollection()[0]);
          UnsignedLong parametersDimension(parameters.getDimension());
          Description parametersDescription(parameters.getDescription());
          String name(parameters.getName());
          for (UnsignedLong i = 0; i < parametersDimension; i++)
            {
            parametersDescription[i] = OSS() << name << "_" << parametersDescription[i];
            }
          parameters.setDescription(parametersDescription);
          inverseTransformation.setParameters(parameters);
          return inverseTransformation;
        }
      // General case, inverse Rosenblatt transformation
      return NumericalMathFunctionImplementation(new InverseRosenblattEvaluation(clone()));
      }

      /* Get the standard distribution */
01214       DistributionImplementation::Implementation DistributionImplementation::getStandardDistribution() const
      {
      return new Normal(dimension_);
      throw NotYetImplementedException(HERE) << "in DistributionImplementation::getStandardDistribution() const";
      }

      /* Compute the radial distribution CDF */
01221       NumericalScalar DistributionImplementation::computeRadialDistributionCDF(const NumericalScalar radius) const
      {
      throw NotYetImplementedException(HERE) << "in DistributionImplementation::computeRadialDistributionCDF(const NumericalScalar radius) const";
      }

      /* Draw the PDF of the distribution when its dimension is 1 */
01227       DistributionImplementation::Graph DistributionImplementation::drawPDF(const NumericalScalar xMin, const NumericalScalar xMax, const UnsignedLong pointNumber) const
      throw(InvalidDimensionException, InvalidArgumentException)
      {
      if (dimension_ != 1) throw InvalidDimensionException(HERE) << "Error: can draw a PDF only if dimension equals 1, here dimension=" << dimension_;
      if (xMax <= xMin) throw InvalidArgumentException(HERE) << "Error: cannot draw a PDF with xMax >= xMin, here xmin=" << xMin << " and xmax=" << xMax;
      if (pointNumber < 2) throw InvalidArgumentException(HERE) << "Error: cannot draw a PDF with a point number < 2";
      // Discretization of the x axis
      String title(OSS() << getName() << " PDF");
      Curve curvePDF(computePDF(xMin, xMax, pointNumber), "red", "solid", 2, title);
      Graph graphPDF(title, "x", "PDF", true, "topright");
      graphPDF.addDrawable(curvePDF);
      return graphPDF;
      }

      /* Draw the PDF of the distribution when its dimension is 1 */
      DistributionImplementation::Graph DistributionImplementation::drawPDF(const UnsignedLong pointNumber) const
      throw(InvalidDimensionException, InvalidArgumentException)
      {
      NumericalScalar xMin(computeQuantile(DistributionImplementation::QMin)[0]);
      NumericalScalar xMax(computeQuantile(DistributionImplementation::QMax)[0]);
      NumericalScalar delta(2.0 * (xMax - xMin) * (1.0 - 0.5 * (DistributionImplementation::QMax - DistributionImplementation::QMin)));
      return drawPDF(xMin - delta, xMax + delta, pointNumber);
      }

      /* Draw the PDF of a 1D marginal */
01252       DistributionImplementation::Graph DistributionImplementation::drawMarginal1DPDF(const UnsignedLong marginalIndex, const NumericalScalar xMin, const NumericalScalar xMax, const UnsignedLong pointNumber) const
      throw(InvalidArgumentException)
      {
      Graph marginalGraph(getMarginal(marginalIndex)->drawPDF(xMin, xMax, pointNumber));
      marginalGraph.setTitle(OSS() << getName() << ", " << description_[marginalIndex] << " component PDF");
      return marginalGraph;
      }

      /* Draw the PDF of the distribution when its dimension is 2 */
01261       DistributionImplementation::Graph DistributionImplementation::drawPDF(const NumericalPoint & xMin, const NumericalPoint & xMax, const NumericalPoint & pointNumber) const
      throw(InvalidDimensionException, InvalidArgumentException)
      {
      if ((xMin.getDimension() != 2) || (xMax.getDimension() != 2) || (pointNumber.getDimension() != 2)) throw InvalidArgumentException(HERE) << "Error: xMin, xMax and PointNumber must be 2D NumericalPoints";
      if ((pointNumber[0] <= 2) || (pointNumber[1] <= 2)) throw InvalidArgumentException(HERE) << "Error: the discretization must have at least 2 points per component";
      NumericalPoint discretization(2);
      NumericalPoint scaling(2);
      NumericalPoint origin(2);
      NumericalScalar nX(pointNumber[0] - 2);
      discretization[0] = nX;
      // Discretization of the first component
      NumericalSample x = Box(NumericalPoint(1, nX)).generate();
      origin[0] = 0.5 * (xMax[0] + xMin[0]);
      scaling[0] = xMax[0] - xMin[0];
      x.scale(NumericalPoint(1, scaling[0]));
      x.translate(NumericalPoint(1, origin[0]));
      NumericalScalar nY(pointNumber[1] - 2);
      discretization[1] = nY;
      // Discretization of the second component
      NumericalSample y = Box(NumericalPoint(1, nY)).generate();
      origin[1] = 0.5 * (xMax[1] + xMin[1]);
      scaling[1] = xMax[1] - xMin[1];
      y.scale(NumericalPoint(1, scaling[1]));
      y.translate(NumericalPoint(1, origin[1]));
      // Discretization of the XY plane
      NumericalSample xy = Box(discretization).generate();
      xy.scale(scaling);
      xy.translate(origin);
      NumericalSample z(xy.getSize(), 1);
      for (UnsignedLong indexZ = 0; indexZ < z.getSize(); ++indexZ)
        {
          z[indexZ][0] = computePDF(xy[indexZ]);
        }
      String title(OSS() << getName() << " iso-PDF");
      Contour isoPDF(Contour(x, y, z, NumericalPoint(0), Description(0), true, title));
      isoPDF.buildDefaultLevels();
      isoPDF.buildDefaultLabels();
      Graph graph(title, description_[0], description_[1], true, "topright");
      graph.addDrawable(isoPDF);
      return graph;
      }

      /* Draw the PDF of the distribution when its dimension is 2 */
      DistributionImplementation::Graph DistributionImplementation::drawPDF(const NumericalPoint & xMin, const NumericalPoint & xMax) const
      throw(InvalidDimensionException, InvalidArgumentException)
      {
      return drawPDF(xMin, xMax, NumericalPoint(2, DefaultPointNumber));
      }

      /* Draw the PDF of the distribution when its dimension is 2 */
      DistributionImplementation::Graph DistributionImplementation::drawPDF(const NumericalPoint & pointNumber) const
      throw(InvalidDimensionException, InvalidArgumentException)
      {
      NumericalPoint xMin(2);
      xMin[0] = (getMarginal(0)->computeQuantile(DistributionImplementation::QMin)[0]);
      xMin[1] = (getMarginal(1)->computeQuantile(DistributionImplementation::QMin)[0]);
      NumericalPoint xMax(2);
      xMax[0] = (getMarginal(0)->computeQuantile(DistributionImplementation::QMax)[0]);
      xMax[1] = (getMarginal(1)->computeQuantile(DistributionImplementation::QMax)[0]);
      NumericalPoint delta(2.0 * (xMax - xMin) * (1.0 - 0.5 * (DistributionImplementation::QMax - DistributionImplementation::QMin)));
      return drawPDF(xMin - delta, xMax + delta, pointNumber);
      }

      /* Draw the PDF of a 2D marginal */
01325       DistributionImplementation::Graph DistributionImplementation::drawMarginal2DPDF(const UnsignedLong firstMarginal, const UnsignedLong secondMarginal, const NumericalPoint & xMin, const NumericalPoint & xMax, const NumericalPoint & pointNumber) const
      throw(InvalidArgumentException)
      {
      Indices indices(2);
      indices[0] = firstMarginal;
      indices[1] = secondMarginal;
      Graph marginalGraph(getMarginal(indices)->drawPDF(xMin, xMax, pointNumber));
      marginalGraph.setTitle(OSS() << getName() << " (" << description_[firstMarginal] << ", " << description_[secondMarginal] << ") components iso-PDF");
      return marginalGraph;
      }

      /* Draw the PDF of the distribution when its dimension is 1 or 2 */
01337       DistributionImplementation::Graph DistributionImplementation::drawPDF() const
      throw(InvalidDimensionException, InvalidArgumentException)
      {
      UnsignedLong dimension(getDimension());
      // Generic interface for the 1D and 2D cases
      if (dimension == 1)
        {
          return drawPDF(DefaultPointNumber);
        }
      if (dimension == 2)
        {
          return drawPDF(NumericalPoint(2, DefaultPointNumber));
        }
      throw InvalidDimensionException(HERE) << "Error: can draw a PDF only if dimension equals 1 or 2, here dimension=" << dimension;
      }

      /* Draw the CDF of the distribution when its dimension is 1 */
01354       DistributionImplementation::Graph DistributionImplementation::drawCDF(const NumericalScalar xMin, const NumericalScalar xMax, const UnsignedLong pointNumber) const
      throw(InvalidDimensionException, InvalidArgumentException)
      {
      if (dimension_ != 1) throw InvalidDimensionException(HERE) << "Error: can draw a CDF only if dimension equals 1, here dimension=" << dimension_;
      if (xMax <= xMin) throw InvalidArgumentException(HERE) << "Error: cannot draw a CDF with xMax >= xMin, here xmin=" << xMin << " and xmax=" << xMax;
      if (pointNumber < 2) throw InvalidArgumentException(HERE) << "Error: cannot draw a CDF with a point number < 2";
      String title(OSS() << getName() << " CDF");
      Curve curveCDF(computeCDF(xMin, xMax, pointNumber), "red", "solid", 2, title);
      Graph graphCDF(title, "x", "CDF", true, "topright");
      graphCDF.addDrawable(curveCDF);
      return graphCDF;
      }

      /* Draw the CDF of the distribution when its dimension is 1 */
      DistributionImplementation::Graph DistributionImplementation::drawCDF(const UnsignedLong pointNumber) const
      throw(InvalidDimensionException, InvalidArgumentException)
      {
      NumericalScalar xMin(computeQuantile(DistributionImplementation::QMin)[0]);
      NumericalScalar xMax(computeQuantile(DistributionImplementation::QMax)[0]);
      NumericalScalar delta(2.0 * (xMax - xMin) * (1.0 - 0.5 * (DistributionImplementation::QMax - DistributionImplementation::QMin)));
      return drawCDF(xMin - delta, xMax + delta, pointNumber);
      }

      /* Draw the CDF of a 1D marginal */
01378       DistributionImplementation::Graph DistributionImplementation::drawMarginal1DCDF(const UnsignedLong marginalIndex, const NumericalScalar xMin, const NumericalScalar xMax, const UnsignedLong pointNumber) const
      throw(InvalidArgumentException)
      {
      Graph marginalGraph(getMarginal(marginalIndex)->drawCDF(xMin, xMax, pointNumber));
      marginalGraph.setTitle(OSS() << getName() << ", " << description_[marginalIndex] << " component CDF");
      return marginalGraph;
      }

      /* Draw the CDF of the distribution when its dimension is 2 */
01387       DistributionImplementation::Graph DistributionImplementation::drawCDF(const NumericalPoint & xMin, const NumericalPoint & xMax, const NumericalPoint & pointNumber) const
      throw(InvalidDimensionException, InvalidArgumentException)
      {
      if ((xMin.getDimension() != 2) || (xMax.getDimension() != 2) || (pointNumber.getDimension() != 2)) throw InvalidArgumentException(HERE) << "Error: xMin, xMax and PointNumber must be 2D NumericalPoints";
      if ((pointNumber[0] <= 2) || (pointNumber[1] <= 2)) throw InvalidArgumentException(HERE) << "Error: the discretization must have at least 2 points per component";
      NumericalPoint discretization(2);
      NumericalPoint scaling(2);
      NumericalPoint origin(2);
      NumericalScalar nX(pointNumber[0] - 2);
      discretization[0] = nX;
      // Discretization of the first component
      NumericalSample x = Box(NumericalPoint(1, nX)).generate();
      origin[0] = 0.5 * (xMax[0] + xMin[0]);
      scaling[0] = xMax[0] - xMin[0];
      x.scale(NumericalPoint(1, scaling[0]));
      x.translate(NumericalPoint(1, origin[0]));
      NumericalScalar nY(pointNumber[1] - 2);
      discretization[1] = nY;
      // Discretization of the second component
      NumericalSample y = Box(NumericalPoint(1, nY)).generate();
      origin[1] = 0.5 * (xMax[1] + xMin[1]);
      scaling[1] = xMax[1] - xMin[1];
      y.scale(NumericalPoint(1, scaling[1]));
      y.translate(NumericalPoint(1, origin[1]));
      // Discretization of the XY plane
      NumericalSample xy = Box(discretization).generate();
      xy.scale(scaling);
      xy.translate(origin);
      NumericalSample z(xy.getSize(), 1);
      for (UnsignedLong indexZ = 0; indexZ < z.getSize(); ++indexZ)
        {
          z[indexZ][0] = computeCDF(xy[indexZ]);
        }
      String title(OSS() << "iso-CDF of " << getName());
      Contour isoCDF(Contour(x, y, z, NumericalPoint(0), Description(0), true, title));
      isoCDF.buildDefaultLevels();
      isoCDF.buildDefaultLabels();
      Graph graph(title, description_[0], description_[1], true, "topright");
      graph.addDrawable(isoCDF);
      return graph;
      }

      /* Draw the CDF of the distribution when its dimension is 2 */
      DistributionImplementation::Graph DistributionImplementation::drawCDF(const NumericalPoint & xMin, const NumericalPoint & xMax) const
      throw(InvalidDimensionException, InvalidArgumentException)
      {
      return drawCDF(xMin, xMax, NumericalPoint(2, DefaultPointNumber));
      }

      /* Draw the CDF of the distribution when its dimension is 2 */
      DistributionImplementation::Graph DistributionImplementation::drawCDF(const NumericalPoint & pointNumber) const
      throw(InvalidDimensionException, InvalidArgumentException)
      {
      NumericalPoint xMin(2);
      xMin[0] = (getMarginal(0)->computeQuantile(DistributionImplementation::QMin)[0]);
      xMin[1] = (getMarginal(1)->computeQuantile(DistributionImplementation::QMin)[0]);
      NumericalPoint xMax(2);
      xMax[0] = (getMarginal(0)->computeQuantile(DistributionImplementation::QMax)[0]);
      xMax[1] = (getMarginal(1)->computeQuantile(DistributionImplementation::QMax)[0]);
      NumericalPoint delta(2.0 * (xMax - xMin) * (1.0 - 0.5 * (DistributionImplementation::QMax - DistributionImplementation::QMin)));
      return drawCDF(xMin - delta, xMax + delta, pointNumber);
      }

      /* Draw the CDF of a 2D marginal */
01451       DistributionImplementation::Graph DistributionImplementation::drawMarginal2DCDF(const UnsignedLong firstMarginal, const UnsignedLong secondMarginal, const NumericalPoint & xMin, const NumericalPoint & xMax, const NumericalPoint & pointNumber) const
      throw(InvalidArgumentException)
      {
      Indices indices(2);
      indices[0] = firstMarginal;
      indices[1] = secondMarginal;
      Graph marginalGraph(getMarginal(indices)->drawCDF(xMin, xMax, pointNumber));
      marginalGraph.setTitle(OSS() << getName() << " (" << description_[firstMarginal] << ", " << description_[secondMarginal] << ") components iso-CDF");
      return marginalGraph;
      }

      /* Draw the CDF of the distribution when its dimension is 1 or 2 */
01463       DistributionImplementation::Graph DistributionImplementation::drawCDF() const
      throw(InvalidDimensionException, InvalidArgumentException)
      {
      UnsignedLong dimension(getDimension());
      // Generic interface for the 1D and 2D cases
      if (dimension == 1)
        {
          return drawCDF(DefaultPointNumber);
        }
      if (dimension == 2)
        {
          return drawCDF(NumericalPoint(2, DefaultPointNumber));
        }
      throw InvalidDimensionException(HERE) << "Error: can draw a CDF only if dimension equals 1 or 2, here dimension=" << dimension;
      }

      /* Parameters value and description accessor */
01480       DistributionImplementation::NumericalPointWithDescriptionCollection DistributionImplementation::getParametersCollection() const
      {
      return DistributionImplementation::NumericalPointWithDescriptionCollection(0);
      }

      /* Parameters number */
01486       UnsignedLong DistributionImplementation::getParametersNumber() const
      {
      NumericalPointWithDescriptionCollection parametersCollection(getParametersCollection());
      const UnsignedLong size(parametersCollection.getSize());
      UnsignedLong parametersNumber(0);
      for (UnsignedLong i = 0; i < size; ++i)
        {
          parametersNumber += parametersCollection[i].getDimension();
        }
      return parametersNumber;
      }

      /* Description accessor */
01499       void DistributionImplementation::setDescription(const Description & description)
      {
      const UnsignedLong size(description.getSize());
      if (size != getDimension()) throw InvalidArgumentException(HERE) << "Error: the description must have the same size than the distribution dimension, here size=" << size << " and dimension=" << getDimension();
      // Check if the description is valid
      // First, copy the description
      Description test(description);
      // Second, sort the copy
      std::sort(test.begin(), test.end());
      // Third, move the duplicates at the end
      Description::const_iterator it = std::unique(test.begin(), test.end());
      // Fourth, check if there was any duplicate
      if (it != test.end())
        {
          Log::Info(OSS() << "Warning! The description of the distribution " << getName() << " is " << description << " and cannot identify uniquely the marginal distribution. Appended unique identifier to fix it:");
          Description newDescription(description);
          for (UnsignedLong i = 0; i < size; ++i)
            {
            newDescription[i] = OSS() << "marginal_" << i + 1 << "_" << description[i];
            }
          Log::Info(OSS() << "the new description is " << newDescription);
          description_ = newDescription;
        }
      else
        {
          description_ = description;
        }
      }

      /* Description accessot */
      DistributionImplementation::Description DistributionImplementation::getDescription() const
      {
      return description_;
      }

      /* Accessor to PDF computation precision */
01535       NumericalScalar DistributionImplementation::getPDFEpsilon() const
      {
      return pdfEpsilon_;
      }

      /* Accessor to CDF computation precision */
01541       NumericalScalar DistributionImplementation::getCDFEpsilon() const
      {
      return cdfEpsilon_;
      }

      /* Get a positon indicator for a 1D distribution */
01547       NumericalScalar DistributionImplementation::getPositionIndicator() const
      {
      if (dimension_ != 1) throw InvalidDimensionException(HERE) << "Error: cannot get the position indicator of a distribution with dimension > 1";
      // Return the median of the distribution
      return computeQuantile(0.5)[0];
      }

      /* Get a dispersion indicator for a 1D distribution */
01555       NumericalScalar DistributionImplementation::getDispersionIndicator() const
      {
      if (dimension_ != 1) throw InvalidDimensionException(HERE) << "Error: cannot get the dispertion indicator of a distribution with dimension > 1";
      // Return the interquartile of the distribution
      return computeQuantile(0.75)[0] - computeQuantile(0.25)[0];
      }



      /* Method save() stores the object through the StorageManager */
01565       void DistributionImplementation::save(const StorageManager::Advocate & adv) const
      {
      PersistentObject::save(adv);
      adv.writeValue(mean_, StorageManager::MemberNameAttribute, "mean_");
      adv.writeValue(covariance_, StorageManager::MemberNameAttribute, "covariance_");
      adv.writeValue(gaussNodesAndWeights_, StorageManager::MemberNameAttribute, "gaussNodesAndWeights_");
      adv.writeValue("integrationNodesNumber_", integrationNodesNumber_);
      adv.writeValue("isAlreadyComputedMean_", isAlreadyComputedMean_);
      adv.writeValue("isAlreadyComputedCovariance_", isAlreadyComputedCovariance_);
      adv.writeValue("isAlreadyComputedGaussNodesAndWeights_", isAlreadyComputedGaussNodesAndWeights_);
      adv.writeValue("dimension_", dimension_);
      adv.writeValue("weight_", weight_);
      adv.writeValue(description_, StorageManager::MemberNameAttribute, "description_");
      }

      /* Method load() reloads the object from the StorageManager */
01581       void DistributionImplementation::load(const StorageManager::Advocate & adv)
      {
      PersistentObject::load(adv);
      adv.readValue(description_, StorageManager::MemberNameAttribute, "description_");
      adv.readValue(mean_, StorageManager::MemberNameAttribute, "mean_");
      adv.readValue(covariance_, StorageManager::MemberNameAttribute, "covariance_");
      adv.readValue(gaussNodesAndWeights_, StorageManager::MemberNameAttribute, "gaussNodesAndWeights_");
      String name;
      UnsignedLong unsignedLongValue;
      StorageManager::List objList = adv.getList(StorageManager::UnsignedLongEntity);
      for(objList.firstValueToRead(); objList.moreValuesToRead(); objList.nextValueToRead()) {
        if (objList.readValue(name, unsignedLongValue)) {
          if (name == "integrationNodesNumber_") integrationNodesNumber_ = unsignedLongValue;
          if (name == "dimension_") dimension_ = unsignedLongValue;
        }
      }
      Bool boolValue;
      objList = adv.getList(StorageManager::BoolEntity);
      for(objList.firstValueToRead(); objList.moreValuesToRead(); objList.nextValueToRead()) {
        if (objList.readValue(name, boolValue)) {
          if (name == "isAlreadyComputedMean_") isAlreadyComputedMean_ = boolValue;
          if (name == "isAlreadyComputedCovariance_") isAlreadyComputedCovariance_ = boolValue;
          if (name == "isAlreadyComputedGaussNodesAndWeights_") isAlreadyComputedCovariance_ = boolValue;
        }
      }
      NumericalScalar scalarValue;
      objList = adv.getList(StorageManager::NumericalScalarEntity);
      for(objList.firstValueToRead(); objList.moreValuesToRead(); objList.nextValueToRead()) {
        if (objList.readValue(name, scalarValue)) {
          if (name == "weight_") weight_ = scalarValue;
        }
      }
      }

    } /* namespace Model */
  } /* namespace Uncertainty */
} /* namespace OpenTURNS */

Generated by  Doxygen 1.6.0   Back to index