Improved DistributionFactory

Now the user can fix the value of an arbitrary subset of parameters and estimate the others by MLE in all the factories.

lib/src/Uncertainty/Distribution/AliMikhailHaqCopulaFactory.cxx

@@ -109,6 +109,7 @@ AliMikhailHaqCopula AliMikhailHaqCopulaFactory::buildAsAliMikhailHaqCopula(const
   }
   AliMikhailHaqCopula result(theta);
   result.setDescription(sample.getDescription());
+  adaptToKnownParameter(sample, &result);
   return result;
 }
 

lib/src/Uncertainty/Distribution/ArcsineFactory.cxx

@@ -76,6 +76,7 @@ Arcsine ArcsineFactory::buildAsArcsine(const Sample & sample) const
   parameters[1] = standardDeviation;
   Arcsine result(buildAsArcsine(ArcsineMuSigma()(parameters)));
   result.setDescription(sample.getDescription());
+  adaptToKnownParameter(sample, &result);
   return result;
 }
 

lib/src/Uncertainty/Distribution/BernoulliFactory.cxx

@@ -80,6 +80,7 @@ Bernoulli BernoulliFactory::buildAsBernoulli(const Sample & sample) const
   }
   Bernoulli result(sum / size);
   result.setDescription(sample.getDescription());
+  adaptToKnownParameter(sample, &result);
   return result;
 }
 

lib/src/Uncertainty/Distribution/BetaFactory.cxx

@@ -75,7 +75,10 @@ Beta BetaFactory::buildAsBeta(const Sample & sample) const
 
   const Scalar mu = sample.computeMean()[0];
   const Scalar sigma = sample.computeStandardDeviation()[0];
-  return buildAsBeta(BetaMuSigma(mu, sigma, a, b).evaluate());
+  Beta result(buildAsBeta(BetaMuSigma(mu, sigma, a, b).evaluate()));
+  result.setDescription(sample.getDescription());
+  adaptToKnownParameter(sample, &result);
+  return result;
 }
 
 Beta BetaFactory::buildAsBeta(const Point & parameters) const

lib/src/Uncertainty/Distribution/BinomialFactory.cxx

@@ -135,6 +135,7 @@ Binomial BinomialFactory::buildAsBinomial(const Sample & sample) const
   }
   Binomial result(maxN, mean / maxN);
   result.setDescription(sample.getDescription());
+  adaptToKnownParameter(sample, &result);
   return result;
 }
 

lib/src/Uncertainty/Distribution/BurrFactory.cxx

@@ -140,6 +140,7 @@ Burr BurrFactory::buildAsBurr(const Sample & sample) const
   const Scalar k = size / sumLogXC;
   Burr result(c, k);
   result.setDescription(sample.getDescription());
+  adaptToKnownParameter(sample, &result);
   return result;
 }
 

lib/src/Uncertainty/Distribution/ChiFactory.cxx

@@ -78,6 +78,7 @@ Chi ChiFactory::buildAsChi(const Sample & sample) const
   {
     Chi result(sumSquares / size);
     result.setDescription(sample.getDescription());
+    adaptToKnownParameter(sample, &result);
     return result;
   }
   catch (const InvalidArgumentException &)

lib/src/Uncertainty/Distribution/ChiSquareFactory.cxx

@@ -75,6 +75,7 @@ ChiSquare ChiSquareFactory::buildAsChiSquare(const Sample & sample) const
   if (xMin == xMax) throw InvalidArgumentException(HERE) << "Error: cannot estimate a ChiSquare distribution from a constant sample.";
   ChiSquare result(mean);
   result.setDescription(sample.getDescription());
+  adaptToKnownParameter(sample, &result);
   return result;
 }
 

lib/src/Uncertainty/Distribution/ClaytonCopulaFactory.cxx

@@ -66,6 +66,7 @@ ClaytonCopula ClaytonCopulaFactory::buildAsClaytonCopula(const Sample & sample)
   if (tau == 1) throw InvalidArgumentException(HERE) << "Error: cannot build a ClaytonCopula distribution from a sample with Kendall tau equal to 1";
   ClaytonCopula result(2.0 * tau / (1.0 - tau));
   result.setDescription(sample.getDescription());
+  adaptToKnownParameter(sample, &result);
   return result;
 }
 

lib/src/Uncertainty/Distribution/DirichletFactory.cxx

@@ -160,6 +160,7 @@ Dirichlet DirichletFactory::buildAsDirichlet(const Sample & sample) const
   }
   Dirichlet result(theta);
   result.setDescription(sample.getDescription());
+  adaptToKnownParameter(sample, &result);
   return result;
 }
 

lib/src/Uncertainty/Distribution/ExponentialFactory.cxx

@@ -19,6 +19,7 @@
  *
  */
 #include "openturns/ExponentialFactory.hxx"
+#include "openturns/MaximumLikelihoodFactory.hxx"
 #include "openturns/SpecFunc.hxx"
 #include "openturns/PersistentObjectFactory.hxx"
 
@@ -73,6 +74,7 @@ Exponential ExponentialFactory::buildAsExponential(const Sample & sample) const
   const Scalar lambda = 1.0 / (mean - gamma);
   Exponential result(lambda, gamma);
   result.setDescription(sample.getDescription());
+  adaptToKnownParameter(sample, &result);
   return result;
 }
 

lib/src/Uncertainty/Distribution/FarlieGumbelMorgensternCopulaFactory.cxx

@@ -76,6 +76,7 @@ FarlieGumbelMorgensternCopula FarlieGumbelMorgensternCopulaFactory::buildAsFarli
   }
   FarlieGumbelMorgensternCopula result(theta);
   result.setDescription(sample.getDescription());
+  adaptToKnownParameter(sample, &result);
   return result;
 }
 

lib/src/Uncertainty/Distribution/FisherSnedecorFactory.cxx

@@ -101,8 +101,10 @@ FisherSnedecor FisherSnedecorFactory::buildMethodOfLikelihoodMaximization(const
   factory.setOptimizationAlgorithm(solver);
 
   // Configure bounds
-  Interval bounds(parametersLowerBound, Point(dimension, SpecFunc::MaxScalar), Interval::BoolCollection(dimension, true), Interval::BoolCollection(dimension, false));
-  factory.setOptimizationBounds(bounds);
+  const Interval bounds(parametersLowerBound, Point(dimension, SpecFunc::MaxScalar), Interval::BoolCollection(dimension, true), Interval::BoolCollection(dimension, false));
+  // Use bounds only for unknown parameters
+  factory.setOptimizationBounds(bounds.getMarginal(knownParameterIndices_.complement(bounds.getDimension())));
+  factory.setKnownParameter(knownParameterValues_, knownParameterIndices_);
 
   return buildAsFisherSnedecor(factory.buildParameter(sample));
 }

lib/src/Uncertainty/Distribution/FrankCopulaFactory.cxx

@@ -102,6 +102,7 @@ FrankCopula FrankCopulaFactory::buildAsFrankCopula(const Sample & sample) const
   theta = solver.solve(f, tau, minTheta, maxTheta, minTau, maxTau);
   FrankCopula result(isTauNegative ? -theta : theta);
   result.setDescription(sample.getDescription());
+  adaptToKnownParameter(sample, &result);
   return result;
 }
 

lib/src/Uncertainty/Distribution/FrechetFactory.cxx

@@ -94,7 +94,10 @@ Frechet FrechetFactory::buildAsFrechet(const Sample & sample) const
   const Scalar margin = std::max(1.0, ResourceMap::GetAsScalar("FrechetFactory-BoundMargin"));
   const Point lower = {betaFrechet / margin, alphaFrechet / margin, gamma - margin * std::abs(gamma)};
   const Point upper = {margin * betaFrechet,  margin * alphaFrechet, gamma + margin * std::abs(gamma)};
-  mleFactory.setOptimizationBounds(Interval(lower, upper));
+  const Interval bounds(lower, upper);
+  // Use bounds only for unknown parameters
+  mleFactory.setOptimizationBounds(bounds.getMarginal(knownParameterIndices_.complement(bounds.getDimension())));
+  mleFactory.setKnownParameter(knownParameterValues_, knownParameterIndices_);
   const Point parameters(mleFactory.buildParameter(sample));
   return buildAsFrechet(parameters);
 }

lib/src/Uncertainty/Distribution/GammaFactory.cxx

@@ -75,6 +75,7 @@ Gamma GammaFactory::buildAsGamma(const Sample & sample) const
   lambda /= sigma;
   Gamma result(k, lambda, gamma);
   result.setDescription(sample.getDescription());
+  adaptToKnownParameter(sample, &result);
   return result;
 }
 

lib/src/Uncertainty/Distribution/GeometricFactory.cxx

@@ -75,6 +75,7 @@ Geometric GeometricFactory::buildAsGeometric(const Sample & sample) const
   }
   Geometric result(size / sum);
   result.setDescription(sample.getDescription());
+  adaptToKnownParameter(sample, &result);
   return result;
 }
 

lib/src/Uncertainty/Distribution/GumbelCopulaFactory.cxx

@@ -66,6 +66,7 @@ GumbelCopula GumbelCopulaFactory::buildAsGumbelCopula(const Sample & sample) con
   if (tau == 1) throw InvalidArgumentException(HERE) << "Error: cannot build a GumbelCopula distribution from a sample with Kendall tau equal to 1";
   GumbelCopula result(1.0 / (1.0 - tau));
   result.setDescription(sample.getDescription());
+  adaptToKnownParameter(sample, &result);
   return result;
 }
 

lib/src/Uncertainty/Distribution/GumbelFactory.cxx

@@ -73,6 +73,7 @@ Gumbel GumbelFactory::buildAsGumbel(const Sample & sample) const
   const Point parameters = {mu, sigma};
   Gumbel result(buildAsGumbel(GumbelMuSigma()(parameters)));
   result.setDescription(sample.getDescription());
+  adaptToKnownParameter(sample, &result);
   return result;
 }
 

lib/src/Uncertainty/Distribution/InverseNormalFactory.cxx

@@ -87,6 +87,7 @@ InverseNormal InverseNormalFactory::buildAsInverseNormal(const Sample & sample)
     lambda = std::pow(mu, 3) / (sigma * sigma);
   InverseNormal result(mu, lambda);
   result.setDescription(sample.getDescription());
+  adaptToKnownParameter(sample, &result);
   return result;
 }
 

lib/src/Uncertainty/Distribution/LaplaceFactory.cxx

@@ -68,6 +68,7 @@ Laplace LaplaceFactory::buildAsLaplace(const Sample & sample) const
   if (tau == 0) throw InvalidArgumentException(HERE) << "Error: cannot build a Laplace distribution with infinite lambda.";
   Laplace result(mu, size / tau);
   result.setDescription(sample.getDescription());
+  adaptToKnownParameter(sample, &result);
   return result;
 }
 

lib/src/Uncertainty/Distribution/LeastSquaresDistributionFactory.cxx

@@ -189,6 +189,9 @@ Point LeastSquaresDistributionFactory::buildParameter(const Sample & sample) con
   if (knownParameterValues_.getSize() != knownParameterIndices_.getSize())
     throw InvalidArgumentException(HERE) << "Error: known values size must match indices";
 
+  // Quick return if all the parameter values are known
+  if (knownParameterValues_.getSize() == effectiveParameterSize) return knownParameterValues_;
+
   LeastSquaresFactoryResidualEvaluation residualEvaluation(sample, distribution_, knownParameterValues_, knownParameterIndices_);
   Function residual(residualEvaluation.clone());
 
@@ -294,46 +297,26 @@ OptimizationAlgorithm LeastSquaresDistributionFactory::getOptimizationAlgorithm(
   return solver_;
 }
 
-void LeastSquaresDistributionFactory::setKnownParameter(const Point & values,
-    const Indices & indices)
-{
-  if (values.getSize() != indices.getSize())
-    throw InvalidArgumentException(HERE) << "Known parameters values and indices must have the same size";
-  if (!indices.check(distribution_.getParameter().getSize()))
-    throw InvalidArgumentException(HERE) << "Know parameters indices must be < parameter dimension";
-  knownParameterValues_ = values;
-  knownParameterIndices_ = indices;
-}
-
-Indices LeastSquaresDistributionFactory::getKnownParameterIndices() const
-{
-  return knownParameterIndices_;
-}
-
-Point LeastSquaresDistributionFactory::getKnownParameterValues() const
-{
-  return knownParameterValues_;
-}
-
-
 /* Method save() stores the object through the StorageManager */
 void LeastSquaresDistributionFactory::save(Advocate & adv) const
 {
   DistributionFactoryImplementation::save(adv);
-  adv.saveAttribute("knownParameterValues_", knownParameterValues_);
-  adv.saveAttribute("knownParameterIndices_", knownParameterIndices_);
   adv.saveAttribute("optimizationBounds_", optimizationBounds_);
   adv.saveAttribute("optimizationInequalityConstraint_", optimizationInequalityConstraint_);
+  adv.saveAttribute("knownParameterValues_", knownParameterValues_);
+  adv.saveAttribute("knownParameterIndices_", knownParameterIndices_);
 }
 
 /* Method load() reloads the object from the StorageManager */
 void LeastSquaresDistributionFactory::load(Advocate & adv)
 {
   DistributionFactoryImplementation::load(adv);
-  adv.loadAttribute("knownParameterValues_", knownParameterValues_);
-  adv.loadAttribute("knownParameterIndices_", knownParameterIndices_);
   adv.loadAttribute("optimizationBounds_", optimizationBounds_);
   adv.loadAttribute("optimizationInequalityConstraint_", optimizationInequalityConstraint_);
+  if (adv.hasAttribute("knownParameterValues_"))
+    adv.loadAttribute("knownParameterValues_", knownParameterValues_);
+  if (adv.hasAttribute("knownParameterIndices_"))
+    adv.loadAttribute("knownParameterIndices_", knownParameterIndices_);
 }
 
 END_NAMESPACE_OPENTURNS

lib/src/Uncertainty/Distribution/LogNormalFactory.cxx

@@ -231,41 +231,48 @@ LogNormal LogNormalFactory::buildAsLogNormal(const Sample & sample,
   const UnsignedInteger size = sample.getSize();
   if (size < 3) throw InvalidArgumentException(HERE) << "Error: cannot build a LogNormal distribution from a sample of size < 3";
   if (sample.getDimension() != 1) throw InvalidArgumentException(HERE) << "Error: can build a LogNormal distribution only from a sample of dimension 1, here dimension=" << sample.getDimension();
+  LogNormal result;
   switch (method)
   {
     case 0:
       try
       {
-        return buildMethodOfLocalLikelihoodMaximization(sample);
+        result = buildMethodOfLocalLikelihoodMaximization(sample);
+        break;
       }
       catch (const InvalidArgumentException & ex)
       {
         // We switch to the moment estimate
         LOGWARN(OSS() << ex.what());
-        return buildAsLogNormal(sample, 1);
+        result = buildAsLogNormal(sample, 1);
+        break;
       }
       break;
     case 1:
       try
       {
-        return buildMethodOfModifiedMoments(sample);
+        result = buildMethodOfModifiedMoments(sample);
+        break;
       }
       catch (const InvalidArgumentException & ex)
       {
         // We switch to the moment estimate
         LOGWARN(OSS() << ex.what());
-        return buildAsLogNormal(sample, 2);
+        result = buildAsLogNormal(sample, 2);
+        break;
       }
       break;
     case 2:
-      return buildMethodOfMoments(sample);
+      result = buildMethodOfMoments(sample);
       break;
     case 3:
-      return buildMethodOfLeastSquares(sample);
+      result = buildMethodOfLeastSquares(sample);
       break;
     default:
       throw InvalidArgumentException(HERE) << "Error: invalid value=" << method << " for the key 'LogNormalFactory-EstimationMethod' in ResourceMap";
   }
+  adaptToKnownParameter(sample, &result);
+  return result;  
 }
 
 LogNormal LogNormalFactory::buildAsLogNormal(const Point & parameters) const

lib/src/Uncertainty/Distribution/LogUniformFactory.cxx

@@ -76,6 +76,7 @@ LogUniform LogUniformFactory::buildAsLogUniform(const Sample & sample) const
   Scalar bLog = std::log(b);
   LogUniform result(aLog, bLog);
   result.setDescription(sample.getDescription());
+  adaptToKnownParameter(sample, &result);
   return result;
 }
 

lib/src/Uncertainty/Distribution/LogisticFactory.cxx

@@ -73,6 +73,7 @@ Logistic LogisticFactory::buildAsLogistic(const Sample & sample) const
   if (!(beta > 0.0)) throw InvalidArgumentException(HERE) << "Error: can build a Logistic distribution only if beta > 0.0, here beta=" << beta;
   Logistic result(mu, beta);
   result.setDescription(sample.getDescription());
+  adaptToKnownParameter(sample, &result);
   return result;
 }
 

lib/src/Uncertainty/Distribution/MaximumLikelihoodFactory.cxx

@@ -287,6 +287,9 @@ Point MaximumLikelihoodFactory::buildParameter(const Sample & sample) const
   if (knownParameterValues_.getSize() != knownParameterIndices_.getSize())
     throw InvalidArgumentException(HERE) << "Error: known values size must match indices";
 
+  // Quick return if all the parameter values are known
+  if (knownParameterValues_.getSize() == effectiveParameterSize) return knownParameterValues_;
+
   // Define evaluation
   LogLikelihoodEvaluation logLikelihoodWrapper(sample, distribution_, knownParameterValues_, knownParameterIndices_);
   Function logLikelihood(logLikelihoodWrapper.clone());
@@ -398,31 +401,10 @@ OptimizationAlgorithm MaximumLikelihoodFactory::getOptimizationAlgorithm() const
   return solver_;
 }
 
-void MaximumLikelihoodFactory::setKnownParameter(const Point & values,
-    const Indices & indices)
-{
-  if (values.getSize() != indices.getSize())
-    throw InvalidArgumentException(HERE) << "Known parameters values and indices must have the same size";
-  knownParameterValues_ = values;
-  knownParameterIndices_ = indices;
-}
-
-Indices MaximumLikelihoodFactory::getKnownParameterIndices() const
-{
-  return knownParameterIndices_;
-}
-
-Point MaximumLikelihoodFactory::getKnownParameterValues() const
-{
-  return knownParameterValues_;
-}
-
 /* Method save() stores the object through the StorageManager */
 void MaximumLikelihoodFactory::save(Advocate & adv) const
 {
   DistributionFactoryImplementation::save(adv);
-  adv.saveAttribute("knownParameterValues_", knownParameterValues_);
-  adv.saveAttribute("knownParameterIndices_", knownParameterIndices_);
   adv.saveAttribute("optimizationBounds_", optimizationBounds_);
   adv.saveAttribute("optimizationInequalityConstraint_", optimizationInequalityConstraint_);
 }
@@ -431,10 +413,12 @@ void MaximumLikelihoodFactory::save(Advocate & adv) const
 void MaximumLikelihoodFactory::load(Advocate & adv)
 {
   DistributionFactoryImplementation::load(adv);
-  adv.loadAttribute("knownParameterValues_", knownParameterValues_);
-  adv.loadAttribute("knownParameterIndices_", knownParameterIndices_);
   adv.loadAttribute("optimizationBounds_", optimizationBounds_);
   adv.loadAttribute("optimizationInequalityConstraint_", optimizationInequalityConstraint_);
+  if (adv.hasAttribute("knownParameterValues_"))
+    adv.loadAttribute("knownParameterValues_", knownParameterValues_);
+  if (adv.hasAttribute("knownParameterIndices_"))
+    adv.loadAttribute("knownParameterIndices_", knownParameterIndices_);
 }
 
 

lib/src/Uncertainty/Distribution/MeixnerDistributionFactory.cxx

@@ -76,6 +76,7 @@ MeixnerDistribution MeixnerDistributionFactory::buildAsMeixnerDistribution(const
   const Scalar mu = m - alpha * delta * std::tan(0.5 * beta);
   MeixnerDistribution result(alpha, beta, delta, mu);
   result.setDescription(sample.getDescription());
+  adaptToKnownParameter(sample, &result);
   return result;
 }