Fixed Dirichlet

The formula for skewness was wrong. Now, the 1D marginal distributions are explicit Beta distributions. It fixes openturns#2631
regislebrun · Apr 28, 2024 · ac5f46b · ac5f46b
1 parent 8eb4958
commit ac5f46b
Show file tree

Hide file tree

Showing 3 changed files with 11 additions and 7 deletions.
diff --git a/lib/src/Uncertainty/Distribution/Dirichlet.cxx b/lib/src/Uncertainty/Distribution/Dirichlet.cxx
@@ -21,6 +21,7 @@
 #include <cmath>
 #include "openturns/Indices.hxx"
 #include "openturns/Dirichlet.hxx"
+#include "openturns/Beta.hxx"
 #include "openturns/RandomGenerator.hxx"
 #include "openturns/SpecFunc.hxx"
 #include "openturns/DistFunc.hxx"
@@ -469,7 +470,7 @@ Point Dirichlet::getSkewness() const
   for (UnsignedInteger i = 0; i < dimension; ++i)
   {
     const Scalar thetaI = theta_[i];
-    skewness[i] = 2.0 * (sumTheta_ - 2.0 * thetaI) / (sumTheta_ + 2.0) * std::sqrt(sumTheta_ + 1.0) / (thetaI * (sumTheta_ - thetaI));
+    skewness[i] = 2.0 * (sumTheta_ - 2.0 * thetaI) / (sumTheta_ + 2.0) * std::sqrt((sumTheta_ + 1.0) / (thetaI * (sumTheta_ - thetaI)));
   }
   return skewness;
 }
@@ -536,11 +537,14 @@ Distribution Dirichlet::getMarginal(const UnsignedInteger i) const
   const UnsignedInteger dimension = getDimension();
   if (i >= dimension) throw InvalidArgumentException(HERE) << "The index of a marginal distribution must be in the range [0, dim-1]";
   if (dimension == 1) return clone();
+  Beta marginal(theta_[i], sumTheta_ - theta_[i], 0.0, 1.0);
+  /*
   Point thetaMarginal(2);
   thetaMarginal[0] = theta_[i];
   thetaMarginal[1] = sumTheta_ - theta_[i];
   Dirichlet::Implementation marginal(new Dirichlet(thetaMarginal));
-  marginal->setDescription(Description(1, getDescription()[i]));
+  */
+  marginal.setDescription({getDescription()[i]});
   return marginal;
 }
 

diff --git a/lib/test/t_Dirichlet_std.expout b/lib/test/t_Dirichlet_std.expout
@@ -18,7 +18,7 @@ entropy=-0.041475
 entropy (MC)=-0.041557
 mean=class=Point name=Unnamed dimension=1 values=[0.454545]
 standard deviation=class=Point name=Unnamed dimension=1 values=[0.352089]
-skewness=class=Point name=Unnamed dimension=1 values=[0.108715]
+skewness=class=Point name=Unnamed dimension=1 values=[0.148865]
 kurtosis=class=Point name=Unnamed dimension=1 values=[1.98398]
 covariance=class=CovarianceMatrix dimension=1 implementation=class=MatrixImplementation name=Unnamed rows=1 columns=1 values=[0.0661157]
 parameters collection=[[theta : 1.25, sum theta : 1.5]]
@@ -53,7 +53,7 @@ entropy=-0.82087
 entropy (MC)=-0.82062
 mean=class=Point name=Unnamed dimension=2 values=[0.277778,0.333333]
 standard deviation=class=Point name=Unnamed dimension=2 values=[0.384946,0.426401]
-skewness=class=Point name=Unnamed dimension=2 values=[0.35525,0.240534]
+skewness=class=Point name=Unnamed dimension=2 values=[0.71603,0.51025]
 kurtosis=class=Point name=Unnamed dimension=2 values=[2.86651,2.53846]
 covariance=class=CovarianceMatrix dimension=2 implementation=class=MatrixImplementation name=Unnamed rows=2 columns=2 values=[0.0364759,-0.016835,-0.016835,0.040404]
 parameters collection=[[theta : 1.25, sum theta : 3.25],[theta : 1.5, sum theta : 3]]
@@ -65,12 +65,12 @@ conditional quantile=0.42284
 sequential conditional PDF=class=Point name=Unnamed dimension=2 values=[2.12095,1.32086]
 sequential conditional CDF(class=Point name=Unnamed dimension=2 values=[0.05,0.15])=class=Point name=Unnamed dimension=2 values=[0.0893938,0.141993]
 sequential conditional quantile(class=Point name=Unnamed dimension=2 values=[0.0893938,0.141993])=class=Point name=Unnamed dimension=2 values=[0.05,0.15]
-margin=class=Dirichlet name=Dirichlet dimension=1 theta=class=Point name=Unnamed dimension=2 values=[1.25,3.25]
+margin=class=Beta name=Beta dimension=1 alpha=1.25 beta=3.25 a=0 b=1
 margin PDF=0.88989
 margin CDF=0.85595
 margin quantile=class=Point name=Unnamed dimension=1 values=[0.643596]
 margin realization=class=Point name=Unnamed dimension=1 values=[0.337255]
-margin=class=Dirichlet name=Dirichlet dimension=1 theta=class=Point name=Unnamed dimension=2 values=[1.5,3]
+margin=class=Beta name=Beta dimension=1 alpha=1.5 beta=3 a=0 b=1
 margin PDF=1.1601
 margin CDF=0.78445
 margin quantile=class=Point name=Unnamed dimension=1 values=[0.704013]

diff --git a/python/test/t_Dirichlet_std.expout b/python/test/t_Dirichlet_std.expout
@@ -30,7 +30,7 @@ Unilateral confidence interval (upper tail)= [1.8464e-05, 1]
 beta= [0.95]
 mean= class=Point name=Unnamed dimension=1 values=[0.333333]
 standard deviation= class=Point name=Unnamed dimension=1 values=[0.125988]
-skewness= class=Point name=Unnamed dimension=1 values=[1.92418]
+skewness= class=Point name=Unnamed dimension=1 values=[0.680301]
 kurtosis= class=Point name=Unnamed dimension=1 values=[1.90909]
 covariance= class=CovarianceMatrix dimension=1 implementation=class=MatrixImplementation name=Unnamed rows=1 columns=1 values=[0.126984]
 conditional PDF=0.442267