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

Normal.cxx

Go to the documentation of this file.
//                                               -*- C++ -*-
/**
 *  @file  Normal.cxx
 *  @brief The Normal distribution
 *
 *  (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: Normal.cxx 995 2008-10-31 10:52:04Z dutka $
 */
#include <cmath>
#include "Normal.hxx"
#include "SpecFunc.hxx"
#include "Mvndstpack.hxx"
#include "Tvpack.hxx"
#include "Log.hxx"
#include "OSS.hxx"
#include "DistFunc.hxx"
#include "PersistentObjectFactory.hxx"
#include "Matrix.hxx"
#include "MatrixImplementation.hxx"
#include "IdentityMatrix.hxx"
#include "NormalCopula.hxx"

namespace OpenTURNS {

  namespace Uncertainty {

    namespace Distribution {

      const UnsignedLong Normal::MaximumNumberOfPoints;
      const UnsignedLong Normal::MinimumNumberOfPoints;
      const NumericalScalar Normal::MaximumCDFEpsilon;
      const NumericalScalar Normal::MinimumCDFEpsilon;

      typedef Base::Type::IdentityMatrix               IdentityMatrix;
      typedef Base::Type::Matrix                       Matrix;
      typedef Base::Type::MatrixImplementation         MatrixImplementation;
      typedef Base::Common::Log                        Log;
      typedef Base::Common::NotYetImplementedException NotYetImplementedException;

00057       CLASSNAMEINIT(Normal);

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

      /* Constructor for multiD standard normal distribution */
      Normal::Normal(const UnsignedLong dimension)
      throw(InvalidArgumentException)
      : EllipticalDistribution(NumericalPoint(dimension, 0.0), NumericalPoint(dimension, 1.0), CorrelationMatrix(dimension), "Normal"),
        normalizationFactor_(1.0 / sqrt(pow(2 * M_PI, dimension))),
        hasIndependentCopula_(true)
      {
00068       // Nothing to do
      computeRange();
      }

      /* Constructor for 1D normal distribution */
      Normal::Normal(const NumericalScalar mu, const NumericalScalar sd)
      throw(InvalidArgumentException)
      : EllipticalDistribution(NumericalPoint(1, mu), NumericalPoint(1, sd), CorrelationMatrix(1), "Normal"),
        normalizationFactor_(1.0 / sqrt(2 * M_PI)),
        hasIndependentCopula_(true)
      {
00079       // Nothing to do
      computeRange();
      }

      /* Constructor for multiD normal distribution */
      Normal::Normal(const NumericalPoint & mean,
                 const NumericalPoint & sigma,
                 const CorrelationMatrix & R)
      throw(InvalidArgumentException)
      : EllipticalDistribution(mean, sigma, R, "Normal"),
        normalizationFactor_(1.0 / sqrt(pow(2 * M_PI, mean.getDimension()))),
        hasIndependentCopula_(false)
      {
      checkIndependentCopula();
      computeRange();
      }

      Normal::Normal(const NumericalPoint & mean,
                 const CovarianceMatrix & C)
      throw(InvalidArgumentException)
      : EllipticalDistribution(mean, NumericalPoint(mean.getDimension(), 1.0), IdentityMatrix(mean.getDimension()), "Normal"),
        normalizationFactor_(1.0 / sqrt(pow(2 * M_PI, mean.getDimension()))),
        hasIndependentCopula_(false)
      {
      UnsignedLong dimension(mean.getDimension());
      if (C.getDimension() != dimension) throw InvalidArgumentException(HERE) << "Error: the mean vector and the covariance matrix have incompatible dimensions";
      if (!C.isPositiveDefinite()) throw InvalidArgumentException(HERE) << "Error: the covariance matrix is not positive definite";
      NumericalPoint sigma(dimension);
      CorrelationMatrix R(dimension);
      for (UnsignedLong i = 0; i < dimension; ++i)
        {
          sigma[i] = sqrt(C(i, i));
          for (UnsignedLong j = 0; j < i; ++j)
            {
            R(i, j) = C(i, j) / (sigma[i] * sigma[j]);
            }
        }
      setSigma(sigma);
      setCorrelationMatrix(R);
00118       checkIndependentCopula();
      computeRange();
      }

      /* String converter */
      String Normal::str() const
      {
      OSS oss;
      oss << "class=" << Normal::GetClassName()
          << " name=" << getName()
          << " dimension=" << getDimension()
          << " mean={" << mean_
          << "} sigma={" << sigma_
          << "} correlationMatrix={" << R_
00132           << "}";
      return oss;
      }
  
      /* Virtual constructor */
      Normal * Normal::clone() const
00138       {
      return new Normal(*this);
      }

      /* Compute the numerical range of the distribution given the parameters values */
      void Normal::computeRange()
      {
      const UnsignedLong dimension(getDimension());
      NumericalPoint lowerBound(getMean());
      NumericalPoint upperBound(lowerBound);
      for (UnsignedLong i = 0; i < dimension; ++i)
        {
          const NumericalScalar delta(0.5 * sigma_[i] * sqrt(4.0 * (log(SpecFunc::ISQRT2PI / sigma_[i]) - SpecFunc::LogMinNumericalScalar)));
          lowerBound[i] -= delta;
          upperBound[i] += delta;
        }
      const Interval::BoolCollection finiteLowerBound(dimension, false);
00155       const Interval::BoolCollection finiteUpperBound(dimension, false);
      setRange(Interval(lowerBound, upperBound, finiteLowerBound, finiteUpperBound));
      }

      /* Get one realization of the distribution */
      Normal::NumericalPoint Normal::getRealization() const
      {
      UnsignedLong dimension(getDimension());
      if (dimension == 1)
        {
          return NumericalPoint(1, mean_[0] + sigma_[0] * DistFunc::rNormal());
        }
      NumericalPoint value(dimension);
      // First, a realization of independant standard coordinates
      for (UnsignedLong i = 0; i < dimension; i++)
        {
          value[i] = DistFunc::rNormal();
        }
      // Then, transform the independant standard coordinates into the needed ones */
      if (hasIndependentCopula_)
        {
          for (UnsignedLong i = 0; i < dimension; i++)
            {
            value[i] *= sigma_[i];
            value[i] += mean_[i];
            }
          return value;
        }
      // General case
      return cholesky_ * value + mean_;
00185       }
     
      /* Compute the density generator of the ellipticalal generator, i.e.
       *  the function phi such that the density of the distribution can
       *  be written as p(x) = phi(t(x-mu)S^(-1)(x-mu))                      */
      NumericalScalar Normal::computeDensityGenerator(const NumericalScalar betaSquare) const
00191       {
      return normalizationFactor_ * exp(-0.5 * betaSquare);
      }

      /* Compute the derivative of the density generator */
      NumericalScalar Normal::computeDensityGeneratorDerivative(const NumericalScalar betaSquare) const
00197       {
      return -0.5 * normalizationFactor_ * exp(-0.5 * betaSquare);
      }

      /* Compute the seconde derivative of the density generator */
      NumericalScalar Normal::computeDensityGeneratorSecondDerivative(const NumericalScalar betaSquare) const
      {
      return 0.25 * normalizationFactor_ * exp(-0.5 * betaSquare);
      }

      /* Get the CDF of the distribution */
      NumericalScalar Normal::computeCDF(const NumericalPoint & point, const Bool tail) const
      {
      const UnsignedLong dimension(getDimension());
      if (point.getDimension() != dimension) throw InvalidArgumentException(HERE) << "Error: the given point has a dimension incompatible with the distribution.";
      // Special case for dimension 1
      if (dimension == 1) return DistFunc::pNormal((point[0] - mean_[0]) / sigma_[0]);
      // Normalize the point to use the standard form of the multivariate normal distribution
      NumericalPoint u(normalize(point));
      /* Special treatment for independent components */
      if (hasIndependentCopula_)
        {
          NumericalScalar value(DistFunc::pNormal(u[0]));
          for (UnsignedLong i = 1; i < dimension; ++i)
            {
            value *= DistFunc::pNormal(u[i], tail);
            }
          return value;
        }
      /* General case */
      // For the bidimensional case, use specialized high precision routine
      if (dimension == 2)
        {
          int zero(0);
          double r(R_(0, 1));
          return BVTL_F77(&zero, &(u[0]), &(u[1]), &r);
        }
      // For the tridimensional case, use specialized high precision routine
      if (dimension == 3)
        {
          int zero(0);
          double r[3];
          r[0] = R_(1, 0);
          r[1] = R_(2, 0);
          r[2] = R_(2, 1);
          double eps(MaximumCDFEpsilon);
          cdfEpsilon_ = eps;
          return TVTL_F77(&zero, &(u[0]), &r[0], &eps);
        }
      // For larger dimension, use an adaptive integration
      if (dimension <= 500)
        {
          NumericalPoint lower(dimension, 0.0);
          std::vector<int> infin(dimension, 0);
          NumericalPoint correl(dimension * (dimension - 1) / 2, 0.0);
          /* Copy the correlation matrix in the proper format for mvndst */
          for (UnsignedLong i = 0; i < dimension; i++)
            {
            for (UnsignedLong j = 0; j < i; j++)
              {
                correl[j + i * (i - 1) / 2] = R_(j, i);
              }
            }
          int maxpts(MinimumNumberOfPoints);
          // Use only relative precision
          double abseps(0.0);
          // Reduce the precision according to the dimension. It ranges from MaximumCDFEpsilon for dimension=4 to MinimumCDFEpsilon for dimension=500
          double releps(MaximumCDFEpsilon * pow(MinimumCDFEpsilon / MaximumCDFEpsilon, (4.0 + dimension) / 496.0));
          double error;
          double value;
          int inform;
          int dim = static_cast<UnsignedLong>( dimension );
          do
            {
            MVNDST_F77(&dim, &lower[0], &u[0], &infin[0], &correl[0], &maxpts, &abseps, &releps, &error, &value, &inform);
            if (inform == 1)
              {
                Log::Warn(OSS() << "Warning, in Normal::computeCDF(), the required precision has not been achieved with maxpts=" << NumericalScalar(maxpts) << ", we only got an absolute error of " << error << " and a relative error of " << error / value << ". We retry with maxpoint=" << 10 * maxpts);
                maxpts *= 10;
              }
            } while ((static_cast<UnsignedLong>(maxpts) <= MaximumNumberOfPoints) && (inform == 1));
          if (inform == 1) Log::Warn(OSS() << "Warning, in Normal::computeCDF(), the required precision has not been achieved with maxpts=" << NumericalScalar(maxpts) << ", we only got an absolute error of " << error << " and a relative error of " << error / value << ". No more retry.");
          if ((value < 0.0) || (value > 1.0)) Log::Warn(OSS() << "Warning, in Normal::computeCDF(), got a value outside of [0, 1], value=" << value << " your dependence structure might be too complex. The value will be truncated.");
          cdfEpsilon_ = error;
          return std::min(std::max(value, 0.0), 1.0);
        }
      // For very large dimension, use a MonteCarlo algorithm
      Log::Warn(OSS() << "Warning, in Normal::computeCDF(), the dimension is very high. We will use a Monte Carlo method for the computation with a relative precision of 0.1% at 99% confidence level and a maximum of " << 10.0 * MaximumNumberOfPoints << " realizations. Expect a long running time and a poor accuracy for low values of the CDF...");
      NumericalScalar value(0.0);
      NumericalScalar variance(0.0);
      NumericalScalar a99(DistFunc::qNormal(0.995));
      UnsignedLong outerMax(10 * MaximumNumberOfPoints / MinimumNumberOfPoints);
      NumericalScalar precision(0.0);
      for (UnsignedLong indexOuter = 0; indexOuter < outerMax; ++indexOuter)
        {
          NumericalScalar valueBlock(0.0);
          NumericalScalar varianceBlock(0.0);
          for (UnsignedLong indexSample = 0; indexSample < MinimumNumberOfPoints; ++indexSample)
            {
            Bool inside(true);
            NumericalPoint realization(getRealization());
            // Check if the realization is in the integration domain
            for (UnsignedLong i = 0; i < dimension; ++i)
              {
                inside = realization[i] <= point[i];
                if (!inside) break;
              }
            // ind value is 1.0 if the realization is inside of the integration domain, 0.0 else.
            NumericalScalar ind(inside);
            NumericalScalar norm(1.0 / (indexSample + 1.0));
            varianceBlock = (varianceBlock * indexSample + (1.0 - norm) * (valueBlock - ind) * (valueBlock - ind)) * norm;
            valueBlock = (valueBlock * indexSample + ind) * norm;
            }
          NumericalScalar norm(1.0 / (indexOuter + 1.0));
          variance = (varianceBlock + indexOuter * variance + (1.0 - norm) * (value - valueBlock) * (value - valueBlock)) * norm;
          value = (value * indexOuter + valueBlock) * norm;
          // Quick return for value = 1
          if ((value >= 1.0 - DefaultQuantileEpsilon) && (variance == 0.0)) return 1.0;
          precision = a99 * sqrt(variance / (indexOuter + 1.0) / MinimumNumberOfPoints);
          if (precision < MinimumCDFEpsilon * value) return value;
          // 0.1 * ((1000 * indexOuter) / outerMax) is to print percents with one figure after the decimal point
          Log::Info(OSS() << 0.1 * ((1000 * indexOuter) / outerMax) << "% value=" << value << " absolute precision(99%)=" << precision << " relative precision(99%)=" << ((value > 0.0) ? precision / value : -1.0));
        }
00320       cdfEpsilon_ = precision;
      return value;
      } // computeCDF

      /* Compute the probability content of an interval */
      NumericalScalar Normal::computeProbability(const Interval & interval) const
      {
      if (interval.isEmpty()) return 0.0;
      const UnsignedLong dimension(getDimension());
      // Decompose and normalize the interval
      NumericalPoint lower(normalize(interval.getLowerBound()));
      NumericalPoint upper(normalize(interval.getUpperBound()));
      Interval::BoolCollection finiteLower(interval.getFiniteLowerBound());
      Interval::BoolCollection finiteUpper(interval.getFiniteUpperBound());
      /* Special treatment for independent components */
      if (hasIndependentCopula_)
        {
          NumericalScalar lowerCDF(0.0);
          if (finiteLower[0]) lowerCDF = DistFunc::pNormal(lower[0]);
          NumericalScalar upperCDF(1.0);
          if (finiteUpper[0]) upperCDF = DistFunc::pNormal(upper[0]);
          NumericalScalar value(upperCDF - lowerCDF);
          for (UnsignedLong i = 1; i < dimension; ++i)
            {
            lowerCDF = 0.0;
            if (finiteLower[0]) lowerCDF = DistFunc::pNormal(lower[0]);
            upperCDF = 1.0;
            if (finiteUpper[0]) upperCDF = DistFunc::pNormal(upper[0]);
            value *= upperCDF - lowerCDF;
            }
          return value;
        }
      /* General case */
      // For large dimension, use an adaptive integration
      if (dimension <= 500)
        {
          // Build finite/infinite flags according to MVNDST documentation
          std::vector<int> infin(dimension, -1);
          for (UnsignedLong i = 0; i < dimension; i++)
            {
            // Infin[i] should be:
            //  2 if finiteLower[i] && finiteUpper[i]
            //  1 if finiteLower[i] && !finiteUpper[i]
            //  0 if !finiteLower[i] && finiteUpper[i]
            // -1 if !finiteLower[i] && !finiteUpper[i]
            infin[i] += 2 * int(finiteLower[i]) + int(finiteUpper[i]);
            }
          NumericalPoint correl(dimension * (dimension - 1) / 2, 0.0);
          /* Copy the correlation matrix in the proper format for mvndst */
          for (UnsignedLong i = 0; i < UnsignedLong(dimension); i++)
            {
            for (UnsignedLong j = 0; j < i; j++)
              {
                correl[j + i * (i - 1) / 2] = R_(j, i);
              }
            }
          int maxpts(MinimumNumberOfPoints);
          // Use only relative precision
          double abseps(0.0);
          // Reduce the precision according to the dimension. It ranges from MaximumCDFEpsilon for dimension=4 to MinimumCDFEpsilon for dimension=500
          double releps(MaximumCDFEpsilon * pow(MinimumCDFEpsilon / MaximumCDFEpsilon, (4.0 + dimension) / 496.0));
          double error;
          double value;
          int inform;
          int dim = static_cast<UnsignedLong>( dimension );
          do
            {
            MVNDST_F77(&dim, &lower[0], &upper[0], &infin[0], &correl[0], &maxpts, &abseps, &releps, &error, &value, &inform);
            if (inform == 1)
              {
                Log::Warn(OSS() << "Warning, in Normal::computeProbability(), the required precision has not been achieved with maxpts=" << NumericalScalar(maxpts) << ", we only got an absolute error of " << error << " and a relative error of " << error / value << ". We retry with maxpoint=" << 10 * maxpts);
                maxpts *= 10;
              }
            } while ((static_cast<UnsignedLong>(maxpts) <= MaximumNumberOfPoints) && (inform == 1));
          if (inform == 1) Log::Warn(OSS() << "Warning, in Normal::computeProbability(), the required precision has not been achieved with maxpts=" << NumericalScalar(maxpts) << ", we only got an absolute error of " << error << " and a relative error of " << error / value << ". No more retry.");
          if ((value < 0.0) || (value > 1.0)) Log::Warn(OSS() << "Warning, in Normal::computeProbability(), got a value outside of [0, 1], value=" << value << " your dependence structure might be too complex. The value will be truncated.");
          return std::min(std::max(value, 0.0), 1.0);
        }
      // For very large dimension, use a MonteCarlo algorithm
      Log::Warn(OSS() << "Warning, in Normal::computeProbability(), the dimension is very high. We will use a Monte Carlo method for the computation with a relative precision of 0.1% at 99% confidence level and a maximum of " << 10.0 * MaximumNumberOfPoints << " realizations. Expect a long running time and a poor accuracy for low values of the CDF...");
      NumericalScalar value(0.0);
      NumericalScalar variance(0.0);
      NumericalScalar a99(DistFunc::qNormal(0.995));
      UnsignedLong outerMax(10 * MaximumNumberOfPoints / MinimumNumberOfPoints);
      NumericalScalar precision(0.0);
      for (UnsignedLong indexOuter = 0; indexOuter < outerMax; ++indexOuter)
        {
          NumericalScalar valueBlock(0.0);
          NumericalScalar varianceBlock(0.0);
          for (UnsignedLong indexSample = 0; indexSample < MinimumNumberOfPoints; ++indexSample)
            {
            // ind value is 1.0 if the realization is inside of the integration domain, 0.0 else.
            NumericalScalar ind(interval.isInside(getRealization()));
            NumericalScalar norm(1.0 / (indexSample + 1.0));
            varianceBlock = (varianceBlock * indexSample + (1.0 - norm) * (valueBlock - ind) * (valueBlock - ind)) * norm;
            valueBlock = (valueBlock * indexSample + ind) * norm;
            }
          NumericalScalar norm(1.0 / (indexOuter + 1.0));
          variance = (varianceBlock + indexOuter * variance + (1.0 - norm) * (value - valueBlock) * (value - valueBlock)) * norm;
          value = (value * indexOuter + valueBlock) * norm;
          // Quick return for value = 1
          if ((value >= 1.0 - DefaultQuantileEpsilon) && (variance == 0.0)) return 1.0;
          precision = a99 * sqrt(variance / (indexOuter + 1.0) / MinimumNumberOfPoints);
          if (precision < MinimumCDFEpsilon * value) return value;
          // 0.1 * ((1000 * indexOuter) / outerMax) is to print percents with one figure after the decimal point
          Log::Info(OSS() << 0.1 * ((1000 * indexOuter) / outerMax) << "% value=" << value << " absolute precision(99%)=" << precision << " relative precision(99%)=" << ((value > 0.0) ? precision / value : -1.0));
00426         }
      return value;
      }

      /* Get the characteristic function of the distribution, i.e. phi(u) = E(exp(I*u*X)) */
      NumericalComplex Normal::computeCharacteristicFunction(const NumericalPoint & point) const
      {
      if (getDimension() !=1) throw NotYetImplementedException(HERE);
00434       NumericalScalar t(point[0]);
      return exp(NumericalComplex(-0.5 * t * t * sigma_[0] * sigma_[0], -mean_[0] * t));
      }

      /* Get the CDF gradient of the distribution */
      Normal::NumericalPoint Normal::computeCDFGradient(const NumericalPoint & point) const
      {
      UnsignedLong dimension(getDimension());
      // Be carefull! integer arithmetic used! don't factorize dimension!
      NumericalPoint gradientCDF(2 * dimension);
      if (dimension == 1)
        {
          NumericalScalar pdf(computePDF(point));
          gradientCDF[0] = -pdf;
          gradientCDF[1] = -pdf * (point[0] - mean_[0]) / sigma_[0];
          return gradientCDF;
        }
00451       // To be implemented
      return gradientCDF;
      }

      /* Compute the scalar quantile of the scalar standard normal distribution */
      NumericalScalar Normal::computeScalarQuantile(const NumericalScalar p, const NumericalScalar initialGuess, const NumericalScalar initialStep) const
      {
      if (p < 0.0 || p > 1.0) throw InvalidArgumentException(HERE) << "Error: cannot compute a quantile for a probability level outside of [0, 1]";
      if (p == 0.0) return -0.5 * sqrt(4.0 * (log(SpecFunc::ISQRT2PI) - SpecFunc::LogMinNumericalScalar));
      if (p == 1.0) return 0.5 * sqrt(4.0 * (log(SpecFunc::ISQRT2PI) - SpecFunc::LogMinNumericalScalar));
      return DistFunc::qNormal(p);
      } // computeScalarQuantile

      /* Compute the PDF of Xi | X1, ..., Xi-1. x = Xi, y = (X1,...,Xi-1)
       For Normal distribution, the conditional distribution is also Normal, with mean and covariance
       such as:
00467        mean_cond = mean(x) + cov(x, y).cov(y, y)^(-1)(y - mean(y))
       cov_cond = cov(x, x) - cov(x, y).cov(y, y)^(-1)cov(x, y)
       This expression simplifies if we use the inverse of the Cholesky factor of the covariance matrix.
       See [Lebrun, Dutfoy, "Rosenblatt and Nataf transformation"]
      */
      NumericalScalar Normal::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
      NumericalScalar meanRos(0.0);
      const NumericalScalar sigmaRos(1.0 / inverseCholesky_(conditioningDimension, conditioningDimension));
      for (UnsignedLong i = 0; i < conditioningDimension; ++i)
        {
          meanRos += inverseCholesky_(conditioningDimension, i) / sqrt(sigma_[i]) * (y[i] - mean_[i]);
        }
00485       meanRos = mean_[conditioningDimension] - sigmaRos * sqrt(sigma_[conditioningDimension]) * meanRos;
      return exp(-0.5 * pow(x - meanRos, 2.0) / (sigmaRos * sigmaRos)) / (sigmaRos * sqrt(2.0 * M_PI));
      }
      
      /* Compute the CDF of Xi | X1, ..., Xi-1. x = Xi, y = (X1,...,Xi-1) */
      NumericalScalar Normal::computeConditionalCDF(const NumericalScalar x, const NumericalPoint & y) const
      {
      const UnsignedLong conditioningDimension(y.getDimension());
      if (conditioningDimension >= getDimension()) throw InvalidArgumentException(HERE) << "Error: cannot compute a conditional CDF 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)->computeCDF(x);
      // General case
      NumericalScalar meanRos(0.0);
      const NumericalScalar sigmaRos(1.0 / inverseCholesky_(conditioningDimension, conditioningDimension));
      for (UnsignedLong i = 0; i < conditioningDimension; ++i)
        {
          meanRos += inverseCholesky_(conditioningDimension, i) / sqrt(sigma_[i]) * (y[i] - mean_[i]);
        }
00503       meanRos = mean_[conditioningDimension] - sigmaRos * sqrt(sigma_[conditioningDimension]) * meanRos;
      return DistFunc::pNormal((x - meanRos) / sigmaRos);
      }
      
      /* Compute the quantile of Xi | X1, ..., Xi-1, i.e. x such that CDF(x|y) = q with x = Xi, y = (X1,...,Xi-1) */
      NumericalScalar Normal::computeConditionalQuantile(const NumericalScalar q, const NumericalPoint & y) const
      {
      const UnsignedLong conditioningDimension(y.getDimension());
      if (conditioningDimension >= getDimension()) 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]";
      // Special case when no contitioning or independent copula
      if ((conditioningDimension == 0) || hasIndependentCopula()) return mean_[conditioningDimension] + sigma_[conditioningDimension] * DistFunc::qNormal(q);
      // General case
      NumericalScalar meanRos(0.0);
      const NumericalScalar sigmaRos(1.0 / inverseCholesky_(conditioningDimension, conditioningDimension));
      for (UnsignedLong i = 0; i < conditioningDimension; ++i)
        {
          meanRos += inverseCholesky_(conditioningDimension, i) / sqrt(sigma_[i]) * (y[i] - mean_[i]);
        }
      meanRos = mean_[conditioningDimension] - sigmaRos * sqrt(sigma_[conditioningDimension]) * meanRos;
      if (q == 0.0) return meanRos - 0.5 * sigmaRos * sqrt(4.0 * (log(SpecFunc::ISQRT2PI / sigmaRos) - SpecFunc::LogMinNumericalScalar));
00524       if (q == 1.0) return meanRos + 0.5 * sigmaRos * sqrt(4.0 * (log(SpecFunc::ISQRT2PI / sigmaRos) - SpecFunc::LogMinNumericalScalar));
      return meanRos + sigmaRos * DistFunc::qNormal(q);
      }

      /* Get the i-th marginal distribution */
      Normal::Implementation Normal::getMarginal(const UnsignedLong i) const throw(InvalidArgumentException)
      {
            if (i >= getDimension()) throw InvalidArgumentException(HERE) << "The index of a marginal distribution must be in the range [0, dim-1]";
      // Special case for dimension 1
      if (getDimension() == 1) return clone();
      // General case
            CorrelationMatrix R(1);
            NumericalPoint sigma(1, sigma_[i]);
            NumericalPoint mean(1, mean_[i]);
      Normal result(mean, sigma, R);
      result.setName("Marginal of " + getName());
00540       result.setWeight(getWeight());
            return new Normal(result);
      }

      /* Get the distribution of the marginal distribution corresponding to indices dimensions */
      Normal::Implementation Normal::getMarginal(const Indices & indices) const throw(InvalidArgumentException)
      {
            UnsignedLong dimension(getDimension());
      if (!indices.check(dimension - 1)) throw InvalidArgumentException(HERE) << "The indices of a marginal distribution must be in the range [0, dim-1] and  must be different";
      // Special case for dimension 1
      if (dimension == 1) return clone();
      // General case
            const UnsignedLong outputDimension(indices.getSize());
            CorrelationMatrix R(outputDimension);
            NumericalPoint sigma(outputDimension);
            NumericalPoint mean(outputDimension);
      // Extract the correlation matrix, the marginal standard deviations and means
      for (UnsignedLong i = 0; i < outputDimension; ++i)
        {
          const UnsignedLong index_i(indices[i]);
          sigma[i] = sigma_[index_i];
          mean[i] = mean_[index_i];
          for (UnsignedLong j = 0; j <= i; j++)
            {
            R(i, j) = R_(index_i, indices[j]);
            }
        }
      Normal result(mean, sigma, R);
      result.setName("Marginal of " + getName());
00569       result.setWeight(getWeight());
            return new Normal(result);
      } // getMarginal(Indices)

      /* Get the skewness of the distribution */
      Normal::NumericalPoint Normal::getSkewness() const throw(NotDefinedException)
00575       {
      return NumericalPoint(getDimension(), 0.0);
      }

      /* Get the roughness, i.e. the L2-norm of the PDF */
      NumericalScalar Normal::getRoughness() const
      {
00582       // 0.2820947917738781434740398 = 1 / (2 * sqrt(Pi))
      return 0.2820947917738781434740398 / getSigma()[0];
      }

      /* Get the kurtosis of the distribution */
      Normal::NumericalPoint Normal::getKurtosis() const throw(NotDefinedException)
      {
00589       return NumericalPoint(getDimension(), 3.0);
      }

      /* Get the standard distribution, i.e. a distribution of the same kind but with zero mean,
       * unit marginal standard distribution and identity correlation */
00594       Normal::Implementation Normal::getStandardDistribution() const
      {
      return new Normal(getDimension());
      }

      Normal::Implementation Normal::getCopula() const
00600       {
      return new NormalCopula(R_);
      }

      /* Compute the radial distribution CDF */
      NumericalScalar Normal::computeRadialDistributionCDF(const NumericalScalar radius) const
      {
      // It should be a chi-squared distribution with dimension degrees of freedom
      // We only have a gamma distribution available for this time (OpenTURNS rc1)
00609       // so we use the relation chisq(n) = gamma(n/2, 1/2)
      return DistFunc::pGamma(0.5 * getDimension(), 0.5 * radius * radius);
      }

      /* Sigma accessor */
      void Normal::setSigma(const NumericalPoint & sigma)
      throw(InvalidArgumentException)
      {
00617       EllipticalDistribution::setSigma(sigma);
      sigma_ = sigma;
      }

      /* Correlation matrix accessor */
      void Normal::setCorrelationMatrix(const CorrelationMatrix & R)
      throw(InvalidArgumentException)
      {
      EllipticalDistribution::setCorrelationMatrix(R);
00626       R_ = R;
      checkIndependentCopula();
      }

      /* Tell if the distribution has independent copula */
      Bool Normal::hasIndependentCopula() const
00632       {
      return hasIndependentCopula_;
      }

      /* Check if the distribution has independent copula */
      void Normal::checkIndependentCopula()
      {
      hasIndependentCopula_ = true;
      UnsignedLong dimension(getDimension());
      for (UnsignedLong i = 0; i < dimension; i++)
        {
          for (UnsignedLong j = 0; j < i; j++)
            {
            if (R_(i, j) != 0.0)
              {
                hasIndependentCopula_ = false;
                return;
              }
00650             }
        }
      } // checkIndependentCopula

      /* Method save() stores the object through the StorageManager */
      void Normal::save(const StorageManager::Advocate & adv) const
      {
      EllipticalDistribution::save(adv);
00658       adv.writeValue("normalizationFactor_", normalizationFactor_);
      adv.writeValue("hasIndependentCopula_", hasIndependentCopula_);
      }

      /* Method load() reloads the object from the StorageManager */
      void Normal::load(const StorageManager::Advocate & adv)
      {
      EllipticalDistribution::load(adv);

      String name;
      NumericalScalar numericalScalarValue;
      StorageManager::List objList = adv.getList(StorageManager::NumericalScalarEntity);
      for(objList.firstValueToRead(); objList.moreValuesToRead(); objList.nextValueToRead()) {
        if (objList.readValue(name, numericalScalarValue)) {
          if (name == "normalizationFactor_") normalizationFactor_ = numericalScalarValue;
        }
      }
      Bool boolValue;
      objList = adv.getList(StorageManager::BoolEntity);
      for(objList.firstValueToRead(); objList.moreValuesToRead(); objList.nextValueToRead()) {
        if (objList.readValue(name, boolValue)) {
          if (name == "hasIndependentCopula_") hasIndependentCopula_ = boolValue;
        }
      }
      computeRange();
      }
      
    } /* namespace Distribution */
  } /* namespace Uncertainty */
} /* namespace OpenTURNS */

Generated by  Doxygen 1.6.0   Back to index