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

RandomMixture.cxx

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

namespace OpenTURNS {

  namespace Base
  {  
    namespace Type
    {
      using Uncertainty::Distribution::RandomMixture;
      TEMPLATE_CLASSNAMEINIT(PersistentCollection<NumericalComplex>);
      static Common::Factory<PersistentCollection<NumericalComplex> > RegisteredFactory2("PersistentCollection<NumericalComplex>");
      using Uncertainty::Model::Distribution;
      TEMPLATE_CLASSNAMEINIT(PersistentCollection<Distribution>);
      static Common::Factory<PersistentCollection<Distribution> > RegisteredFactory3("PersistentCollection<Distribution>");
    } /* namespace Type */
  } /* namespace Base */

  namespace Uncertainty {

    namespace Distribution {

      typedef Base::Common::NotYetImplementedException NotYetImplementedException;
      typedef Base::Common::InternalException          InternalException;

      CLASSNAMEINIT(RandomMixture);

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

      const UnsignedLong RandomMixture::DefaultBlockMin;
      const UnsignedLong RandomMixture::DefaultBlockMax;
      const UnsignedLong RandomMixture::DefaultMaxSize;
      const NumericalScalar RandomMixture::DefaultAlpha;
      const NumericalScalar RandomMixture::DefaultPDFEpsilon;
      const NumericalScalar RandomMixture::DefaultCDFEpsilon;
      const NumericalScalar RandomMixture::GraphPDFEpsilon;
      const NumericalScalar RandomMixture::GraphCDFEpsilon;
00065 
      /* Default constructor */
      RandomMixture::RandomMixture(const DistributionCollection & coll, const NumericalScalar constant)
      throw (InvalidArgumentException)
      : DistributionImplementation("RandomMixture"),
        distributionCollection_(),
        constant_(constant),
        blockMin_(DefaultBlockMin),
        blockMax_(DefaultBlockMax),
        maxSize_(DefaultMaxSize),
        storedSize_(0),
        characteristicValuesCache_(0),
        alpha_(DefaultAlpha),
        beta_(0.0),
        pdfPrecision_(DefaultPDFEpsilon),
        cdfPrecision_(DefaultCDFEpsilon)
      {
      // We could NOT set distributionCollection_ in the member area of the constructor
      // because we must check before if the collection is valid (ie, if all the
      // distributions of the collection have the same dimension). We do this by calling
      // the setDistributionCollection() method that do it for us.
      setDistributionCollection( coll );
      computePositionIndicator();
      computeDispersionIndicator();
      computeBeta();
      computeReferenceBandwidth();
      }
00092 
      /* Compute the numerical range of the distribution given the parameters values */
      void RandomMixture::computeRange()
      {
      const UnsignedLong size(distributionCollection_.getSize());
      Interval range(constant_, constant_);
      for (UnsignedLong i = 0; i < size; ++i) range += distributionCollection_[i].getRange() * distributionCollection_[i].getWeight();
      setRange(range);
      }
00101 
      /* Comparison operator */
      Bool RandomMixture::operator ==(const RandomMixture & other) const {
      Bool sameObject = false;

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

      return sameObject;
      }
00112   
      /* String converter */
      String RandomMixture::str() const {
      OSS oss;
      oss << "class=" << RandomMixture::GetClassName()
          << " name=" << getName()
          << " distribution collection=" << distributionCollection_
          << " constant=" << constant_;
      return oss;
      }
00122   
      /* Weights distribution accessor */
      void RandomMixture::setWeights(const NumericalPoint & weights)
      {
      if (weights.getDimension() != distributionCollection_.getSize()) throw InvalidArgumentException(HERE) << "Error: the weights collection must have the same size as the distribution collection";
      const UnsignedLong size(distributionCollection_.getSize());
      for (UnsignedLong i = 0; i < size; ++i)
        {
          distributionCollection_[i].setWeight(weights[i]);
        }
      }

      RandomMixture::NumericalPoint RandomMixture::getWeights() const
      {
      const UnsignedLong size(distributionCollection_.getSize());
      NumericalPoint weights(size);
      for (UnsignedLong i = 0; i < size; ++i)
        {
          weights[i] = distributionCollection_[i].getWeight();
        }
      return weights;
      }

00145       
      /* Distribution collection accessor */
      void RandomMixture::setDistributionCollection(const DistributionCollection & coll)
      throw (InvalidArgumentException)
      {
      distributionCollection_ = coll;
      const UnsignedLong size(coll.getSize());
      for(UnsignedLong i = 0; i < size; ++i)
        {
          if (coll[i].getDimension() != 1) throw InvalidArgumentException(HERE) << "Error: a RandomMixture cannot be built from a collection of distributions of dimension not equal to 1";
        }
      setDimension(1);
      }
00158 
      /* Constant accessor */
      void RandomMixture::setConstant(const NumericalScalar constant)
      {
      constant_ = constant;
      }

      NumericalScalar RandomMixture::getConstant() const
      {
      return constant_;
      }

      /* Distribution collection accessor */
      const RandomMixture::DistributionCollection & RandomMixture::getDistributionCollection() const
      {
      return distributionCollection_;
      }
00175 
      /* Virtual constructor */
      RandomMixture * RandomMixture::clone() const
      {
      return new RandomMixture(*this);
      }
00181 
      /* Get one realization of the RandomMixture */
      RandomMixture::NumericalPoint RandomMixture::getRealization() const
      {
      const UnsignedLong size(distributionCollection_.getSize());
      NumericalPoint realization(1, constant_);
      for (UnsignedLong i = 0; i < size; ++i)
        {
          realization += distributionCollection_[i].getWeight() * distributionCollection_[i].getRealization();
        }
      return realization;
      }
     

      /* Get the DDF of the RandomMixture */
      RandomMixture::NumericalPoint RandomMixture::computeDDF(const NumericalPoint & point) const 
      {
      throw NotYetImplementedException(HERE);
      }

      /* Get the PDF of the RandomMixture. It uses the Poisson inversion formula as described in the reference:
       "Abate, J. and Whitt, W. (1992). The Fourier-series method for inverting
       transforms of probability distributions. Queueing Systems 10, 5--88., 1992",
       formula 5.5.
       We use an incremental update of the trigonometric functions and reduce the complex arithmetic to a real
00206        arithmetic for performance purpose.
      */
      NumericalScalar RandomMixture::computePDF(const NumericalPoint & point) const 
      {
      const UnsignedLong size(distributionCollection_.getSize());
      const NumericalScalar x(point[0]);
      // Special case for combination containing only one contributor
      if (size == 1) return distributionCollection_[0].computePDF((x - constant_) / distributionCollection_[0].getWeight());
      // Check range
      const Interval range(getRange());
      const NumericalPoint lowerBound(range.getLowerBound());
      const NumericalPoint upperBound(range.getUpperBound());
      if ((x <= lowerBound[0]) || (x >= upperBound[0])) return 0.0;
      const NumericalScalar position(getPositionIndicator());
      const NumericalScalar dispersion(getDispersionIndicator());
      // TODO: use range information
      const NumericalScalar bandwidth(adjustBandwidth(2.0 * M_PI / (2.0 * fabs(x - position) + beta_ * dispersion)));
      const NumericalScalar cosHX(cos(bandwidth * x));
      const NumericalScalar sinHX(sin(bandwidth * x));
      NumericalComplex cfValue(computeCharacteristicFunction(1, bandwidth));
      NumericalScalar value(0.5 + cfValue.real() * cosHX + cfValue.imag() * sinHX);
      NumericalScalar cosMHX(cosHX);
      NumericalScalar sinMHX(sinHX);
      UnsignedLong k(1);
      const NumericalScalar precision(pdfPrecision_ * M_PI / bandwidth);
      const UnsignedLong kmin(1 << blockMin_);
      const UnsignedLong kmax(1 << blockMax_);
      NumericalScalar error(2.0 * precision);
      while ( (k < kmin) || ( (k < kmax) && (error > precision) ) )
        {
          error = 0.0;
          for (UnsignedLong m = k + 1; m <= 2 * k; ++m)
            {
            const NumericalScalar tmp(cosMHX);
            cosMHX = tmp * cosHX - sinMHX * sinHX;
            sinMHX = sinMHX * cosHX + tmp * sinHX;
            cfValue = computeCharacteristicFunction(m, bandwidth);
            error += cfValue.real() * cosMHX + cfValue.imag() * sinMHX;
            }
          value += error;
          error = fabs(error);
          k *= 2;
        }
      value *= bandwidth / M_PI;
      // For very low level of PDF, the computed value can be slightly negative. Round it up to zero.
      if (value < 0.0) value = 0.0;
      pdfEpsilon_ = error * bandwidth / M_PI;
      return value;
      }
00255 
      /*  Compute the PDF of 1D distributions over a regular grid. The precision is reduced as this method is for drawing purpose only. */
      RandomMixture::NumericalSample RandomMixture::computePDF(const NumericalScalar xMin, const NumericalScalar xMax, const UnsignedLong pointNumber) const
      {
      if (getDimension() != 1) throw InvalidArgumentException(HERE) << "Error: cannot compute the PDF over a regular 1D grid if the dimension is > 1";
      NumericalSample result(pointNumber, 2);
      NumericalScalar x(xMin);
      const NumericalScalar step((xMax - xMin) / NumericalScalar(pointNumber - 1.0));
      const NumericalScalar oldPDFPrecision(pdfPrecision_);
      pdfPrecision_ = GraphPDFEpsilon;
      for (UnsignedLong i = 0; i < pointNumber; ++i)
        {
          result[i][0] = x;
          result[i][1] = computePDF(x);
          x += step;
        }
      pdfPrecision_ = oldPDFPrecision;
      return result;
      }
00274 
      /* Get the CDF of the RandomMixture */
      NumericalScalar RandomMixture::computeCDF(const NumericalPoint & point, const Bool tail) const 
      {
      const UnsignedLong size(distributionCollection_.getSize());
      const NumericalScalar x(point[0]);
      // Special case for combination containing only one contributor
      if (size == 1) return distributionCollection_[0].computeCDF((x - constant_) / distributionCollection_[0].getWeight(), tail);
      // Check range
      const Interval range(getRange());
      const NumericalPoint lowerBound(range.getLowerBound());
      const NumericalPoint upperBound(range.getUpperBound());
      if (x <= lowerBound[0]) return 0.0;
      if (x >= upperBound[0]) return 1.0;
      // Here we call computeProbability with a ]-inf, x] interval
      const Interval interval(lowerBound, point, getRange().getFiniteLowerBound(), Interval::BoolCollection(1, true));
      return computeProbability(interval);
      }
00292 
      /*  Compute the CDF of 1D distributions over a regular grid. The precision is reduced as this method is for drawing purpose only. */
      RandomMixture::NumericalSample RandomMixture::computeCDF(const NumericalScalar xMin, const NumericalScalar xMax, const UnsignedLong pointNumber) const
      {
      if (getDimension() != 1) throw InvalidArgumentException(HERE) << "Error: cannot compute the CDF over a regular 1D grid if the dimension is > 1";
      NumericalSample result(pointNumber, 2);
      NumericalScalar x(xMin);
      const NumericalScalar step((xMax - xMin) / NumericalScalar(pointNumber - 1.0));
      const NumericalScalar oldCDFPrecision(cdfPrecision_);
      cdfPrecision_ = GraphPDFEpsilon;
      for (UnsignedLong i = 0; i < pointNumber; ++i)
        {
          result[i][0] = x;
          result[i][1] = computeCDF(x);
          x += step;
        }
      cdfPrecision_ = oldCDFPrecision;
      return result;
      }

      /* Get the probability content of an interval. It uses the Poisson inversion formula as described in the reference:
       "Abate, J. and Whitt, W. (1992). The Fourier-series method for inverting
       transforms of probability distributions. Queueing Systems 10, 5--88., 1992",
       formula 5.14.
       We use an incremental update of the trigonometric functions and reduce the complex arithmetic to a real
00317        arithmetic for performance purpose.
      */
      NumericalScalar RandomMixture::computeProbability(const Interval & interval) const 
      {
      if (interval.isEmpty()) return 0.0;
      const UnsignedLong size(distributionCollection_.getSize());
      // Special case for combination containing only one contributor
      if (size == 1)
        {
          const NumericalScalar lower(interval.getLowerBound()[0]);
          const NumericalScalar upper(interval.getUpperBound()[0]);
          const NumericalScalar weight(getWeight());
          // Negative weight, swap upper and lower bound flags
          if (weight < 0.0)
            {
            const Interval adjustedInterval(NumericalPoint(1, (upper - constant_) / weight), NumericalPoint(1, (lower - constant_) / weight), interval.getFiniteUpperBound(), interval.getFiniteUpperBound());
            return distributionCollection_[0].computeProbability(adjustedInterval.intersect(getRange()));
            }
          else
            {
            const Interval adjustedInterval(NumericalPoint(1, (lower - constant_) / weight), NumericalPoint(1, (upper - constant_) / weight), interval.getFiniteUpperBound(), interval.getFiniteUpperBound());
            return distributionCollection_[0].computeProbability(interval.intersect(getRange()));
            }
        }
      const Bool finiteLowerBound(interval.getFiniteLowerBound()[0]);
      const Bool finiteUpperBound(interval.getFiniteUpperBound()[0]);
      // Quick return for integral over the whole real line
      if (!finiteLowerBound && !finiteUpperBound) return 1.0;
      NumericalScalar lowerBound(interval.getLowerBound()[0]);
        NumericalScalar upperBound(interval.getUpperBound()[0]);
      const NumericalScalar position(getPositionIndicator());
      const NumericalScalar dispersion(getDispersionIndicator());
      // TODO: use range information
      const NumericalScalar bandwidth(adjustBandwidth(M_PI / (beta_ * dispersion)));
      // First, truncate infinite intervals
      if (!finiteLowerBound) lowerBound = position - beta_ * dispersion;
      if (!finiteUpperBound) upperBound = position + beta_ * dispersion;
      NumericalComplex cfValue(computeCharacteristicFunction(1, bandwidth));
      const NumericalScalar precision(cdfPrecision_ * M_PI / bandwidth);
      NumericalScalar error(2.0 * precision);
      NumericalScalar mH(bandwidth);
      const NumericalScalar cosHUpper(cos(bandwidth * upperBound));
      const NumericalScalar sinHUpper(sin(bandwidth * upperBound));
      const NumericalScalar cosHLower(cos(bandwidth * lowerBound));
      const NumericalScalar sinHLower(sin(bandwidth * lowerBound));
      NumericalScalar cosMHUpper(cosHUpper);
      NumericalScalar sinMHUpper(sinHUpper);
      NumericalScalar cosMHLower(cosHLower);
      NumericalScalar sinMHLower(sinHLower);
      NumericalScalar value(0.5 * (upperBound - lowerBound) + ( cfValue.real() * (-sinMHLower + sinMHUpper) + cfValue.imag() * (cosMHLower - cosMHUpper) ) / mH);
      UnsignedLong k(1);
      const UnsignedLong kmin(1 << blockMin_);
      const UnsignedLong kmax(1 << blockMax_);
      while ( (k < kmax) && (error > precision || k < kmin) )
        {
          error = 0.0;
          for (UnsignedLong m = k+1; m <= 2 * k; ++m)
            {
            mH += bandwidth;
            NumericalScalar tmp(cosMHUpper * cosHUpper - sinMHUpper * sinHUpper);
            sinMHUpper = sinMHUpper * cosHUpper + cosMHUpper * sinHUpper;
            cosMHUpper = tmp;
            tmp = cosMHLower * cosHLower - sinMHLower * sinHLower;
            sinMHLower = sinMHLower * cosHLower + cosMHLower * sinHLower;
            cosMHLower = tmp;
            cfValue = computeCharacteristicFunction(m, bandwidth);
            error += (cfValue.real() * (-sinMHLower + sinMHUpper) + cfValue.imag() * (cosMHLower - cosMHUpper)) / mH;
            }
          value += error;
          error = fabs(error);
          k *= 2;
        }
      value *= bandwidth / M_PI;
      // For extrem values of the argument, the computed value can be slightly outside of [0,1]. Truncate it.
      return (value < 0.0 ? 0.0 : (value > 1.0 ? 1.0 : value));
      }
00393 
       /* Get the characteristic function of the distribution, i.e. phi(u) = E(exp(I*u*X)) */
      NumericalComplex RandomMixture::computeCharacteristicFunction(const NumericalPoint & point) const
      {
      NumericalComplex cfValue(exp(NumericalComplex(0.0, constant_ * point[0])));
      UnsignedLong size(distributionCollection_.getSize());
      for(UnsignedLong i = 0; i < size; ++i)
        {
          cfValue *= distributionCollection_[i].computeCharacteristicFunction(distributionCollection_[i].getWeight() * point);
        } /* end for */
      return cfValue;
      }
00405 
      // Compute a value of the characteristic function on a prescribed discretization
      NumericalComplex RandomMixture::computeCharacteristicFunction(const UnsignedLong index, const NumericalScalar step) const
      {
      UnsignedLong realIndex = static_cast<UnsignedLong>(round(index * step / referenceBandwidth_ + 0.5));
      if (realIndex == 0) return 1.0;
      // The cached values are computed and stored in an ascending order without hole: this function is always called on sequence starting from 0 to n-1
      // Usual case first: the index is within the already computed values
      if (realIndex <= storedSize_) return characteristicValuesCache_[realIndex - 1];
      // If the index is higher than the maximum allowed storage
      if (realIndex > maxSize_) return computeCharacteristicFunction(NumericalPoint(1, realIndex * step));
      // Here, the index has not been computed so far, fill-in the gap
      if (realIndex > storedSize_)
        {
          for (UnsignedLong i = storedSize_; i < realIndex; ++i)
            {
            characteristicValuesCache_.add(computeCharacteristicFunction(NumericalPoint(1, i * referenceBandwidth_)));
            }
          storedSize_ = realIndex;
          return characteristicValuesCache_[storedSize_ - 1];
        }
      // Should never go there
      throw InternalException(HERE) << "Error: trying to access to a cached characteristic value in an incorrect pattern.";
      }
00429 
      /* Get the PDF gradient of the distribution */
      RandomMixture::NumericalPoint RandomMixture::computePDFGradient(const NumericalPoint & point) const
      {
      throw NotYetImplementedException(HERE);
      }
00435 
      /* Get the CDF gradient of the distribution */
      RandomMixture::NumericalPoint RandomMixture::computeCDFGradient(const NumericalPoint & point) const
      {
      throw NotYetImplementedException(HERE);
      }
00441 
      /* Compute the mean of the RandomMixture */
      void RandomMixture::computeMean() const throw(NotDefinedException)
      {
      mean_ = NumericalPoint(1, constant_);
      const UnsignedLong size(distributionCollection_.getSize());
      for(UnsignedLong i = 0; i < size; ++i)
        {
          mean_[0] += distributionCollection_[i].getWeight() * distributionCollection_[i].getMean()[0];
        } /* end for */
      isAlreadyComputedMean_ = true;
      }
00453 
      /* Compute the covariance of the RandomMixture */
      void RandomMixture::computeCovariance() const throw(NotDefinedException)
      {
      covariance_ = CovarianceMatrix(1);
      const UnsignedLong size(distributionCollection_.getSize());
      for(UnsignedLong i = 0; i < size; ++i)
        {
          covariance_(0, 0) += distributionCollection_[i].getWeight() * distributionCollection_[i].getCovariance().operator()(0, 0);
        }
      isAlreadyComputedCovariance_ = true;
      }
00465 
      /* Get the standard deviation of the RandomMixture */
      RandomMixture::NumericalPoint RandomMixture::getStandardDeviation() const throw(NotDefinedException)
      {
      return NumericalPoint(1, sqrt(getCovariance().operator()(0, 0)));
      }
00471 
      /* Get the skewness of the RandomMixture */
      RandomMixture::NumericalPoint RandomMixture::getSkewness() const throw(NotDefinedException)
      {
      NumericalPoint skewness(1, 0.0);
      const UnsignedLong size(distributionCollection_.getSize());
      for(UnsignedLong i = 0; i < size; ++i)
        {
          skewness[0] += pow(distributionCollection_[i].getWeight(), 3.0) * distributionCollection_[i].getSkewness()[0];
        } /* end for */
      skewness[0] /= pow(getStandardDeviation()[0], 3.0);
      return skewness;
      }
00484 
      /* Get the kurtosis of the RandomMixture */
      RandomMixture::NumericalPoint RandomMixture::getKurtosis() const throw(NotDefinedException)
      {
      NumericalPoint kurtosis(1, 0.0);
      NumericalScalar sum2(0.0);
      NumericalScalar sum4(0.0);
      const UnsignedLong size(distributionCollection_.getSize());
      for(UnsignedLong i = 0; i < size; ++i)
        {
          const NumericalScalar w(distributionCollection_[i].getWeight());
          const NumericalScalar w2(w * w);
          const NumericalScalar var(distributionCollection_[i].getCovariance().operator()(0, 0));
          const NumericalScalar term2(w2 * var);
          sum2 += term2;
          sum4 += term2 * term2;
          kurtosis[0] += w2 * w2 * distributionCollection_[i].getKurtosis()[0];
        } /* end for */
      kurtosis[0] = 3.0 + (kurtosis[0] - 3.0 * sum4) / (sum2 * sum2);
      return kurtosis;
      }
00505 
      /** Parameters value and description accessor */
      RandomMixture::NumericalPointWithDescriptionCollection RandomMixture::getParametersCollection() const
      {
      const UnsignedLong dimension(getDimension());
      const UnsignedLong size(distributionCollection_.getSize());
      // Special case for dimension=1
      if (dimension == 1)
        {
          NumericalPointWithDescriptionCollection parameters(1);
          Description description;
          // Form a big NumericalPoint from the parameters of each atom
          for (UnsignedLong i = 0; i < size; ++i)
            {
            const NumericalPointWithDescription atomParameters(distributionCollection_[i].getParametersCollection()[0]);
            const Description atomDescription(atomParameters.getDescription());
            const UnsignedLong atomParameterDimension(atomParameters.getDimension());
            // Add the current atom parameters
            for (UnsignedLong j = 0; j < atomParameterDimension; ++j)
              {
                parameters[0].add(atomParameters[i]);
                description.add(atomDescription[i]);
              }
            }
          parameters[0].setDescription(description);
          parameters[0].setName(getName());
          return parameters;
        }
      // General case
      NumericalPointWithDescriptionCollection parameters(dimension + 1);
      Description description;
      // 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(distributionCollection_[marginalIndex].getName());
          parameters[marginalIndex] = point;
        } // marginalIndex

      // Form a big NumericalPoint from the dependence parameters of each atom
      for (UnsignedLong i = 0; i < size; ++i)
        {
          const NumericalPointWithDescription atomDependenceParameters(distributionCollection_[i].getParametersCollection()[dimension]);
          const Description atomDescription(atomDependenceParameters.getDescription());
          const UnsignedLong atomParameterDimension(atomDependenceParameters.getDimension());
          // Add the current atom dependence parameters
          for (UnsignedLong j = 0; j < atomParameterDimension; ++j)
            {
            parameters[dimension].add(atomDependenceParameters[i]);
            description.add(atomDescription[i]);
            }
        }
      parameters[dimension].setDescription(description);
      parameters[dimension].setName("dependence");
      return parameters;
      } // getParametersCollection
00563       
      /* Get a positon indicator for a 1D distribution */
      NumericalScalar RandomMixture::getPositionIndicator() const
      {
      if (!isAlreadyComputedPositionIndicator_)
        {
          computePositionIndicator();
        }
      return positionIndicator_;
      }
00573 
      /* Compute a positon indicator for a 1D distribution */
      void RandomMixture::computePositionIndicator() const
      {
        positionIndicator_ = 0.0;
      const UnsignedLong size(distributionCollection_.getSize());
      // Assume an additive behaviour of the position indicator. It is true for the mean value, and almost true for the median of moderatly skewed distributions
      for(UnsignedLong i = 0; i < size; ++i)
        {
          positionIndicator_ += distributionCollection_[i].getWeight() * distributionCollection_[i].getPositionIndicator();
        } /* end for */
      isAlreadyComputedPositionIndicator_ = true;
      }
      
00587       
      /* Get a dispersion indicator for a 1D distribution */
      NumericalScalar RandomMixture::getDispersionIndicator() const
      {
      if (!isAlreadyComputedDispersionIndicator_)
        {
          computeDispersionIndicator();
        }
      return dispersionIndicator_;
      }
00597 
      /* Compute a dispersion indicator for a 1D distribution */
      void RandomMixture::computeDispersionIndicator() const
      {
      dispersionIndicator_ = 0.0;
      const UnsignedLong size(distributionCollection_.getSize());
      // Assume a quadratic additive behaviour of the dispersion indicator. It is true for the standard deviation value, and almost true for the interquartile of moderatly skewed distributions
      for(UnsignedLong i = 0; i < size; ++i)
        {
          dispersionIndicator_ += pow(distributionCollection_[i].getWeight() * distributionCollection_[i].getDispersionIndicator(), 2.0);
        } /* end for */
      dispersionIndicator_ = sqrt(dispersionIndicator_);
      isAlreadyComputedDispersionIndicator_ = true;
      }
00611 
      /* BlockMin accessor */
      void RandomMixture::setBlockMin(const UnsignedLong blockMin)
      {
      blockMin_ = blockMin;
      }

      UnsignedLong RandomMixture::getBlockMin() const
      {
      return blockMin_;
      }
00622 
      /* BlockMax accessor */
      void RandomMixture::setBlockMax(const UnsignedLong blockMax)
      {
      blockMax_ = blockMax;
      }
      
      UnsignedLong RandomMixture::getBlockMax() const
      {
      return blockMax_;
      }
00633 
      /* MaxSize accessor */
      void RandomMixture::setMaxSize(const UnsignedLong maxSize)
      {
      maxSize_ = maxSize;
      characteristicValuesCache_.resize(maxSize);
      }

      UnsignedLong RandomMixture::getMaxSize() const
      {
      return maxSize_;
      }
00645 
      /* Alpha accessor */
      void RandomMixture::setAlpha(const NumericalScalar alpha)
      {
      if (alpha <= 0.0) throw InvalidArgumentException(HERE) << "Error: the alpha parameter must be strictly positive";
      alpha_ = alpha;
      computeReferenceBandwidth();
      }

      NumericalScalar RandomMixture::getAlpha() const
      {
      return alpha_;
00657       }

      void RandomMixture::computeBeta()
      {
      const NumericalScalar step(getDispersionIndicator());
      NumericalScalar x(getPositionIndicator() + alpha_ * step);
      while (abs(computeCharacteristicFunction(x)) > pdfPrecision_)
        {
          x += step;
        }
      beta_ = 2.0 * (x - getPositionIndicator()) / getDispersionIndicator();
00668       }

      NumericalScalar RandomMixture::getBeta() const
      {
      return beta_;
      }

      /* Compute the reference bandwidth. It is defined as the maximum bandwidth
00676        that allow a precise computation of the PDF over the range
       [positionIndicator_ +/- beta * dispersionIndicator_] */
      void RandomMixture::computeReferenceBandwidth()
      {
      referenceBandwidth_ = 2.0 * M_PI / ((beta_ + 4.0 * alpha_) * getDispersionIndicator());
      // Reset the cached values
      storedSize_ = 0;
      characteristicValuesCache_ = NumericalComplexPersistentCollection(0);
      }

00686       /* Adjust a given bandwidth with respect to a reference bandwidth,
       in order to be an integral divisor or multiple of the reference bandwidth */
      NumericalScalar RandomMixture::adjustBandwidth(const NumericalScalar bandwidth) const
      {
      // If the given bandwidth is a strict divisor of the reference bandwidth
      if (bandwidth < referenceBandwidth_) return referenceBandwidth_ / floor(referenceBandwidth_ / bandwidth);
      // If the given bandwidth is a strict multiple of the reference bandwidth
      if (bandwidth > referenceBandwidth_) return floor(bandwidth / referenceBandwidth_) * referenceBandwidth_;
      // else the given bandwidth is equal to the reference bandwidth
      return referenceBandwidth_;
      }
00697 
      /* Method save() stores the object through the StorageManager */
      void RandomMixture::save(const StorageManager::Advocate & adv) const
      {
      DistributionImplementation::save(adv);
      adv.writeValue(distributionCollection_, StorageManager::MemberNameAttribute, "distributionCollection_");
      adv.writeValue("constant_", constant_);
      adv.writeValue("positionIndicator_", positionIndicator_);
      adv.writeValue("isAlreadyComputedPositionIndicator_", isAlreadyComputedPositionIndicator_);
      adv.writeValue("dispersionIndicator_", dispersionIndicator_);
      adv.writeValue("isAlreadyComputedDispersionIndicator_", isAlreadyComputedDispersionIndicator_);
      adv.writeValue("blockMin_", blockMin_);
      adv.writeValue("blockMax_", blockMax_);
      adv.writeValue("referenceBandwidth_", referenceBandwidth_);
      adv.writeValue("maxSize_", maxSize_);
      adv.writeValue("storedSize_", storedSize_);
      adv.writeValue(characteristicValuesCache_, StorageManager::MemberNameAttribute, "characteristicValuesCache_");
      adv.writeValue("alpha_", alpha_);
      adv.writeValue("beta_", beta_);
      adv.writeValue("pdfPrecision_", pdfPrecision_);
      adv.writeValue("cdfPrecision_", cdfPrecision_);
      } // save
00719 
      /* Method load() reloads the object from the StorageManager */
      void RandomMixture::load(const StorageManager::Advocate & adv)
      {
      DistributionImplementation::load(adv);
      adv.readValue(distributionCollection_, StorageManager::MemberNameAttribute, "distributionCollection_");
      adv.readValue(characteristicValuesCache_, StorageManager::MemberNameAttribute, "characteristicValuesCache_");
      // Read the Bool
      {
        String name;
        Bool value;
        StorageManager::List objList = adv.getList(StorageManager::BoolEntity);
        for(objList.firstValueToRead(); objList.moreValuesToRead(); objList.nextValueToRead()) {
          if (objList.readValue(name, value)) {
            if (name == "isAlreadyComputedPositionIndicator_") isAlreadyComputedPositionIndicator_ = value;
            if (name == "isAlreadyComputedDispersionIndicator_") isAlreadyComputedDispersionIndicator_ = value;
          }
        }
      } // Bool
      // Read the UnsignedLong
      {
        String name;
        UnsignedLong value;
        StorageManager::List objList = adv.getList(StorageManager::UnsignedLongEntity);
        for(objList.firstValueToRead(); objList.moreValuesToRead(); objList.nextValueToRead()) {
          if (objList.readValue(name, value)) {
            if (name == "blockMin_") blockMin_ = value;
            if (name == "blockMax_") blockMax_ = value;
            if (name == "maxSize_") maxSize_ = value;
            if (name == "storedSize_") storedSize_ = value;
          }
        }
      } // UnsignedLong
      // Read the NumericalScalar
      {
        String name;
        NumericalScalar value;
        StorageManager::List objList = adv.getList(StorageManager::NumericalScalarEntity);
        for(objList.firstValueToRead(); objList.moreValuesToRead(); objList.nextValueToRead()) {
          if (objList.readValue(name, value)) {
            if (name == "constant_") constant_ = value;
            if (name == "positionIndicator_") positionIndicator_ = value;
            if (name == "dispersionIndicator_") dispersionIndicator_ = value;
            if (name == "referenceBandwidth_") referenceBandwidth_ = value;
            if (name == "alpha_") alpha_ = value;
            if (name == "beta_") beta_ = value;
            if (name == "pdfPrecision_") pdfPrecision_ = value;
            if (name == "cdfPrecision_") cdfPrecision_ = value;
          }
        }
      } // NumericalScalar
      computeRange();
      } // load

    } /* namespace Distribution */
  } /* namespace Uncertainty */
} /* namespace OpenTURNS */

Generated by  Doxygen 1.6.0   Back to index