[hist] Improve stability of ComputeVariance()

Author: hahnjo

hist/histv7/inc/ROOT/RHistStats.hxx

@@ -132,7 +132,7 @@ public:
    /// \f]
    /// With some rewriting, this is equivalent to:
    /// \f[
-   /// \sigma^2 = \frac{\sum w_i \cdot x_i^2}{\sum w_i} - \frac{(\sum w_i \cdot x_i)^2}{(\sum w_i)^2}
+   /// \sigma^2 = \frac{\sum w_i \cdot x_i^2}{\sum w_i} - \mu^2
    /// \f]
    ///
    /// This function does not include Bessel's correction needed for an unbiased estimator of population variance.
@@ -147,7 +147,8 @@ public:
       if (fSumW == 0) {
          return 0;
       }
-      return (stats.fSumWX2 - stats.fSumWX * stats.fSumWX / fSumW) / fSumW;
+      double mean = ComputeMean(dim);
+      return stats.fSumWX2 / fSumW - mean * mean;
    }
 
    /// Compute the standard deviation of unbinned values.

hist/histv7/test/hist_stats.cxx

@@ -342,6 +342,16 @@ TEST(RHistStats, ComputeKurtosisWeighted)
    EXPECT_FLOAT_EQ(stats.ComputeKurtosis(2), -1.2483086);
 }
 
+TEST(RHistStats, ComputeSkewnessKurtosisVar0)
+{
+   RHistStats stats(1);
+   stats.Fill(1, RWeight(0.1));
+   ASSERT_EQ(stats.GetNEntries(), 1);
+   EXPECT_EQ(stats.ComputeVariance(), 0);
+   EXPECT_EQ(stats.ComputeSkewness(), 0);
+   EXPECT_EQ(stats.ComputeKurtosis(), 0);
+}
+
 TEST(RHistStats, FillInvalidNumberOfArguments)
 {
    RHistStats stats1(1);