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

ComposedDistribution.cxx

Go to the documentation of this file.
//                                               -*- C++ -*-
/**
 *  @file  ComposedDistribution.cxx
 *  @brief Abstract top-level class for all ComposedDistributions
 *
 *  (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: ComposedDistribution.cxx 995 2008-10-31 10:52:04Z dutka $
 */
#include <cmath>

#include "ComposedDistribution.hxx"
#include "MarginalTransformationEvaluation.hxx"
#include "MarginalTransformationGradient.hxx"
#include "MarginalTransformationHessian.hxx"
#include "InverseMarginalTransformationEvaluation.hxx"
#include "InverseMarginalTransformationGradient.hxx"
#include "InverseMarginalTransformationHessian.hxx"
#include "NumericalMathFunctionImplementation.hxx"
#include "Indices.hxx"
#include "PersistentObjectFactory.hxx"
#include "Uniform.hxx"
#include "IndependentCopula.hxx"
#include "Log.hxx"

namespace OpenTURNS {

  namespace Uncertainty {

    namespace Distribution {

      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 Base::Func::NumericalMathFunctionImplementation               NumericalMathFunctionImplementation;
      typedef NumericalMathFunctionImplementation::EvaluationImplementation EvaluationImplementation;
      typedef NumericalMathFunctionImplementation::GradientImplementation   GradientImplementation;
      typedef NumericalMathFunctionImplementation::HessianImplementation    HessianImplementation;
      typedef Base::Type::Indices                                           Indices;
00059       typedef Base::Common::Log                                             Log;

      CLASSNAMEINIT(ComposedDistribution);

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

      /* Default constructor */
      ComposedDistribution::ComposedDistribution()
      : DistributionImplementation("ComposedDistribution"),
        distributionCollection_(),
        copula_(IndependentCopula(1))
      {
00071       setDimension(1);
      DistributionCollection coll(1);
      coll[0] = Uniform();
      setDistributionCollection( coll );
      }

      /* Default constructor */
      ComposedDistribution::ComposedDistribution(const DistributionCollection & coll)
      : DistributionImplementation("ComposedDistribution"),
        distributionCollection_(),
        copula_(IndependentCopula(coll.getSize()))
      {
      setDimension(coll.getSize());
      // We can NOT set distributionCollection_ in the member area of the constructor
00085       // because we must check before if the collection is valid (ie, if all the
      // distributions of the collection have the dimension 1). We do this by calling
      // the setDistributionCollection() method that do it for us.
      setDistributionCollection( coll );
      }

      /* Default constructor */
      ComposedDistribution::ComposedDistribution(const DistributionCollection & coll,
                                     const Copula & copula)
      throw (InvalidArgumentException)
      : DistributionImplementation("ComposedDistribution"),
        distributionCollection_(),
        copula_(copula)
      {
      setDimension(copula.getDimension());
      // We can NOT set distributionCollection_ in the member area of the constructor
00101       // because we must check before if the collection is valid (ie, if all the
      // distributions of the collection have the dimension 1). We do this by calling
      // the setDistributionCollection() method that do it for us.
      setDistributionCollection( coll );
      }

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

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

      return sameObject;
      }
  
      /* String converter */
      String ComposedDistribution::str() const
      {
      OSS oss;
      oss << "class=" << ComposedDistribution::GetClassName()
          << " name=" << getName()
          << " dimension=" << getDimension()
          << " copula=" << copula_;
      for (UnsignedLong i = 0; i < getDimension(); ++i)
00128         oss << " marginal[" << i << "]=" << distributionCollection_[i];
      return oss;
      }
  


      /* Distribution collection accessor */
      void ComposedDistribution::setDistributionCollection(const DistributionCollection & coll)
      throw (InvalidArgumentException)
      {
      // Check if the collection is not empty
      const UnsignedLong size(coll.getSize());
      if ((getDimension() != 0) && (size != getDimension())) throw InvalidArgumentException(HERE) << "The distribution collection must have a size equals to the distribution dimension";
      Description description(size);
      NumericalPoint lowerBound(size);
      NumericalPoint upperBound(size);
      Interval::BoolCollection finiteLowerBound(size);
      Interval::BoolCollection finiteUpperBound(size);
      if (size == 0) throw InvalidArgumentException(HERE) << "Collection of distributions is empty";
      // First, check that all the marginal distributions are of dimension 1
      for (UnsignedLong i = 0; i < size; ++i)
        {
          if (coll[i].getDimension() != 1) throw InvalidArgumentException(HERE) << "The marginal distribution " << i << " is of dimension " << coll[i].getDimension() << ", which is different from 1.";
          const Interval marginalRange(coll[i].getRange());
          lowerBound[i] = marginalRange.getLowerBound()[0];
          upperBound[i] = marginalRange.getUpperBound()[0];
          finiteLowerBound[i] = marginalRange.getFiniteLowerBound()[0];
          finiteUpperBound[i] = marginalRange.getFiniteUpperBound()[0];
          // The description of the ComposedDistribution is built first by using the marginal names
          // then by using the marginal description if the name is equal to OT::DefaultName.
          // Any distribution has a nonempty description by construction.
          const String marginalName(coll[i].getName());
          if (marginalName != OT::DefaultName)
            {
            description[i] = marginalName;
            }
          else
            {
            Log::Info(OSS() << "Warning: using description of the marginal " << i << " instead of its name for building the description of the ComposedDistribution, because the name is " << marginalName);
            description[i] = coll[i].getDescription()[0];
            }
        }
      // Everything is ok, store the collection
        distributionCollection_ = coll;
      isAlreadyComputedMean_ = false;
      isAlreadyComputedCovariance_ = false;
      setDescription(description);
      setRange(Interval(lowerBound, upperBound, finiteLowerBound, finiteUpperBound));
      }


      /* Distribution collection accessor */
00180       const ComposedDistribution::DistributionCollection & ComposedDistribution::getDistributionCollection() const
      {
      return distributionCollection_;
      }


      /* Copula accessor */
      void ComposedDistribution::setCopula(const Copula & copula)
      throw (InvalidArgumentException)
      {
      // We check if the copula has a dimension compatible with the one of the object,
      // especially if the object has already been created and has a collection of distribution
      if (getDimension() != 0) {
        if (getDimension() != copula.getDimension())
          throw InvalidArgumentException(HERE) << "Copula has a dimension different from the ComposedDistribution's";
      } else
        setDimension(copula.getDimension());
      
      copula_ = copula;
00199       isAlreadyComputedCovariance_ = false;
      // We ensure that the copula has the same description than the ComposedDistribution
      copula_.setDescription(getDescription());
      }


00205       /* Copula accessor */
      ComposedDistribution::Implementation ComposedDistribution::getCopula() const
      {
      return copula_.getImplementation();
      }

00211       /* Virtual constructor */
      ComposedDistribution * ComposedDistribution::clone() const
      {
      return new ComposedDistribution(*this);
      }

      /* Get one realization of the ComposedDistribution */
      ComposedDistribution::NumericalPoint ComposedDistribution::getRealization() const
      {
      const UnsignedLong dimension(getDimension());
      // Special case for independent copula
      NumericalPoint result(dimension);
      if (hasIndependentCopula())
        {
          for (UnsignedLong i = 0; i < dimension; ++i)
            {
            result[i] = distributionCollection_[i].getRealization()[0];
            }
          return result;
        }
      // General case
      NumericalPoint realization(copula_.getRealization());
      for (UnsignedLong i = 0; i < dimension; ++i)
        {
          result[i] = distributionCollection_[i].computeQuantile(realization[i])[0];
        }
      return result;
      }
     
      /* Get the DDF of the ComposedDistribution */
      ComposedDistribution::NumericalPoint ComposedDistribution::computeDDF(const NumericalPoint & point) const 
      {
      /* PDF = PDF_copula(CDF_dist1(p1), ..., CDF_distn(pn))xPDF_dist1(p1)x...xPDF_distn(pn) */
      const UnsignedLong dimension(getDimension());
      if(point.getDimension() != dimension) throw InvalidArgumentException(HERE) << "The given point has a wrong dimension";
      NumericalPoint uPoint(dimension);
      NumericalPoint pdfMarginal(dimension);
      NumericalPoint ddfMarginal(dimension);
      NumericalScalar productPDF(1.0);
      for (UnsignedLong i = 0; i < dimension; ++i)
        {
          const NumericalPoint component(NumericalPoint(1, point[i]));
          uPoint[i] = distributionCollection_[i].computeCDF(component);
          pdfMarginal[i] = distributionCollection_[i].computePDF(component);
          ddfMarginal[i] = distributionCollection_[i].computeDDF(component)[0];
          productPDF *= pdfMarginal[i];
        }
      // Initialization with the values of an independent copula
      NumericalScalar pdfCopula(1.0);
      NumericalPoint ddfCopula(dimension, 0.0);
      // If the distribution does not have an independent copula
      if (!hasIndependentCopula())
        {
          pdfCopula = copula_.computePDF(uPoint);
          ddfCopula = copula_.computeDDF(uPoint);
        }
      // Compute the ddf
      NumericalPoint ddf(dimension);
      for (UnsignedLong i = 0; i < dimension; ++i)
        {
          ddf[i] = productPDF * (ddfCopula[i] * pdfMarginal[i] + pdfCopula * ddfMarginal[i] / pdfMarginal[i]);
        }
      return ddf;
      }

      /* Get the PDF of the ComposedDistribution */
      NumericalScalar ComposedDistribution::computePDF(const NumericalPoint & point) const 
      {
      /* PDF = PDF_copula(CDF_dist1(p1), ..., CDF_distn(pn))xPDF_dist1(p1)x...xPDF_distn(pn) */
      const UnsignedLong dimension(getDimension());
      if(point.getDimension() != dimension) throw InvalidArgumentException(HERE) << "The given point has a wrong dimension";
      NumericalPoint uPoint(dimension);
      NumericalScalar productPDF(1.0);
      for (UnsignedLong i = 0; i < dimension; ++i)
        {
          const NumericalPoint component(NumericalPoint(1, point[i]));
          uPoint[i] = distributionCollection_[i].computeCDF(component);
          productPDF *= distributionCollection_[i].computePDF(component);
        }
      // Special case for the independent case, to boost performances
      if (hasIndependentCopula()) return productPDF;
      // General case
      return copula_.computePDF(uPoint) * productPDF;
      }

      /* Get the CDF of the ComposedDistribution */
      NumericalScalar ComposedDistribution::computeCDF(const NumericalPoint & point, const Bool tail) const
      {
      /* CDF = CDF_copula(CDF_dist1(p1), ..., CDF_distn(pn)) */
      const UnsignedLong dimension(getDimension());
      if(point.getDimension() != dimension) throw InvalidArgumentException(HERE) << "The given point has a wrong dimension";
      NumericalPoint uPoint(dimension);
      for (UnsignedLong i = 0; i < dimension; ++i)
        {
          uPoint[i] = distributionCollection_[i].computeCDF(NumericalPoint(1, point[i]), tail);
        }
00307       // Special case for the independent case, to boost performances
      if (hasIndependentCopula()) return IndependentCopula(dimension).computeCDF(uPoint);
      // General case
      return copula_.computeCDF(uPoint, tail);
      }

      /* Compute the probability content of an interval */
      NumericalScalar ComposedDistribution::computeProbability(const Interval & interval) const
      {
      // If the interval is empty
      if (interval.isEmpty()) return 0.0;
      const UnsignedLong dimension(getDimension());
      const NumericalPoint lower(interval.getLowerBound());
      const NumericalPoint upper(interval.getUpperBound());
      const Interval::BoolCollection finiteLower(interval.getFiniteLowerBound());
      const Interval::BoolCollection finiteUpper(interval.getFiniteUpperBound());
      NumericalPoint lowerCopula(dimension);
      NumericalPoint upperCopula(dimension);
      for (UnsignedLong i = 0; i < dimension; ++i)
        {
          if (finiteLower[i])
            {
            lowerCopula[i] = distributionCollection_[i].computeCDF(lower[i]);
            }
          else
            {
            lowerCopula[i] = 0.0;
            }
          if (finiteUpper[i])
            {
            upperCopula[i] = distributionCollection_[i].computeCDF(upper[i]);
            }
          else
            {
00341             upperCopula[i] = 1.0;
            }
        }
      return copula_.computeProbability(Interval(lowerCopula, upperCopula));
      }

      /* Get the PDF gradient of the distribution */
      ComposedDistribution::NumericalPoint ComposedDistribution::computePDFGradient(const NumericalPoint & point) const
      {
      const UnsignedLong dimension(getDimension());
      NumericalPoint gradient;
      // First, put the gradient according to marginal parameters
      // The marginal parameters are supposed to be independent from one marginal distribution
      // to the others
      for (UnsignedLong i = 0; i < dimension; ++i)
        {
          const NumericalPoint marginalGradient(distributionCollection_[i].computePDFGradient(point));
          const UnsignedLong marginalParameterDimension(marginalGradient.getDimension());
          for (UnsignedLong j = 0; j < marginalParameterDimension; ++j)
            {
            gradient.add(marginalGradient[j]);
            }
        }
      const NumericalPoint copulaGradient(copula_.computePDFGradient(point));
      const UnsignedLong copulaParameterDimension(copulaGradient.getDimension());
      // Then, put the gradient according to the copula parameters
      for (UnsignedLong j = 0; j < copulaParameterDimension; ++j)
00368         {
          gradient.add(copulaGradient[j]);
        }
      return gradient;
      }

      /* Get the CDF gradient of the distribution */
      ComposedDistribution::NumericalPoint ComposedDistribution::computeCDFGradient(const NumericalPoint & point) const
      {
      const UnsignedLong dimension(getDimension());
      NumericalPoint gradient;
      // First, put the gradient according to marginal parameters
      // The marginal parameters are supposed to be independent from one marginal distribution
      // to the others
      for (UnsignedLong i = 0; i < dimension; ++i)
        {
          const NumericalPoint marginalGradient(distributionCollection_[i].computeCDFGradient(point));
          const UnsignedLong marginalParameterDimension(marginalGradient.getDimension());
          for (UnsignedLong j = 0; j < marginalParameterDimension; ++j)
            {
            gradient.add(marginalGradient[j]);
            }
        }
      const NumericalPoint copulaGradient(copula_.computeCDFGradient(point));
      const UnsignedLong copulaParameterDimension(copulaGradient.getDimension());
      // Then, put the gradient according to the copula parameters
      for (UnsignedLong j = 0; j < copulaParameterDimension; ++j)
00395         {
          gradient.add(copulaGradient[j]);
        }
      return gradient;
      }

      /* Get the quantile of the ComposedDistribution */
      ComposedDistribution::NumericalPoint ComposedDistribution::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]";
      NumericalPoint quantile(copula_.computeQuantile(prob));
      for (UnsignedLong i = 0; i < getDimension(); ++i)
00407         {
          quantile[i] = distributionCollection_[i].computeQuantile(quantile[i])[0];
        }
      return quantile;
      }

      /* Compute the PDF of Xi | X1, ..., Xi-1. x = Xi, y = (X1,...,Xi-1) */
      NumericalScalar ComposedDistribution::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
      NumericalPoint u(conditioningDimension);
      for (UnsignedLong i = 0; i < conditioningDimension; ++i)
00423         {
          u[i] = distributionCollection_[i].computeCDF(y[i]);
        }
      return distributionCollection_[conditioningDimension].computePDF(x) * copula_.computeConditionalPDF(distributionCollection_[conditioningDimension].computeCDF(x), u);
      }

      /* Compute the CDF of Xi | X1, ..., Xi-1. x = Xi, y = (X1,...,Xi-1) */
      NumericalScalar ComposedDistribution::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 distributionCollection_[conditioningDimension].computeCDF(x);
      // General case
      NumericalPoint u(conditioningDimension);
      for (UnsignedLong i = 0; i < conditioningDimension; ++i)
00439         {
          u[i] = distributionCollection_[i].computeCDF(y[i]);
        }
      return copula_.computeConditionalCDF(distributionCollection_[conditioningDimension].computeCDF(x), u);
      }

      /* Get the mean of the distribution. It is cheap if the marginal means are cheap */
      ComposedDistribution::NumericalPoint ComposedDistribution::getMean() const throw(NotDefinedException)
      {
      const UnsignedLong dimension(getDimension());
      NumericalPoint mean(dimension);
      for (UnsignedLong i = 0; i < dimension; ++i)
00451         {
          mean[i] = distributionCollection_[i].getMean()[0];
        }
      return mean;
      }

      /* Get the standard deviation of the distribution */
      ComposedDistribution::NumericalPoint ComposedDistribution::getStandardDeviation() const throw(NotDefinedException)
      {
      const UnsignedLong dimension(getDimension());
      NumericalPoint standardDeviation(dimension);
      for (UnsignedLong i = 0; i < dimension; ++i)
00463         {
          standardDeviation[i] = distributionCollection_[i].getStandardDeviation()[0];
        }
      return standardDeviation;
      }

      /* Compute the covariance of the distribution */
      void ComposedDistribution::computeCovariance() const throw(NotDefinedException)
      {
      const 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();
      // Compute the marginal standard deviation
      const NumericalPoint sigma(getStandardDeviation());
      // Compute the weights and nodes of the 1D gauss quadrature over [-1, 1]
      NumericalSample nodesAndWeights(getGaussNodesAndWeights());
      // Convert the nodes and weights for the interval [0, 1]
      for (UnsignedLong i = 0; i < integrationNodesNumber_; ++i)
        {
          nodesAndWeights[0][i] = 0.5 * (nodesAndWeights[0][i] + 1.0);
          nodesAndWeights[1][i] *= 0.5;
        }
      // Compute the marginal quantiles at the nodes
      NumericalSample marginalQuantiles(integrationNodesNumber_, dimension);
      NumericalSample marginalPDF(integrationNodesNumber_, dimension);
      for(UnsignedLong component = 0; component < dimension; ++component)
        {
          Distribution marginalDistribution(getMarginal(component));
          for(UnsignedLong nodeIndex = 0; nodeIndex < integrationNodesNumber_; ++nodeIndex)
            {
            const NumericalScalar node(nodesAndWeights[0][nodeIndex]);
            const NumericalPoint q(marginalDistribution.computeQuantile(node));
            marginalQuantiles[nodeIndex][component] = q[0];
            marginalPDF[nodeIndex][component] = marginalDistribution.computePDF(q);
            }
          covariance_(component, component) = sigma[component] * sigma[component];
        }
      // Off-diagonal terms if the copula is not the independent copula
      if (!hasIndependentCopula())
        {
          // Performs the integration for each covariance in the strictly lower triangle of the covariance matrix
          // We simply use a product gauss quadrature
          // We first loop over the coeeficients because the most expensive task is to get the 2D marginal copulas
            Indices indices(2);
          for(UnsignedLong rowIndex = 0; rowIndex < dimension; ++rowIndex)
            {
            indices[0] = rowIndex;
            // We must fill the upper triangle of the covariance matrix in order to access the 2D marginal distributions
            // of the copula in the correct order for the ComposedCopula
            for(UnsignedLong columnIndex = rowIndex + 1; columnIndex < dimension; ++columnIndex)
              {
                indices[1] = columnIndex;
                const Distribution marginalCopula(copula_.getMarginal(indices));
                if (!marginalCopula.getImplementation()->hasIndependentCopula())
                  {
                  NumericalScalar covarianceIJ(0.0);
                  // Then we loop over the integration points
                  for(UnsignedLong rowNodeIndex = 0; rowNodeIndex < integrationNodesNumber_; ++rowNodeIndex)
                    {
                      const NumericalScalar nodeI(nodesAndWeights[0][rowNodeIndex]);
                      const NumericalScalar weightI(nodesAndWeights[1][rowNodeIndex]);
                      for(UnsignedLong columnNodeIndex = 0; columnNodeIndex < integrationNodesNumber_; ++columnNodeIndex)
                        {
                        const NumericalScalar nodeJ(nodesAndWeights[0][columnNodeIndex]);
                        const NumericalScalar weightJ(nodesAndWeights[1][columnNodeIndex]);
                        NumericalPoint in(2);
                        in[0] = nodeI;
                        in[1] = nodeJ;
                        covarianceIJ += weightI * weightJ * (marginalQuantiles[rowNodeIndex][rowIndex] - mean_[rowIndex]) * 
                          (marginalQuantiles[columnNodeIndex][columnIndex] - mean_[columnIndex]) * marginalCopula.computePDF(in);
                        } // loop over J integration nodes
                    } // loop over I integration nodes
                  covariance_(rowIndex, columnIndex) = covarianceIJ;
                  }
00542               } // loop over column indices
            } // loop over row indices
        } // if !hasIndependentCopula
      isAlreadyComputedCovariance_ = true;
      } // computeCovariance

      /* Get the skewness of the distribution */
      ComposedDistribution::NumericalPoint ComposedDistribution::getSkewness() const throw(NotDefinedException)
      {
      const UnsignedLong dimension(getDimension());
      NumericalPoint skewness(getDimension());
      for (UnsignedLong i = 0; i < dimension; ++i)
00554         {
          skewness[i] = distributionCollection_[i].getSkewness()[0];
        }
      return skewness;
      }

      /* Get the kurtosis of the distribution */
      ComposedDistribution::NumericalPoint ComposedDistribution::getKurtosis() const throw(NotDefinedException)
      {
      const UnsignedLong dimension(getDimension());
      NumericalPoint kurtosis(dimension);
      for (UnsignedLong i = 0; i < dimension; ++i)
00566         {
          kurtosis[i] = distributionCollection_[i].getKurtosis()[0];
        }
      return kurtosis;
      }

      /* Get the i-th marginal distribution */
00573       ComposedDistribution::Implementation ComposedDistribution::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]";
      return distributionCollection_[i].getImplementation();
      }

      /* Get the distribution of the marginal distribution corresponding to indices dimensions */
      ComposedDistribution::Implementation ComposedDistribution::getMarginal(const Indices & indices) const throw(InvalidArgumentException)
      {
      // This call will check that indices are correct
      const Copula marginalCopula(copula_.getMarginal(indices));
      DistributionCollection marginalDistributions(0);
      const UnsignedLong size(indices.getSize());
      for (UnsignedLong i = 0; i < size; ++i)
00587         {
          marginalDistributions.add(distributionCollection_[indices[i]]);
        }
      return new ComposedDistribution(marginalDistributions, marginalCopula);
      }

      /* Get the isoprobabilist transformation */
      ComposedDistribution::IsoProbabilisticTransformation ComposedDistribution::getIsoProbabilisticTransformation() const
      {
      // Set the parameters values and descriptions
      NumericalPointWithDescriptionCollection parametersCollection(getParametersCollection());
      // First, compute the dimension of the marginal parameters space
      UnsignedLong dimension(0);
      const UnsignedLong size(parametersCollection.getSize());
      for (UnsignedLong i = 0; i < size; ++i)
        {
          dimension += parametersCollection[i].getDimension();
        }
      NumericalPointWithDescription parameters(dimension);
      Description description(dimension);
      UnsignedLong index(0);
      for (UnsignedLong i = 0; i < size; ++i)
        {
          const NumericalPointWithDescription marginalParameters(parametersCollection[i]);
          const Description marginalDescription(marginalParameters.getDescription());
          const UnsignedLong marginalDimension(marginalParameters.getDimension());
          const String marginalName(marginalParameters.getName());
          for (UnsignedLong j = 0; j < marginalDimension; ++j)
            {
            parameters[index] = marginalParameters[j];
            description[index] = OSS() << marginalName << "_" << marginalDescription[j];
            ++index;
            }
        }
      parameters.setDescription(description);
      // Get the IsoProbabilisticTransformation from the copula
      const IsoProbabilisticTransformation isoprobabilistic(copula_.getIsoProbabilisticTransformation());
      // Get the right function implementations
      const EvaluationImplementation p_rightFunction = new MarginalTransformationEvaluation(distributionCollection_);

      // Get the right gradient implementations
      const GradientImplementation p_rightGradient = new MarginalTransformationGradient(distributionCollection_);

      // Get the right hessian implementations
      const HessianImplementation p_rightHessian = new MarginalTransformationHessian(distributionCollection_);

      IsoProbabilisticTransformation right(p_rightFunction, p_rightGradient, p_rightHessian);
00634       right.setParameters(parameters);
      IsoProbabilisticTransformation transformation(isoprobabilistic, right);

      return transformation;
      }

      /* Get the inverse isoprobabilist transformation */
      ComposedDistribution::InverseIsoProbabilisticTransformation ComposedDistribution::getInverseIsoProbabilisticTransformation() const
      {
      // Set the parameters values and descriptions
      NumericalPointWithDescriptionCollection parametersCollection(getParametersCollection());
      // First, compute the dimension of the marginal parameters space
      UnsignedLong dimension(0);
      const UnsignedLong size(parametersCollection.getSize());
      for (UnsignedLong i = 0; i < size; ++i)
        {
          dimension += parametersCollection[i].getDimension();
        }
      NumericalPointWithDescription parameters(dimension);
      Description description(dimension);
      UnsignedLong index(0);
      for (UnsignedLong i = 0; i < size; ++i)
        {
          const NumericalPointWithDescription marginalParameters(parametersCollection[i]);
          const Description marginalDescription(marginalParameters.getDescription());
          const UnsignedLong marginalDimension(marginalParameters.getDimension());
          const String marginalName(marginalParameters.getName());
          for (UnsignedLong j = 0; j < marginalDimension; ++j)
            {
            parameters[index] = marginalParameters[j];
            description[index] = OSS() << marginalName << "_" << marginalDescription[j];
            ++index;
            }
        }
      parameters.setDescription(description);
      // Get the IsoProbabilisticTransformation from the copula
      const InverseIsoProbabilisticTransformation inverseIsoprobabilistic(copula_.getInverseIsoProbabilisticTransformation());
      // Get the left and right function implementations
      const EvaluationImplementation p_leftFunction = new InverseMarginalTransformationEvaluation(distributionCollection_);

      // Get the left and right gradient implementations
      const GradientImplementation p_leftGradient = new InverseMarginalTransformationGradient(distributionCollection_);

      // Get the left and right hessian implementations
      const HessianImplementation p_leftHessian = new InverseMarginalTransformationHessian(distributionCollection_);

      InverseIsoProbabilisticTransformation left(p_leftFunction, p_leftGradient, p_leftHessian);
      left.setParameters(parameters);
00682 
      const InverseIsoProbabilisticTransformation transformation(left, inverseIsoprobabilistic);

      return transformation;
      }

00688       /* Get the standard distribution */
      ComposedDistribution::Implementation ComposedDistribution::getStandardDistribution() const
      {
      return copula_.getStandardDistribution().getImplementation();
      }

      /* Parameters value and description accessor */
      ComposedDistribution::NumericalPointWithDescriptionCollection ComposedDistribution::getParametersCollection() const
      {
      const UnsignedLong dimension(getDimension());
      NumericalPointWithDescriptionCollection parameters(dimension + (dimension > 1 ? 1 : 0));
      const Description description(getDescription());
      // First put the marginal parameters
      for (UnsignedLong marginalIndex = 0; marginalIndex < dimension; ++marginalIndex)
        {
          // Each marginal distribution must output a collection of parameters of size 1, even if it contains an empty NumericalPoint
          const NumericalPointWithDescriptionCollection marginalParameters(distributionCollection_[marginalIndex].getParametersCollection());
          NumericalPointWithDescription point(marginalParameters[0]);
          point.setName(description[marginalIndex]);
          parameters[marginalIndex] = point;
        } // marginalIndex
#ifdef WITH_DEPENDENCE_SENSITIVITY
      if (dimension > 1)
        {
          // Second put the dependence parameters
          const NumericalPointWithDescriptionCollection copulaParameters(copula_.getParametersCollection());
          NumericalPointWithDescription point(0);
          if (copulaParameters.getSize() != 0)
            {
            point = copulaParameters[0];
            }
          point.setName(copula_.getName());
00720           parameters[dimension] = point;
        } // dimension > 1
#endif
      return parameters;
      } // getParametersCollection

00726       /* Tell if the distribution has independent copula */
      Bool ComposedDistribution::hasIndependentCopula() const
      {
      return copula_.getImplementation()->hasIndependentCopula();
      }

      /* Method save() stores the object through the StorageManager */
      void ComposedDistribution::save(const StorageManager::Advocate & adv) const
00734       {
      DistributionImplementation::save(adv);
      adv.writeValue(distributionCollection_, StorageManager::MemberNameAttribute, "distributionCollection_");
      adv.writeValue(copula_, StorageManager::MemberNameAttribute, "copula_");
      }

      /* Method load() reloads the object from the StorageManager */
      void ComposedDistribution::load(const StorageManager::Advocate & adv)
      {
      DistributionImplementation::load(adv);
      adv.readValue(distributionCollection_, StorageManager::MemberNameAttribute, "distributionCollection_");
      adv.readValue(copula_, StorageManager::MemberNameAttribute, "copula_");
      // To compute the range
      setDistributionCollection(distributionCollection_);
      }
      
    } /* namespace Distribution */
  } /* namespace Uncertainty */
} /* namespace OpenTURNS */

Generated by  Doxygen 1.6.0   Back to index