Re #7672. Added moment calculations to Kernel/Statistics.

The two methods are getMomentsAboutOrigin and getMomentsAboutMean.
They work for densities and histograms.

Code/Mantid/Framework/Kernel/inc/MantidKernel/Statistics.h

@@ -77,7 +77,13 @@ namespace Mantid
     std::vector<double> getModifiedZscore(const std::vector<TYPE>& data, const bool sorted=false);
     /// Return the R-factors (Rwp) of a diffraction pattern data
     Rfactor MANTID_KERNEL_DLL getRFactor(const std::vector<double>& obsI, const std::vector<double>& calI, const std::vector<double>& obsE);
-  
+
+    /// Return the first n-moments of the supplied data.
+    template<typename TYPE>
+    std::vector<double> getMomentsAboutOrigin(const std::vector<TYPE>& x, const std::vector<TYPE>& y, const int maxMoment=3);
+    /// Return the first n-moments of the supplied data.
+    template<typename TYPE>
+    std::vector<double> getMomentsAboutMean(const std::vector<TYPE>& x, const std::vector<TYPE>& y, const int maxMoment=3);
   } // namespace Kernel
 } // namespace Mantid
 #endif /* STATISTICS_H_ */

Code/Mantid/Framework/Kernel/src/Statistics.cpp

@@ -267,12 +267,138 @@ namespace Mantid
       return rfactor;
     }
 
+    /**
+     * This will calculate the first n-moments (inclusive) about the origin. For example
+     * if maxMoment=2 then this will return 3 values: 0th (total weight), 1st (mean), 2nd (deviation).
+     *
+     * @param x The independent values
+     * @param y The dependent values
+     * @param maxMoment The number of moments to calculate
+     * @returns The first n-moments.
+     */
+    template<typename TYPE>
+    std::vector<double> getMomentsAboutOrigin(const std::vector<TYPE>& x, const std::vector<TYPE>& y, const int maxMoment)
+    {
+      // densities have the same number of x and y
+      bool isDensity(x.size() == y.size());
+
+      // if it isn't a density then check for histogram
+      if ((!isDensity) && (x.size() != y.size()+1))
+      {
+        std::stringstream msg;
+        msg << "length of x (" << x.size() << ") and y (" << y.size() << ")do not match";
+        throw std::out_of_range(msg.str());
+      }
+
+      // initialize a result vector with all zeros
+      std::vector<double> result(maxMoment+1, 0.);
+
+      // as backwards as it sounds, the outer loop should be the points rather than the moments
+
+      // densities are calculated using Newton's method for numerical integration
+      if (isDensity)
+      {
+        for (size_t j = 0; j < x.size()-1; ++j)
+        {
+          const double xVal = .5*static_cast<double>(x[j]+x[j+1]); // reduce item lookup
+          const double xDelta = static_cast<double>(x[j+1]-x[j]);
+          double temp = .5*static_cast<double>(y[j]+y[j+1])*xDelta; // this variable will be (x^n)*y
+          result[0] += temp;
+          for (size_t i = 1; i < result.size(); ++i)
+          {
+            temp *= xVal;
+            result[i] += temp;
+          }
+        }
+      }
+      else
+      {
+        for (size_t j = 0; j < y.size(); ++j)
+        {
+          // central x in histogram
+          const double xVal = .5*static_cast<double>(x[j]+x[j+1]);
+          double temp = static_cast<double>(y[j]);
+          result[0] += temp;
+          for (size_t i = 1; i < result.size(); ++i)
+          {
+            temp *= xVal;
+            result[i] += temp;
+          }
+        }
+      }
+
+      return result;
+    }
+
+    /**
+     * This will calculate the first n-moments (inclusive) about the mean (1st moment). For example
+     * if maxMoment=2 then this will return 3 values: 0th (total weight), 1st (mean), 2nd (deviation).
+     *
+     * @param x The independent values
+     * @param y The dependent values
+     * @param maxMoment The number of moments to calculate
+     * @returns The first n-moments.
+     */
+    template<typename TYPE>
+    std::vector<double> getMomentsAboutMean(const std::vector<TYPE>& x, const std::vector<TYPE>& y, const int maxMoment)
+    {
+      // get the zeroth (integrated value) and first moment (mean)
+      std::vector<double> momentsAboutOrigin = getMomentsAboutOrigin(x, y, 1);
+      const double mean = momentsAboutOrigin[1];
+
+      // initialize a result vector with all zeros
+      std::vector<double> result(maxMoment+1, 0.);
+      result[0] = momentsAboutOrigin[0];
+
+      // escape early if we need to
+      if (maxMoment == 0)
+        return result;
+
+      // as backwards as it sounds, the outer loop should be the points rather than the moments
+
+      // densities are calculated using Newton's method for numerical integration
+      if (x.size() == y.size())
+      {
+        for (size_t j = 0; j < x.size(); ++j)
+        {
+          const double xVal = .5*static_cast<double>(x[j]+x[j+1]) - mean; // change of variables
+          const double xDelta = static_cast<double>(x[j+1]-x[j]);
+          double temp = xVal * .5*static_cast<double>(y[j]+y[j+1])*xDelta; // this variable will be (x^n)*y
+          result[1] += temp;
+          for (size_t i = 2; i < result.size(); ++i)
+          {
+            temp *= xVal;
+            result[i] += temp;
+          }
+        }
+      }
+      else
+      {
+        for (size_t j = 0; j < y.size(); ++j)
+        {
+          // central x in histogram with a change of variables
+          const double xVal = .5*static_cast<double>(x[j]+x[j+1]) - mean;
+          double temp = xVal * static_cast<double>(y[j]);
+          result[1] += temp;
+          for (size_t i = 2; i < result.size(); ++i)
+          {
+            temp *= xVal;
+            result[i] += temp;
+          }
+        }
+      }
+
+      return result;
+    }
 
     // -------------------------- Macro to instantiation concrete types --------------------------------
 #define INSTANTIATE(TYPE) \
     template MANTID_KERNEL_DLL Statistics getStatistics<TYPE> (const vector<TYPE> &, const bool); \
     template MANTID_KERNEL_DLL std::vector<double> getZscore<TYPE> (const vector<TYPE> &, const bool); \
-    template MANTID_KERNEL_DLL std::vector<double> getModifiedZscore<TYPE> (const vector<TYPE> &, const bool);
+    template MANTID_KERNEL_DLL std::vector<double> getModifiedZscore<TYPE> (const vector<TYPE> &, const bool); \
+    template MANTID_KERNEL_DLL std::vector<double> getMomentsAboutOrigin<TYPE> (const std::vector<TYPE>& x, const std::vector<TYPE>& y, const int maxMoment); \
+    template MANTID_KERNEL_DLL std::vector<double> getMomentsAboutMean<TYPE> (const std::vector<TYPE>& x, const std::vector<TYPE>& y, const int maxMoment);
+
 
     // --------------------------- Concrete instantiations ---------------------------------------------
     INSTANTIATE(float);

Code/Mantid/Framework/Kernel/test/StatisticsTest.h

@@ -2,6 +2,7 @@
 #define STATISTICSTEST_H_
 
 #include <cxxtest/TestSuite.h>
+#include <cmath>
 #include <vector>
 #include <string>
 #include "MantidKernel/Statistics.h"
@@ -182,6 +183,88 @@ public:
     TS_ASSERT_THROWS_ANYTHING(getRFactor(obsY, calY, obsE));
   }
 
+  /// Test moment calculations about origin and mean
+  void test_getMoments()
+  {
+    const double mean = 5.;
+    const double sigma = 4.;
+    const double deltaX = .2;
+    const size_t numX = 200;
+    // calculate to have same number of points left and right of function
+    const double offsetX = mean - (.5 * deltaX * static_cast<double>(numX));
+    // variance about origin
+    double expVar = mean*mean+sigma*sigma;
+    // skew about origin
+    double expSkew = mean*mean*mean+3.*mean*sigma*sigma;
+
+    // x-values to try out
+    vector<double> x;
+    for (size_t i = 0; i < numX; ++i)
+      x.push_back(static_cast<double>(i) * deltaX + offsetX);
+
+    // just declare so we can have test of exception handling
+    vector<double> y;
+
+    TS_ASSERT_THROWS(getMomentsAboutOrigin(x, y), std::out_of_range);
+
+    // now calculate the y-values
+    for (size_t i = 0; i < numX; ++i)
+    {
+      double temp = (x[i]-mean)/sigma;
+      y.push_back(exp(-.5*temp*temp)/(sigma * sqrt(2.*M_PI)));
+    }
+
+    // Normal distribution values are taken from the wikipedia page
+    {
+      std::cout << "Normal distribution about origin" << std::endl;
+      vector<double> aboutOrigin = getMomentsAboutOrigin(x, y);
+      TS_ASSERT_EQUALS(aboutOrigin.size(), 4);
+      TS_ASSERT_DELTA(aboutOrigin[0], 1., .0001);
+      TS_ASSERT_DELTA(aboutOrigin[1], mean, .0001);
+      TS_ASSERT_DELTA(aboutOrigin[2], expVar, .001*expVar);
+      TS_ASSERT_DELTA(aboutOrigin[3], expSkew, .001*expSkew);
+
+      std::cout << "Normal distribution about mean" << std::endl;
+      vector<double> aboutMean = getMomentsAboutMean(x, y);
+      TS_ASSERT_EQUALS(aboutMean.size(), 4);
+      TS_ASSERT_DELTA(aboutMean[0], 1., .0001);
+      TS_ASSERT_DELTA(aboutMean[1], 0., .0001);
+      TS_ASSERT_DELTA(aboutMean[2], sigma*sigma, .001*expVar);
+      TS_ASSERT_DELTA(aboutMean[3], 0., .0001*expSkew);
+    }
+
+    // Now a gaussian function as a histogram
+    y.clear();
+    for (size_t i = 0; i < numX-1; ++i) // one less y than x makes it a histogram
+    {
+      double templeft = (x[i]-mean)/sigma;
+      templeft = exp(-.5*templeft*templeft)/(sigma * sqrt(2.*M_PI));
+      double tempright = (x[i+1]-mean)/sigma;
+      tempright = exp(-.5*tempright*tempright)/(sigma * sqrt(2.*M_PI));
+      y.push_back(.5*deltaX*(templeft+tempright));
+//      std::cout << i << ":\t" << x[i] << "\t" << y[i] << std::endl;
+    }
+
+    // Normal distribution values are taken from the wikipedia page
+    {
+      std::cout << "Normal distribution about origin" << std::endl;
+      vector<double> aboutOrigin = getMomentsAboutOrigin(x, y);
+      TS_ASSERT_EQUALS(aboutOrigin.size(), 4);
+      TS_ASSERT_DELTA(aboutOrigin[0], 1., .0001);
+      TS_ASSERT_DELTA(aboutOrigin[1], mean, .0001);
+      TS_ASSERT_DELTA(aboutOrigin[2], expVar, .001*expVar);
+      TS_ASSERT_DELTA(aboutOrigin[3], expSkew, .001*expSkew);
+
+      std::cout << "Normal distribution about mean" << std::endl;
+      vector<double> aboutMean = getMomentsAboutMean(x, y);
+      TS_ASSERT_EQUALS(aboutMean.size(), 4);
+      TS_ASSERT_DELTA(aboutMean[0], 1., .0001);
+      TS_ASSERT_DELTA(aboutMean[1], 0., .0001);
+      TS_ASSERT_DELTA(aboutMean[2], sigma*sigma, .001*expVar);
+      TS_ASSERT_DELTA(aboutMean[3], 0., .0001*expSkew);
+    }
+  }
+
 };
 
 #endif // STATISTICSTEST_H_