diff --git a/lib/src/Uncertainty/Model/DistributionImplementation.cxx b/lib/src/Uncertainty/Model/DistributionImplementation.cxx index 733336312f2..d801948b688 100644 --- a/lib/src/Uncertainty/Model/DistributionImplementation.cxx +++ b/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(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 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(functionWrapper, &DistributionImplementationKendallTauWrapper::kernelForCopula, 2, 1)); - else - function = (bindMethod(functionWrapper, &DistributionImplementationKendallTauWrapper::kernelForDistribution, 2, 1)); - tau(rowIndex, columnIndex) = integrator.integrate(function, square)[0]; - } + const Function function(bindMethod(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; diff --git a/lib/test/t_BlockIndependentDistribution_std.expout b/lib/test/t_BlockIndependentDistribution_std.expout index 52a232089e6..588770c00e5 100644 --- a/lib/test/t_BlockIndependentDistribution_std.expout +++ b/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 diff --git a/lib/test/t_EmpiricalBernsteinCopula_std.expout b/lib/test/t_EmpiricalBernsteinCopula_std.expout index 9a72cbd43cd..852b5bf2e2c 100644 --- a/lib/test/t_EmpiricalBernsteinCopula_std.expout +++ b/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 diff --git a/lib/test/t_ExtremeValueCopula_std.expout b/lib/test/t_ExtremeValueCopula_std.expout index d0c0b5d7855..c07e14a2641 100644 --- a/lib/test/t_ExtremeValueCopula_std.expout +++ b/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 diff --git a/lib/test/t_GalambosCopula_std.expout b/lib/test/t_GalambosCopula_std.expout index a784f71e281..c7cdf2a585e 100644 --- a/lib/test/t_GalambosCopula_std.expout +++ b/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 diff --git a/lib/test/t_JoeCopula_std.expout b/lib/test/t_JoeCopula_std.expout index ecdc0c616c3..4cb18f32b24 100644 --- a/lib/test/t_JoeCopula_std.expout +++ b/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 diff --git a/lib/test/t_PlackettCopula_std.expout b/lib/test/t_PlackettCopula_std.expout index 81c125ac39c..0e56d242ed5 100644 --- a/lib/test/t_PlackettCopula_std.expout +++ b/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 diff --git a/python/test/t_BlockIndependentDistribution_std.expout b/python/test/t_BlockIndependentDistribution_std.expout index 8d04e37a89a..f6748aa71fb 100644 --- a/python/test/t_BlockIndependentDistribution_std.expout +++ b/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 diff --git a/python/test/t_GalambosCopula_std.py b/python/test/t_GalambosCopula_std.py index bfc178eee88..13486dac39d 100755 --- a/python/test/t_GalambosCopula_std.py +++ b/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, )