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

EllipticalDistribution.cxx

Go to the documentation of this file.
//                                               -*- C++ -*-
/**
 *  @file  EllipticalDistribution.cxx
 *  @brief Abstract top-level class for elliptical 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: EllipticalDistribution.cxx 995 2008-10-31 10:52:04Z dutka $
 */
#include <cmath>
#include "EllipticalDistribution.hxx"
#include "NatafEllipticalDistributionEvaluation.hxx"
#include "NatafEllipticalDistributionGradient.hxx"
#include "NatafEllipticalDistributionHessian.hxx"
#include "InverseNatafEllipticalDistributionEvaluation.hxx"
#include "InverseNatafEllipticalDistributionGradient.hxx"
#include "InverseNatafEllipticalDistributionHessian.hxx"
#include "Description.hxx"
#include "PersistentObjectFactory.hxx"
#include "Lapack.hxx"
#include "IdentityMatrix.hxx"

namespace OpenTURNS {

  namespace Uncertainty {

    namespace Model {

      typedef Base::Type::IdentityMatrix                              IdentityMatrix;
00046       typedef Base::Type::Description                                 Description;
      typedef Base::Common::NotYetImplementedException                NotYetImplementedException;
      typedef Base::Common::NotDefinedException                       NotDefinedException;
      typedef Algorithm::NatafEllipticalDistributionEvaluation        NatafEllipticalDistributionEvaluation;
      typedef Algorithm::NatafEllipticalDistributionGradient          NatafEllipticalDistributionGradient;
      typedef Algorithm::NatafEllipticalDistributionHessian           NatafEllipticalDistributionHessian;
      typedef Algorithm::InverseNatafEllipticalDistributionEvaluation InverseNatafEllipticalDistributionEvaluation;
      typedef Algorithm::InverseNatafEllipticalDistributionGradient   InverseNatafEllipticalDistributionGradient;
00054       typedef Algorithm::InverseNatafEllipticalDistributionHessian    InverseNatafEllipticalDistributionHessian;

      CLASSNAMEINIT(EllipticalDistribution);

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

      /* Default constructor */
      EllipticalDistribution::EllipticalDistribution(const String & name)
      : ContinuousDistribution(name)
      {
      // Initialize any other class members here
      // At last, allocate memory space if needed, but go to destructor to free it
      }

      /* Parameter constructor */
      EllipticalDistribution::EllipticalDistribution(const NumericalPoint & mean,
                                         const NumericalPoint & sigma,
                                         const CorrelationMatrix & R,
                                         const String & name)
      : ContinuousDistribution(name),
        R_(R),
        sigma_(sigma),
        mean_(mean)
      {
      const UnsignedLong dimension(R.getDimension());
      // We check if the inputs have the same dimension. If not, we throw an exception
      if ( (dimension    != mean.getDimension()) ||
           (dimension != sigma.getDimension()) )
        throw InvalidArgumentException(HERE)
          << "Arguments have incompatible dimensions: R dimension=" << dimension
          << " sigma dimension=" << sigma.getDimension()
          << " mean dimension=" << mean.getDimension();     
      // We check that the given correlation matrix is definite positive
      if ( !R.isPositiveDefinite()) throw InvalidArgumentException(HERE) << "The correlation matrix must be definite positive R=" << R;
      // We check that the marginal standard deviations are > 0
      for(UnsignedLong i = 0; i < dimension; ++i)
        {
          if (sigma[i] <= 0.0) throw InvalidArgumentException(HERE) << "The marginal standard deviations must be > 0 sigma=" << sigma;
        }
00093       // Then we set the dimension of the Elliptical distribution
      setDimension(dimension);
      // We initialize the description
      Description description(dimension);
      // Set default description
      for (UnsignedLong i = 0; i < dimension; ++i)
00099         {
          description[i] = OSS() << "marginal " << i + 1;
        }
      setDescription(description);
      // Compute the auxiliary attributes
      update();
      }

      /* Destructor */
      EllipticalDistribution::~EllipticalDistribution() {
      // Remember to free any allocated memory space
      // Remember NOT to throw any exception here
      }

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

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

      } else sameObject = true;

00125       return sameObject;
      }
  
      /* Centers and reduces a value u = Diag(sigma_)^(-1) * (x - mean_) */
      EllipticalDistribution::NumericalPoint EllipticalDistribution::normalize(const NumericalPoint & x) const
      {
      NumericalPoint u(x - mean_);
      for (UnsignedLong i = 0; i < getDimension(); ++i)
        {
          u[i] /= sigma_[i];
        }
00136       return u;
      }

      /* Decenters and scales a value x = mean_ + Diag(sigma_) * u */
      EllipticalDistribution::NumericalPoint EllipticalDistribution::denormalize(const NumericalPoint & u) const
      {
      NumericalPoint x(mean_);
      for (UnsignedLong i = 0; i < getDimension(); ++i)
00144         {
          x[i] += sigma_[i] * u[i];
        }
      return x;
      }

00150       /* String converter */
      String EllipticalDistribution::str() const {
      OSS oss;
      oss << "class=" << EllipticalDistribution::GetClassName()
          << "parameters collection=" << getParametersCollection();
      return oss;
00156       }
  
      /* Compute the numerical range of the distribution given the parameters values */
      void EllipticalDistribution::computeRange()
      {
      throw NotYetImplementedException(HERE);
      }

      /* Tell if the distribution is elliptical */
00165       Bool EllipticalDistribution::isElliptical() const
      {
      return true;
      }
    
      /* Tell if the distribution has elliptical copula */
00171       Bool EllipticalDistribution::hasEllipticalCopula() const
      {
      return true;
      }


      /* Compute the density generator of the elliptical distribution, i.e.
00178        *  the function phi such that the density of the distribution can
       *  be written as p(x) = phi((x-mean_).C^{-1} * (x-mean_))                      */
      NumericalScalar EllipticalDistribution::computeDensityGenerator(const NumericalScalar betaSquare) const
      {
      throw NotYetImplementedException(HERE);
      }

      /* Compute the derivative of the density generator */
      NumericalScalar EllipticalDistribution::computeDensityGeneratorDerivative(const NumericalScalar betaSquare) const
      {
      // Use centered finite difference
      return (computeDensityGenerator(betaSquare + DefaultPDFEpsilon) - computeDensityGenerator(betaSquare - DefaultPDFEpsilon)) / DefaultPDFEpsilon;
      }

      /* Compute the seconde derivative of the density generator */
      NumericalScalar EllipticalDistribution::computeDensityGeneratorSecondDerivative(const NumericalScalar betaSquare) const
      {
      // Use centered finite difference
      return (computeDensityGenerator(betaSquare + DefaultPDFEpsilon) - 2.0 * computeDensityGenerator(betaSquare) + computeDensityGenerator(betaSquare - DefaultPDFEpsilon)) / (DefaultPDFEpsilon * DefaultPDFEpsilon);
      }

      /* Get the DDF of the distribution */
      EllipticalDistribution::NumericalPoint EllipticalDistribution::computeDDF(const NumericalPoint & point) const
00201       {
      const NumericalPoint iLx(inverseCholesky_ * (point - mean_));
      const NumericalScalar betaSquare(iLx.norm2());
      return 2.0 * normalizationFactor_ * computeDensityGeneratorDerivative(betaSquare) * inverseCholesky_.transpose() * iLx;
      }

      /* Get the PDF of the distribution */
      NumericalScalar EllipticalDistribution::computePDF(const NumericalPoint & point) const
      {
      const NumericalPoint iLx(inverseCholesky_ * (point - mean_));
      const NumericalScalar betaSquare(iLx.norm2());
      return normalizationFactor_ * computeDensityGenerator(betaSquare);
      }

      /* Get the PDF gradient of the distribution */
      EllipticalDistribution::NumericalPoint EllipticalDistribution::computePDFGradient(const NumericalPoint & point) const
      {
      const NumericalPoint minusGardientMean(computeDDF(point));
      const UnsignedLong dimension(getDimension());
      const NumericalPoint u(normalize(point));
      const NumericalPoint iRu(inverseR_ * u);
      const NumericalScalar betaSquare(NumericalPoint::dot(u, iRu));
00223       const NumericalScalar phi(computeDensityGenerator(betaSquare));
      const NumericalScalar phiDerivative(computeDensityGeneratorDerivative(betaSquare));
      NumericalPoint pdfGradient(2 * dimension);
      for (UnsignedLong i = 0; i < dimension; ++i)
        {
          NumericalScalar iSigma(1.0 / sigma_[i]);
          // dPDF / dmu_i
          pdfGradient[i] = -2.0 * normalizationFactor_ * phiDerivative * iRu[i] * iSigma;
          // dPDF / dsigma_i
          pdfGradient[dimension + i] = pdfGradient[i] * u[i] - normalizationFactor_ * phi * iSigma;
        }
      return pdfGradient;
      }

      /* Update the derivative attributes */
      void EllipticalDistribution::update()
      {
      const UnsignedLong dim(getDimension());
      covariance_ = CovarianceMatrix(dim);
      // Build the covariance matrix
      covariance_ = R_;
      for (UnsignedLong i = 0; i < dim; ++i)
        {
          for (UnsignedLong j = 0; j <= i; ++j)
            {
            covariance_(i, j) *= sigma_[i] * sigma_[j];
            }
        }
      // Compute its Cholesky factor
      cholesky_ = covariance_.computeCholesky();
      // Inversion of the Cholesky factor using the dtrsm blas level 3 routine
      // side tells if we solve M.X = alpha.B or X.M = alpha.B
      char side('L');
      int lside(1);
      // M must be triangular. uplo tells if it is upper or lower triangular
      char uplo('L');
      int luplo(1);
      // transa tells if M is transposed or not
      char transa('N');
      int ltransa(1);
      // diag tells if M is unit diagonal or not
      char diag('N');
      int ldiag(1);
      // the row dimension of M
      int m(dim);
      // the column dimension of M
      int n(dim);
      // we solve the case alpha=1
      double alpha(1.0);
      // leading dimension of M
      int lda(dim);
      // leading dimension of B
      int ldb(dim);
      // As we want to inverse M, we set B = Id
      inverseCholesky_ = IdentityMatrix(dim);
      // B stores the result of the routine
      DTRSM_F77(&side, &uplo, &transa, &diag, &m, &n, &alpha, const_cast<double*>(&((*cholesky_.getImplementation())[0])), &lda, const_cast<double*>(&((*inverseCholesky_.getImplementation())[0])), &ldb, &lside, &luplo, &ltransa, &ldiag);

      // Inverse the correlation matrix R = D^(-1).L.L'.D^(-1)
      // R^(-1) = D.L^(-1).L^(-1)'.D
      inverseR_ = SymmetricMatrix(dim);
      const SquareMatrix inverseCovariance(inverseCholesky_.transpose() * inverseCholesky_);
00285       for (UnsignedLong i = 0; i < dim; ++i)
        {
          for (UnsignedLong j = 0; j <= i; ++j)
            {
            inverseR_(i, j) = sigma_[i] * inverseCovariance(i, j) * sigma_[j];
            }
        }
      normalizationFactor_ = 1.0;
      for (UnsignedLong i = 0; i < dim; ++i)
        {
          normalizationFactor_ /= cholesky_(i, i);
        }
      }

      /* Get the quantile of the distribution */
      EllipticalDistribution::NumericalPoint EllipticalDistribution::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 boarding values
      if (prob == 0.0) return getRange().getLowerBound();
      if (prob == 1.0) return getRange().getUpperBound();
      const UnsignedLong dim(getDimension());
      // Special case for dimension 1
      if (dim == 1) return denormalize(NumericalPoint(1, computeScalarQuantile(prob)));

      NumericalScalar left(computeScalarQuantile(prob));
      NumericalScalar right(computeScalarQuantile(1.0 - (1.0 - prob) / dim));

      NumericalScalar leftCDF(computeCDF(denormalize(NumericalPoint(dim, left))));
      NumericalScalar rightCDF(computeCDF(denormalize(NumericalPoint(dim, right))));
      while ((rightCDF - leftCDF > 2.0 * cdfEpsilon_) && (right - left > 2.0 * cdfEpsilon_))
        {
          const NumericalScalar middle(0.5 * (left + right));
          const NumericalScalar middleCDF(computeCDF(denormalize(NumericalPoint(dim, middle))));
00319           if (middleCDF > prob)
            {
            right = middle;
            rightCDF = middleCDF;
            }
          else
            {
            left = middle;
            leftCDF = middleCDF;
            }
        }
      return denormalize(NumericalPoint(dim, 0.5 * (left + right)));
      } // computeQuantile
00332 
      /* Mean point accessor */
      void EllipticalDistribution::setMean(const NumericalPoint & mean)
      throw(InvalidArgumentException)
      {
      if (mean.getDimension() != getDimension())
00338         throw InvalidArgumentException(HERE)
          << "Mean point dimension (" << mean.getDimension()
          << ") differ from distribution dimension(" << getDimension()
          << "). Unable to construct EllipticalDistribution distribution object.";
      mean_ = mean;
      computeRange();
00344       }

      /* Mean point accessor */
      EllipticalDistribution::NumericalPoint EllipticalDistribution::getMean() const throw(NotDefinedException)
      {
      return mean_;
      }

      /* Covariance matrix accessor */
      EllipticalDistribution::CovarianceMatrix EllipticalDistribution::getCovariance() const throw(NotDefinedException)
      {
      return covariance_;
      }

      /* Sigma accessor */
      void EllipticalDistribution::setSigma(const NumericalPoint & sigma)
      throw(InvalidArgumentException)
      {
      if (sigma.getDimension() != getDimension())
        throw InvalidArgumentException(HERE)
00364           << "Sigma dimension (" << sigma.getDimension()
          << ") differ from distribution dimension(" << getDimension()
          << "). Unable to construct elliptical distribution object.";
      
      // We check that the marginal standard deviations are > 0
      for(UnsignedLong i = 0; i < sigma.getDimension(); ++i)
        {
          if (sigma[i] <= 0.0) throw InvalidArgumentException(HERE) << "The marginal standard deviations must be > 0 sigma=" << sigma;
        }
00373       sigma_ = sigma;
      update();
      computeRange();
      }

      /* Sigma accessor */
00379       EllipticalDistribution::NumericalPoint EllipticalDistribution::getSigma() const
      {
      return sigma_;
      }

      /* Get the standard deviation of the distribution.
         Warning! This method MUST be overloaded for elliptical distributions without finite second moment:
       it is possible to have a well-defined sigma vector but no standard deviation, think about Stundent
       distribution with nu < 2 */
      EllipticalDistribution::NumericalPoint EllipticalDistribution::getStandardDeviation() const throw(NotDefinedException)
      {
      return sigma_;
      }

      /* Correlation matrix accessor */
      void EllipticalDistribution::setCorrelationMatrix(const CorrelationMatrix & R)
      throw(InvalidArgumentException)
00396       {
      if (R.getDimension() != getDimension())
        throw InvalidArgumentException(HERE)
          << "Correlation Matrix dimension (" << R.getDimension()
          << ") differ from distribution dimension(" << getDimension()
          << "). Unable to construct elliptical distribution object.";
00402       
      // We check that the given correlation matrix is definite positive
      if ( !R.isPositiveDefinite()) throw InvalidArgumentException(HERE) << "The correlation matrix must be definite positive R=" << R;

      R_ = R;
      update();
00408       }

      /* Correlation matrix accessor */
      EllipticalDistribution::CorrelationMatrix EllipticalDistribution::getCorrelationMatrix() const
      {
      return R_;
00414       }

      /* Inverse correlation matrix accessor */
      EllipticalDistribution::SquareMatrix EllipticalDistribution::getInverseCorrelation() const
      {
      return inverseR_;
00420       }

      /* Cholesky factor of the correlation matrix accessor */
      EllipticalDistribution::SquareMatrix EllipticalDistribution::getCholesky() const
      {
      return cholesky_;
00426       }

      /* Inverse of the Cholesky factor of the correlation matrix accessor */
      EllipticalDistribution::SquareMatrix EllipticalDistribution::getInverseCholesky() const
      {
      return inverseCholesky_;
      }

      /* Virtual copy constructor */
      EllipticalDistribution * EllipticalDistribution::clone() const
      {
      return new EllipticalDistribution(*this);
      }

      /* Get the isoprobabilist transformation */
      EllipticalDistribution::IsoProbabilisticTransformation EllipticalDistribution::getIsoProbabilisticTransformation() const
      {
      IsoProbabilisticTransformation transform;
      transform.setEvaluationImplementation(new NatafEllipticalDistributionEvaluation(mean_, inverseCholesky_));
      transform.setGradientImplementation(new NatafEllipticalDistributionGradient(inverseCholesky_));
      transform.setHessianImplementation(new NatafEllipticalDistributionHessian(getDimension()));
      // Set the parameters values and descriptions
      // The result of parametersGradient is given
      // in the following form:
      // (d/dmu, d/dsigma)
      // There is no gradient according to the dependence parameters yet (28/10/2006)
      const UnsignedLong dimension(getDimension());
      NumericalPointWithDescription parameters(2 * dimension);
      Description description(parameters.getDimension());
      NumericalPointWithDescriptionCollection parametersCollection(getParametersCollection());
00456       for (UnsignedLong i = 0; i < dimension; ++i)
        {
          const Description parametersDescription(parametersCollection[i].getDescription());
          const String marginalName(parametersCollection[i].getName());
          parameters[i] = parametersCollection[i][0];
          parameters[dimension + i] = parametersCollection[i][1];
          description[i] = OSS() << marginalName << "_" << parametersDescription[0];
          description[dimension + i] = OSS() << marginalName << "_" << parametersDescription[1];
        }
      parameters.setDescription(description);
      transform.setParameters(parameters);
      return transform;
      }

      /* Get the inverse isoprobabilist transformation */
      EllipticalDistribution::InverseIsoProbabilisticTransformation EllipticalDistribution::getInverseIsoProbabilisticTransformation() const
      {
      InverseIsoProbabilisticTransformation inverseTransform;
      inverseTransform.setEvaluationImplementation(new InverseNatafEllipticalDistributionEvaluation(mean_, cholesky_));
      inverseTransform.setGradientImplementation(new InverseNatafEllipticalDistributionGradient(cholesky_));
      inverseTransform.setHessianImplementation(new InverseNatafEllipticalDistributionHessian(getDimension()));
      // Set the parameters values and descriptions
      // The result of parametersGradient is given
      // in the following form:
      // (d/dmu, d/dsigma)
      // There is no gradient according to the dependence parameters yet (28/10/2006)
      const UnsignedLong dimension(getDimension());
      NumericalPointWithDescription parameters(2 * dimension);
      Description description(parameters.getDimension());
      NumericalPointWithDescriptionCollection parametersCollection(getParametersCollection());
00486       for (UnsignedLong i = 0; i < dimension; ++i)
        {
          const Description parametersDescription(parametersCollection[i].getDescription());
          const String marginalName(parametersCollection[i].getName());
          parameters[i] = parametersCollection[i][0];
          parameters[dimension + i] = parametersCollection[i][1];
          description[i] = OSS() << marginalName << "_" << parametersDescription[0];
          description[dimension + i] = OSS() << marginalName << "_" << parametersDescription[1];
        }
      parameters.setDescription(description);
      inverseTransform.setParameters(parameters);
      return inverseTransform;
      }

      /* Parameters value and description accessor */
      EllipticalDistribution::NumericalPointWithDescriptionCollection EllipticalDistribution::getParametersCollection() const
      {
      const UnsignedLong dimension(getDimension());
      NumericalPointWithDescriptionCollection parameters(dimension + (dimension > 1 ? 1 : 0));
      // First put the marginal parameters
      Description description(getDescription());
      for (UnsignedLong marginalIndex = 0; marginalIndex < dimension; ++marginalIndex)
        {
          NumericalPointWithDescription point(2);
          Description marginalDescription(point.getDimension());
          point[0] = mean_[marginalIndex];
          point[1] = sigma_[marginalIndex];
          marginalDescription[0] = "mean";
          marginalDescription[1] = "standard deviation";
          point.setDescription(marginalDescription);
          point.setName(description[marginalIndex]);
          parameters[marginalIndex] = point;
        } // marginalIndex
#ifdef WITH_DEPENDENCE_SENSITIVITY
      if (dimension > 1)
        {
          // Second put the dependence parameters
          NumericalPoint point(dimension * (dimension - 1) / 2);
          Description marginalDescription(point.getDimension());
          OSS oss;
          oss << "dependence";
          point.setName(oss);
          UnsignedLong dependenceIndex(0);
          for (UnsignedLong i = 0; i < dimension; ++i)
            {
00531             for (UnsignedLong j = 0; j < i; ++j)
              {
                point[dependenceIndex] = R_(i, j);
                marginalDescription[dependenceIndex] = OSS() << "R_" << i + 1 << "_" << j + 1;
                ++dependenceIndex;
              }
            }
          point.setDescription(marginalDescription);
          parameters[dimension] = point;
        } // dimension > 1
#endif
      return parameters;
      } // getParametersCollection

00545       /* Method save() stores the object through the StorageManager */
      void EllipticalDistribution::save(const StorageManager::Advocate & adv) const
      {
      DistributionImplementation::save(adv);
      adv.writeValue(R_, StorageManager::MemberNameAttribute, "R_");
      adv.writeValue(sigma_, StorageManager::MemberNameAttribute, "sigma_");
      adv.writeValue(mean_, StorageManager::MemberNameAttribute, "mean_");
      adv.writeValue(covariance_, StorageManager::MemberNameAttribute, "covariance_");
      adv.writeValue(inverseR_, StorageManager::MemberNameAttribute, "inverseR_");
      adv.writeValue(cholesky_, StorageManager::MemberNameAttribute, "cholesky_");
      adv.writeValue(inverseCholesky_, StorageManager::MemberNameAttribute, "inverseCholesky_");
      adv.writeValue("normalizationFactor_", normalizationFactor_);
      }

      /* Method load() reloads the object from the StorageManager */
      void EllipticalDistribution::load(const StorageManager::Advocate & adv)
      {
      DistributionImplementation::load(adv);
      adv.readValue(R_, StorageManager::MemberNameAttribute, "R_");
      adv.readValue(sigma_, StorageManager::MemberNameAttribute, "sigma_");
      adv.readValue(mean_, StorageManager::MemberNameAttribute, "mean_");
      adv.readValue(covariance_, StorageManager::MemberNameAttribute, "covariance_");
      adv.readValue(inverseR_, StorageManager::MemberNameAttribute, "inverseR_");
      adv.readValue(cholesky_, StorageManager::MemberNameAttribute, "cholesky_");
      adv.readValue(inverseCholesky_, StorageManager::MemberNameAttribute, "inverseCholesky_");
      String name;
      NumericalScalar scalarValue;
      StorageManager::List objList = adv.getList(StorageManager::NumericalScalarEntity);
      for(objList.firstValueToRead(); objList.moreValuesToRead(); objList.nextValueToRead()) {
        if (objList.readValue(name, scalarValue)) {
          if (name == "normalizationFactor_") normalizationFactor_ = scalarValue;
        }
      }
      }

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

Generated by  Doxygen 1.6.0   Back to index