* Add check and unit test for negative weights.

* Optimise polynomial construction.

Author: rca22

src/Numerics.Tests/StatisticsTests/DescriptiveStatisticsTests.cs

@@ -526,6 +526,12 @@ public void ShortSequences()
             Assert.That(stats9.Kurtosis, Is.Not.NaN);
         }
 
+        [Test]
+        public void NegativeWeightsThrow()
+        {
+            Assert.That(() => new DescriptiveStatistics(new[] { Tuple.Create(-1.0, 1.0) }), Throws.TypeOf<ArgumentOutOfRangeException>());
+        }
+
         [Test]
         public void ZeroVarianceSequence()
         {

src/Numerics.Tests/StatisticsTests/RunningWeightedStatisticsTests.cs

@@ -135,6 +135,15 @@ public void KurtosisConsistentWithR_e1071(string dataSet, double kurtosisType1,
             Assert.That(stats.PopulationKurtosis, Is.EqualTo(kurtosisType1).Within(1e-6), "PopulationKurtosis");
         }
 
+        [Test]
+        public void NegativeWeightsThrow()
+        {
+            Assert.That(() => new RunningWeightedStatistics(new[] { System.Tuple.Create(-1.0, 1.0) }), Throws.TypeOf<System.ArgumentOutOfRangeException>());
+            var stats0 = new RunningWeightedStatistics(new System.Tuple<double,double>[0]);
+            Assert.That(() =>stats0.Push(-1.0, 1.0), Throws.TypeOf<System.ArgumentOutOfRangeException>());
+            Assert.That(() => stats0.PushRange(new[] { System.Tuple.Create(-1.0, 1.0) }), Throws.TypeOf<System.ArgumentOutOfRangeException>());
+        }
+
         [Test]
         public void ShortSequences()
         {

src/Numerics/Statistics/DescriptiveStatistics.cs

@@ -475,7 +475,7 @@ void SetStatistics(double mean, double variance, double skewness, double kurtosi
         }
 
         /// <summary>
-        /// Computes descriptive statistics from a stream of data values.
+        /// Computes descriptive statistics from a stream of data values and weights.
         /// </summary>
         /// <param name="data">A sequence of datapoints.</param>
         void Compute(IEnumerable<Tuple<double, double>> data)
@@ -496,7 +496,11 @@ void Compute(IEnumerable<Tuple<double, double>> data)
 
             foreach (var (w, xi) in data)
             {
-                if (w > 0)
+                if (w < 0)
+                {
+                    throw new ArgumentOutOfRangeException(nameof(data), w, "Expected non-negative weighting of sample");
+                }
+                else if (w > 0)
                 {
                     ++n;
                     double delta = xi - mean;
@@ -537,7 +541,7 @@ void Compute(IEnumerable<Tuple<double, double>> data)
         }
 
         /// <summary>
-        /// Computes descriptive statistics from a stream of data values.
+        /// Computes descriptive statistics from a stream of data values and weights.
         /// </summary>
         /// <param name="data">A sequence of datapoints.</param>
         void ComputeDecimal(IEnumerable<Tuple<double, double>> data)
@@ -558,7 +562,11 @@ void ComputeDecimal(IEnumerable<Tuple<double, double>> data)
 
             foreach (var (w, x) in data)
             {
-                if (w > 0)
+                if (w < 0)
+                {
+                    throw new ArgumentOutOfRangeException(nameof(data), w, "Expected non-negative weighting of sample");
+                }
+                else if (w > 0)
                 {
 
                     decimal xi = (decimal)x;
@@ -643,10 +651,12 @@ void SetStatisticsWeighted(double mean, double variance, double skewness, double
 
                     if (n > 3)
                     {
-                        // common denominator
-                        double poly = w1 * w1 * w1 * w1 - 6.0 * w1 * w1 * w2 + 8.0 * w1 * w3 + 3.0 * w2 * w2 - 6.0 * w4;
-                        double a = w1 * w1 * w1 * w1 - 4.0 * w1 * w3 + 3.0 * w2 * w2;
-                        double b = 3.0 * (w1 * w1 * w1 * w1 - 2.0 * w1 * w1 * w2 + 4.0 * w1 * w3 - 3.0 * w2 * w2);
+                        double p2 = w1 * w1;
+                        double p4 = p2 * p2;
+                        double w2p2 = w2 * w2;
+                        double poly = p4 - 6.0 * p2 * w2 + 8.0 * w1 * w3 + 3.0 * w2p2 - 6.0 * w4;
+                        double a = p4 - 4.0 * w1 * w3 + 3.0 * w2p2;
+                        double b = 3.0 * (p4 - 2.0 * p2 * w2 + 4.0 * w1 * w3 - 3.0 * w2p2);
                         Kurtosis = (a * w1 * kurtosis / (variance * variance) - b) * (den / (w1 * poly));
                     }
                 }

src/Numerics/Statistics/RunningWeightedStatistics.cs

@@ -210,9 +210,12 @@ public double Kurtosis
                     return double.NaN;
                 else
                 {
-                    double poly = _w1 * _w1 * _w1 * _w1 - 6.0 * _w1 * _w1 * _w2 + 8.0 * _w1 * _w3 + 3.0 * _w2 * _w2 - 6.0 * _w4;
-                    double a = _w1 * _w1 * _w1 * _w1 - 4.0 * _w1 * _w3 + 3.0 * _w2 * _w2;
-                    double b = 3.0 * (_w1 * _w1 * _w1 * _w1 - 2.0 * _w1 * _w1 * _w2 + 4.0 * _w1 * _w3 - 3.0 * _w2 * _w2);
+                    double p2 = _w1 * _w1;
+                    double p4 = p2 * p2;
+                    double w2p2 = _w2 * _w2;
+                    double poly = p4 - 6.0 * p2 * _w2 + 8.0 * _w1 * _w3 + 3.0 * w2p2 - 6.0 * _w4;
+                    double a = p4 - 4.0 * _w1 * _w3 + 3.0 * w2p2;
+                    double b = 3.0 * (p4 - 2.0 * p2 * _w2 + 4.0 * _w1 * _w3 - 3.0 * w2p2);
                     return (a * _w1 * _m4 / (_m2 * _m2) - b) * (_den / (_w1 * poly));
                 }
             }
@@ -238,7 +241,7 @@ public void Push(double weight, double value)
             if (weight == 0.0)
                 return;
             if (weight < 0.0)
-                throw new ArgumentException("Expected non-negative weighting of sample", nameof(weight));
+                throw new ArgumentOutOfRangeException(nameof(weight), weight, "Expected non-negative weighting of sample");
 
             _n++;
             double prevW = _w1;