Refs #21752 weighted Zscores

Framework/Crystal/src/PeakStatisticsTools.cpp

@@ -67,7 +67,8 @@ UniqueReflection UniqueReflection::removeOutliers(double sigmaCritical) const {
 
   if (m_peaks.size() > 2) {
     auto intensities = getIntensities();
-    auto zScores = Kernel::getZscore(intensities);
+    auto sigmas = getSigmas();
+    auto zScores = Kernel::getWeightedZscore(intensities, sigmas);
 
     for (size_t i = 0; i < zScores.size(); ++i) {
       if (zScores[i] <= sigmaCritical) {

Framework/Crystal/src/SortHKL.cpp

@@ -127,12 +127,13 @@ void SortHKL::exec() {
       counter++;
       auto wavelengths = unique.second.getWavelengths();
       auto intensities = unique.second.getIntensities();
+      auto sigmas = unique.second.getSigmas();
       std::cout << "HKL " << unique.second.getHKL() << "\n";
       std::cout << "Intensities ";
       for (const auto &e : intensities)
         std::cout << e << "  ";
       std::cout << "\n";
-      auto zScores = Kernel::getZscore(intensities);
+      auto zScores = Kernel::getWeightedZscore(intensities, sigmas);
       if (zScores.size() > maxPeaks)
         maxPeaks = zScores.size();
       // Possibly remove outliers.

Framework/Kernel/inc/MantidKernel/Statistics.h

@@ -93,6 +93,8 @@ Statistics getStatistics(const std::vector<TYPE> &data,
 /// Return the Z score values for a dataset
 template <typename TYPE>
 std::vector<double> getZscore(const std::vector<TYPE> &data);
+template <typename TYPE>
+std::vector<double> getWeightedZscore(const std::vector<TYPE> &data,const std::vector<TYPE> &weights);
 /// Return the modified Z score values for a dataset
 template <typename TYPE>
 std::vector<double> getModifiedZscore(const std::vector<TYPE> &data,

Framework/Kernel/src/Statistics.cpp

@@ -106,6 +106,38 @@ std::vector<double> getZscore(const vector<TYPE> &data) {
   }
   return Zscore;
 }
+/**
+ * There are enough special cases in determining the Z score where it useful to
+ * put it in a single function.
+ */
+template <typename TYPE>
+std::vector<double> getWeightedZscore(const vector<TYPE> &data, const vector<TYPE> &weights) {
+  if (data.size() < 3) {
+    std::vector<double> Zscore(data.size(), 0.);
+    return Zscore;
+  }
+  std::vector<double> Zscore;
+  Statistics stats = getStatistics(data);
+  if (stats.standard_deviation == 0.) {
+    std::vector<double> Zscore(data.size(), 0.);
+    return Zscore;
+  }
+  double sumWeights = 0.0;
+  double sumWeightedData = 0.0;
+  double weightedVariance = 0.0;
+  for (size_t it = 0; it != data.size(); ++it) {
+    sumWeights += static_cast<double>(weights[it]);
+    sumWeightedData += static_cast<double>(weights[it] * data[it]);
+  }
+  double weightedMean = sumWeightedData/sumWeights;
+  for (size_t it = 0; it != data.size(); ++it) {
+    weightedVariance += std::pow(data[it] - weightedMean, 2) * std::pow(weights[it]/sumWeights, 2);
+  }
+  for (auto it = data.cbegin(); it != data.cend(); ++it) {
+    Zscore.push_back(fabs((*it - weightedMean) / std::sqrt(weightedVariance)));
+  }
+  return Zscore;
+}
 /**
  * There are enough special cases in determining the modified Z score where it
  * useful to
@@ -406,6 +438,8 @@ std::vector<double> getMomentsAboutMean(const std::vector<TYPE> &x,
   getStatistics<TYPE>(const vector<TYPE> &, const unsigned int);               \
   template MANTID_KERNEL_DLL std::vector<double> getZscore<TYPE>(              \
       const vector<TYPE> &);                                                   \
+  template MANTID_KERNEL_DLL std::vector<double> getWeightedZscore<TYPE>(              \
+      const vector<TYPE> &, const vector<TYPE> &);                                                   \
   template MANTID_KERNEL_DLL std::vector<double> getModifiedZscore<TYPE>(      \
       const vector<TYPE> &, const bool);                                       \
   template MANTID_KERNEL_DLL std::vector<double> getMomentsAboutOrigin<TYPE>(  \