Fix Wishart and InverseWishart covariance matrix

Fixes openturns#1241

lib/src/Uncertainty/Distribution/InverseWishart.cxx

@@ -238,22 +238,33 @@ void InverseWishart::computeCovariance() const
 {
   const UnsignedInteger p = cholesky_.getDimension();
   const Scalar den = (nu_ - p) * std::pow(nu_ - p - 1.0, 2) * (nu_ - p - 3.0);
-  if (!(den > 0.0)) throw NotDefinedException(HERE) << "Error: the covariance of the inverse Wishart distribution is defined only if nu > p+3";
-  const CovarianceMatrix V(getV());
-  covariance_ = CovarianceMatrix(getDimension());
-  UnsignedInteger indexRow = 0;
+  if (!(den > 0.0)) throw NotDefinedException(HERE) << "Error: the covariance of the inverse Wishart distribution is defined only if nu > p+3, here nu = " << nu_ << " and p = " << p;
+
+  // Get the collection of Indices of the random matrix
+  // in the order of the corresponding flattened random vector
+  Collection<Indices>  matrixIndices;
   for (UnsignedInteger i = 0; i < p; ++i)
+  {
     for (UnsignedInteger j = 0; j <= i; ++j)
     {
-      UnsignedInteger indexColumn = 0;
-      for (UnsignedInteger m = 0; m <= i; ++m)
-        for (UnsignedInteger n = 0; n <= j; ++n)
-        {
-          covariance_(indexRow, indexColumn) = (2.0 * V(i, j) * V(m, n) + (nu_ - p - 1.0) * (V(i, m) * V(j, n) + V(i, n) * V(m, j))) / den;
-          ++indexColumn;
-        }
-      ++indexRow;
+      matrixIndices.add(Indices({i, j}));
+    }
+  }
+
+  // Populate the covariance matrix of the flattened random vector
+  const CovarianceMatrix V(getV());
+  covariance_ = CovarianceMatrix(getDimension());
+  for (UnsignedInteger row = 0; row < matrixIndices.getSize(); ++row)
+  {
+    for (UnsignedInteger col = 0; col <= row; ++col)
+    {
+      const UnsignedInteger irow = matrixIndices[row][0];
+      const UnsignedInteger icol = matrixIndices[col][0];
+      const UnsignedInteger jrow = matrixIndices[row][1];
+      const UnsignedInteger jcol = matrixIndices[col][1];
+      covariance_(row, col) = ((nu_ - p + 1.0) * V(irow, jcol) * V(icol, jrow) + (nu_ - p - 1.0) * V(irow, icol) * V(jrow, jcol)) / den;
     }
+  }
   isAlreadyComputedCovariance_ = true;
 }
 

lib/src/Uncertainty/Distribution/Wishart.cxx

@@ -229,6 +229,38 @@ void Wishart::computeMean() const
   isAlreadyComputedMean_ = true;
 }
 
+/* Compute the covariance of the distribution */
+void Wishart::computeCovariance() const
+{
+  // Get the collection of Indices of the random matrix
+  // in the order of the corresponding flattened random vector
+  const UnsignedInteger p = cholesky_.getDimension();
+  Collection<Indices>  matrixIndices;
+  for (UnsignedInteger i = 0; i < p; ++i)
+  {
+    for (UnsignedInteger j = 0; j <= i; ++j)
+    {
+      matrixIndices.add(Indices({i, j}));
+    }
+  }
+
+  // Populate the covariance matrix of the flattened random vector
+  const CovarianceMatrix V(getV());
+  covariance_ = CovarianceMatrix(getDimension());
+  for (UnsignedInteger row = 0; row < matrixIndices.getSize(); ++row)
+  {
+    for (UnsignedInteger col = 0; col <= row; ++col)
+    {
+      const UnsignedInteger irow = matrixIndices[row][0];
+      const UnsignedInteger icol = matrixIndices[col][0];
+      const UnsignedInteger jrow = matrixIndices[row][1];
+      const UnsignedInteger jcol = matrixIndices[col][1];
+      covariance_(row, col) = nu_ * (V(irow, jcol) *  V(icol, jrow) + V(irow, icol) * V(jrow, jcol));
+    }
+  }
+  isAlreadyComputedCovariance_ = true;
+}
+
 /* Compute the entropy of the distribution */
 Scalar Wishart::computeEntropy() const
 {

lib/src/Uncertainty/Distribution/openturns/Wishart.hxx

@@ -124,6 +124,9 @@ private:
   /** Compute the mean of the distribution */
   void computeMean() const override;
 
+  /** Compute the covariance of the distribution */
+  void computeCovariance() const override;
+
   /** Compute the numerical range of the distribution given the parameters values */
   void computeRange() override;
 

lib/test/t_InverseWishart_std.cxx

@@ -144,6 +144,16 @@ int main(int, char *[])
     Arcsine::PointWithDescriptionCollection parameters = distribution.getParametersCollection();
     fullprint << "parameters=" << parameters << std::endl;
     fullprint << "Standard representative=" << distribution.getStandardRepresentative().__str__() << std::endl;
+
+    // Verify covariance
+    CovarianceMatrix matrix(2);
+    matrix(0, 0) = 100.0;
+    matrix(1, 0) = 50.0;
+    matrix(1, 1) = 80.0;
+    // Empirical estimators are typically bad for this distribution, so we use a high degree of freedom.
+    InverseWishart multidimensional(matrix, 40.0);
+    CovarianceMatrix empirical(multidimensional.getSample(100000).computeCovariance());
+    assert_almost_equal(multidimensional.getCovariance(), empirical, 0.1, 0.0);
   }
   catch (TestFailed & ex)
   {

lib/test/t_Wishart_std.cxx

@@ -143,6 +143,15 @@ int main(int, char *[])
     Arcsine::PointWithDescriptionCollection parameters = distribution.getParametersCollection();
     fullprint << "parameters=" << parameters << std::endl;
     fullprint << "Standard representative=" << distribution.getStandardRepresentative().__str__() << std::endl;
+
+    // Verify covariance
+    CovarianceMatrix matrix(2);
+    matrix(0, 0) = 10.0;
+    matrix(1, 0) = 5.0;
+    matrix(1, 1) = 8.0;
+    Wishart multidimensional(matrix, 3.0);
+    CovarianceMatrix empirical(multidimensional.getSample(100000).computeCovariance());
+    assert_almost_equal(multidimensional.getCovariance(), empirical, 0.02, 0.0);
   }
   catch (TestFailed & ex)
   {