Improved DistributionImplementation

The computation of Kendall tau has been simplified and fixed. It fixes openturns#2631

lib/src/Uncertainty/Model/DistributionImplementation.cxx

@@ -3242,50 +3242,53 @@ CorrelationMatrix DistributionImplementation::getSpearmanCorrelation() const
   return getCopula().getSpearmanCorrelation();
 }
 
+/* This helper class is here to compute the Kendall tau of bivariate distributions
+
+   If the distribution is a copula with CDF C and PDF c, then:
+
+   tau = 4\int_0^1\int_0^1 C(u,v)c(u,v)dudv - 1
+
+   If the distribution is general, with CDF F, PDF f, marginal quantiles Qx and Qy, marginal PDF fx and fy, then:
+
+   C(u,v)=F(Qx(u),Qy(v))
+   c(u,v)=d²C(u,v)/dudv=Qx'(u)Qy'(v)f(Qx(u),Qy(v))
+
+   and:
+
+   tau = 4\int_0^1\int_0^1 F(Qx(u),Qy(v))f(Qx(u),Qy(v))Qx'(u)Qy'(v)dudv - 1
+
+   then, with x=Qx(u) and y=Qy(u), dxdy=Qx'(u)duQy'(v)dv and
+
+       = 4\int_R\int_R F(x,y)f(x,y)dxdy - 1
+
+   So in all the cases, we have to integrate the product CDFxPDF over the range
+   of the distribution
+ */
 struct DistributionImplementationKendallTauWrapper
 {
   DistributionImplementationKendallTauWrapper(const Distribution & distribution)
     : distribution_(distribution)
   {
-    if (!distribution.isCopula())
-    {
-      const UnsignedInteger dimension = distribution.getDimension();
-      marginalCollection_ = Collection<Distribution>(dimension);
-      for (UnsignedInteger i = 0; i < dimension; ++i)
-        marginalCollection_[i] = distribution.getMarginal(i);
-    }
-  }
-
-  Point kernelForCopula(const Point & point) const
-  {
-    return Point(1, distribution_.computeCDF(point) * distribution_.computePDF(point));
+    // Nothing to do
   }
 
-  Point kernelForDistribution(const Point & point) const
+  Point kernel(const Point & point) const
   {
-    const UnsignedInteger dimension = distribution_.getDimension();
-    Point x(dimension);
-    Scalar factor = 1.0;
-    for (UnsignedInteger i = 0; i < dimension; ++i)
-    {
-      const Point xi(marginalCollection_[i].computeQuantile(point[i]));
-      x[i] = xi[0];
-      factor *= marginalCollection_[i].computePDF(xi);
-      if (std::abs(factor) < SpecFunc::Precision) return Point(1, 0.0);
-    }
-    return Point(1, distribution_.computeCDF(point) * distribution_.computePDF(x) / factor);
+    const Scalar pdf = distribution_.computePDF(point);
+    if (std::abs(pdf) < SpecFunc::Precision) return Point(1, 0.0);
+    return Point(1, distribution_.computeCDF(point) * pdf);
   }
 
   const Distribution & distribution_;
-  Collection<Distribution> marginalCollection_;
-}; // DistributionImplementationKendallTauWrapperx
+}; // DistributionImplementationKendallTauWrapper
 
 /* Get the Kendall concordance of the distribution */
 CorrelationMatrix DistributionImplementation::getKendallTau() const
 {
   CorrelationMatrix tau(dimension_);
   // First special case: independent marginals
-  if (hasIndependentCopula()) return tau;
+  if (hasIndependentCopula())
+    return tau;
   // Second special case: elliptical distribution
   if (hasEllipticalCopula())
   {
@@ -3297,7 +3300,6 @@ CorrelationMatrix DistributionImplementation::getKendallTau() const
   }
   // General case
   const IteratedQuadrature integrator;
-  const Interval square(2);
   // Performs the integration in the strictly lower triangle of the tau matrix
   Indices indices(2);
   for(UnsignedInteger rowIndex = 0; rowIndex < dimension_; ++rowIndex)
@@ -3306,18 +3308,14 @@ CorrelationMatrix DistributionImplementation::getKendallTau() const
     for (UnsignedInteger columnIndex = rowIndex + 1; columnIndex < dimension_; ++columnIndex)
     {
       indices[1] = columnIndex;
-      const Distribution marginalDistribution(getMarginal(indices).getImplementation());
+      const Distribution marginalDistribution(getMarginal(indices));
       if (!marginalDistribution.hasIndependentCopula())
       {
         // Build the integrand
         const DistributionImplementationKendallTauWrapper functionWrapper(marginalDistribution);
-        Function function;
-        if (isCopula())
-          function = (bindMethod<DistributionImplementationKendallTauWrapper, Point, Point>(functionWrapper, &DistributionImplementationKendallTauWrapper::kernelForCopula, 2, 1));
-        else
-          function = (bindMethod<DistributionImplementationKendallTauWrapper, Point, Point>(functionWrapper, &DistributionImplementationKendallTauWrapper::kernelForDistribution, 2, 1));
-        tau(rowIndex, columnIndex) = integrator.integrate(function, square)[0];
-      }
+        const Function function(bindMethod<DistributionImplementationKendallTauWrapper, Point, Point>(functionWrapper, &DistributionImplementationKendallTauWrapper::kernel, 2, 1));
+        tau(rowIndex, columnIndex) = 4.0 * integrator.integrate(function, marginalDistribution.getRange())[0] - 1.0;
+      } // !independent margins
     } // loop over column indices
   } // loop over row indices
   return tau;

lib/test/t_BlockIndependentDistribution_std.expout

@@ -41,7 +41,7 @@ entropy (MC)=10.9656
 covariance=class=CovarianceMatrix dimension=7 implementation=class=MatrixImplementation name=Unnamed rows=7 columns=7 values=[1,0.184,0,0,0,0,0,0.184,1,0,0,0,0,0,0,0,4,2,1,0,0,0,0,2,4,0,0,0,0,0,1,0,4,0,0,0,0,0,0,0,1,0.06227,0,0,0,0,0,0.06227,1]
 correlation=class=CorrelationMatrix dimension=7 implementation=class=MatrixImplementation name=Unnamed rows=7 columns=7 values=[1,0.184,0,0,0,0,0,0.184,1,0,0,0,0,0,0,0,1,0.5,0.25,0,0,0,0,0.5,1,0,0,0,0,0,0.25,0,1,0,0,0,0,0,0,0,1,0.06227,0,0,0,0,0,0.06227,1]
 spearman=class=CorrelationMatrix dimension=7 implementation=class=MatrixImplementation name=Unnamed rows=7 columns=7 values=[1,0.1924,0,0,0,0,0,0.1924,1,0,0,0,0,0,0,0,1,0.4826,0.2394,0,0,0,0,0.4826,1,0,0,0,0,0,0.2394,0,1,0,0,0,0,0,0,0,1,0.08306,0,0,0,0,0,0.08306,1]
-kendall=class=CorrelationMatrix dimension=7 implementation=class=MatrixImplementation name=Unnamed rows=7 columns=7 values=[1,0.493,0,0,0,0,0,0.493,1,0,0,0,0,0,0,0,1,0.3333,0.1609,0,0,0,0,0.3333,1,0,0,0,0,0,0.1609,0,1,0,0,0,0,0,0,0,1,0.148,0,0,0,0,0,0.148,1]
+kendall=class=CorrelationMatrix dimension=7 implementation=class=MatrixImplementation name=Unnamed rows=7 columns=7 values=[1,0.1288,0,0,0,0,0,0.1288,1,0,0,0,0,0,0,0,1,0.3333,0.1609,0,0,0,0,0.3333,1,0,0,0,0,0,0.1609,0,1,0,0,0,0,0,0,0,1,0.05542,0,0,0,0,0,0.05542,1]
 conditional PDF=0.556752
 conditional CDF=0.490621
 conditional quantile=0.823038

lib/test/t_EmpiricalBernsteinCopula_std.expout

@@ -78,7 +78,7 @@ beta=0.974525
 covariance=class=CovarianceMatrix dimension=2 implementation=class=MatrixImplementation name=Unnamed rows=2 columns=2 values=[0.0833,0.0104,0.0104,0.0833]
 correlation=class=CorrelationMatrix dimension=2 implementation=class=MatrixImplementation name=Unnamed rows=2 columns=2 values=[1,0.125,0.125,1]
 spearman=class=CorrelationMatrix dimension=2 implementation=class=MatrixImplementation name=Unnamed rows=2 columns=2 values=[1,0.125,0.125,1]
-kendall=class=CorrelationMatrix dimension=2 implementation=class=MatrixImplementation name=Unnamed rows=2 columns=2 values=[1,0.271,0.271,1]
+kendall=class=CorrelationMatrix dimension=2 implementation=class=MatrixImplementation name=Unnamed rows=2 columns=2 values=[1,0.085,0.085,1]
 margin=class=Uniform name=Uniform dimension=1 a=0 b=1
 margin PDF=1
 margin CDF=0.25

lib/test/t_ExtremeValueCopula_std.expout

@@ -18,7 +18,7 @@ CDF(quantile)=0.5
 covariance=class=CovarianceMatrix dimension=2 implementation=class=MatrixImplementation name=Unnamed rows=2 columns=2 values=[0.0833333,0.0489531,0.0489531,0.0833333]
 correlation=class=CorrelationMatrix dimension=2 implementation=class=MatrixImplementation name=Unnamed rows=2 columns=2 values=[1,0.587437,0.587437,1]
 spearman=class=CorrelationMatrix dimension=2 implementation=class=MatrixImplementation name=Unnamed rows=2 columns=2 values=[1,0.587437,0.587437,1]
-kendall=class=CorrelationMatrix dimension=2 implementation=class=MatrixImplementation name=Unnamed rows=2 columns=2 values=[1,0.3546,0.3546,1]
+kendall=class=CorrelationMatrix dimension=2 implementation=class=MatrixImplementation name=Unnamed rows=2 columns=2 values=[1,0.418399,0.418399,1]
 CDF(x|y)=0.660204
 Quantile(p|y)=0.6
 margin=class=IndependentCopula name=IndependentCopula dimension=1

lib/test/t_GalambosCopula_std.expout

@@ -18,7 +18,7 @@ CDF(quantile)=0.5
 covariance=class=CovarianceMatrix dimension=2 implementation=class=MatrixImplementation name=Unnamed rows=2 columns=2 values=[0.0833333,0.0241064,0.0241064,0.0833333]
 correlation=class=CorrelationMatrix dimension=2 implementation=class=MatrixImplementation name=Unnamed rows=2 columns=2 values=[1,0.289277,0.289277,1]
 spearman=class=CorrelationMatrix dimension=2 implementation=class=MatrixImplementation name=Unnamed rows=2 columns=2 values=[1,0.289277,0.289277,1]
-kendall=class=CorrelationMatrix dimension=2 implementation=class=MatrixImplementation name=Unnamed rows=2 columns=2 values=[1,0.299108,0.299108,1]
+kendall=class=CorrelationMatrix dimension=2 implementation=class=MatrixImplementation name=Unnamed rows=2 columns=2 values=[1,0.19643,0.19643,1]
 CDF(x|y)=0.627074
 Quantile(p|y)=0.6
 margin=class=IndependentCopula name=IndependentCopula dimension=1

lib/test/t_JoeCopula_std.expout

@@ -18,7 +18,7 @@ CDF(quantile)=0.5
 covariance=class=CovarianceMatrix dimension=2 implementation=class=MatrixImplementation name=Unnamed rows=2 columns=2 values=[0.0833333,0.0302653,0.0302653,0.0833333]
 correlation=class=CorrelationMatrix dimension=2 implementation=class=MatrixImplementation name=Unnamed rows=2 columns=2 values=[1,0.363184,0.363184,1]
 spearman=class=CorrelationMatrix dimension=2 implementation=class=MatrixImplementation name=Unnamed rows=2 columns=2 values=[1,0.363184,0.363184,1]
-kendall=class=CorrelationMatrix dimension=2 implementation=class=MatrixImplementation name=Unnamed rows=2 columns=2 values=[1,0.314282,0.314282,1]
+kendall=class=CorrelationMatrix dimension=2 implementation=class=MatrixImplementation name=Unnamed rows=2 columns=2 values=[1,0.257128,0.257128,1]
 CDF(x|y)=0.660023
 Quantile(p|y)=0.6
 margin=class=IndependentCopula name=IndependentCopula dimension=1

lib/test/t_PlackettCopula_std.expout

@@ -30,7 +30,7 @@ beta=0.974228
 covariance=class=CovarianceMatrix dimension=2 implementation=class=MatrixImplementation name=Unnamed rows=2 columns=2 values=[0.0833333,0.024761,0.024761,0.0833333]
 correlation=class=CorrelationMatrix dimension=2 implementation=class=MatrixImplementation name=Unnamed rows=2 columns=2 values=[1,0.297132,0.297132,1]
 spearman=class=CorrelationMatrix dimension=2 implementation=class=MatrixImplementation name=Unnamed rows=2 columns=2 values=[1,0.297132,0.297132,1]
-kendall=class=CorrelationMatrix dimension=2 implementation=class=MatrixImplementation name=Unnamed rows=2 columns=2 values=[1,0.300338,0.300338,1]
+kendall=class=CorrelationMatrix dimension=2 implementation=class=MatrixImplementation name=Unnamed rows=2 columns=2 values=[1,0.201351,0.201351,1]
 margin=class=IndependentCopula name=IndependentCopula dimension=1
 margin PDF=1
 margin CDF=0.25

python/test/t_BlockIndependentDistribution_std.expout

@@ -70,13 +70,13 @@ spearman= 7x7
  [ 0       0       0       0       0       1       0.08306 ]
  [ 0       0       0       0       0       0.08306 1       ]]
 kendall= 7x7
-[[ 1      0.493  0      0      0      0      0      ]
- [ 0.493  1      0      0      0      0      0      ]
- [ 0      0      1      0.3333 0.1609 0      0      ]
- [ 0      0      0.3333 1      0      0      0      ]
- [ 0      0      0.1609 0      1      0      0      ]
- [ 0      0      0      0      0      1      0.148  ]
- [ 0      0      0      0      0      0.148  1      ]]
+[[ 1       0.1288  0       0       0       0       0       ]
+ [ 0.1288  1       0       0       0       0       0       ]
+ [ 0       0       1       0.3333  0.1609  0       0       ]
+ [ 0       0       0.3333  1       0       0       0       ]
+ [ 0       0       0.1609  0       1       0       0       ]
+ [ 0       0       0       0       0       1       0.05542 ]
+ [ 0       0       0       0       0       0.05542 1       ]]
 conditional PDF=0.55675
 conditional CDF=0.49062
 conditional quantile=0.82304

python/test/t_GalambosCopula_std.py

@@ -47,7 +47,7 @@
 )
 ott.assert_almost_equal(
     copula.getKendallTau(),
-    ot.CovarianceMatrix(2, [1, 0.299108, 0.299108, 1]),
+    ot.CovarianceMatrix(2, [1, 0.19643, 0.19643, 1]),
     1e-5,
     0.0,
 )