Changeset 15442 for branches/MathNetNumerics-Exploration-2789/Test/UnitTest1.cs

Timestamp:

10/31/17 20:01:02 (6 years ago)

Author:

gkronber

Message:

#2789: added cross-validation and automatic determination of regularization parameter (naive CV is slow!).

Research produced ACM Algorithm 642.f (Fortran) for univariate cubic spline smoothing with generalized cross-validation in O(n) described in HUTCHINSON, M. F., AND DE HOOG, F.R. Smoothing noisy data with spline functions. Numer. Math. 4 7 (1985), 99-106. (TODO: translate to C#)

File:

• branches/MathNetNumerics-Exploration-2789/Test/UnitTest1.cs (modified) (2 diffs)

branches/MathNetNumerics-Exploration-2789/Test/UnitTest1.cs

@@ -1 +1 @@
 using System;
+using System.Diagnostics;
+using System.Linq;
 using HeuristicLab.Algorithms.DataAnalysis.Experimental;
 using HeuristicLab.Algorithms.GeneticAlgorithm;
+using HeuristicLab.Common;
+using HeuristicLab.Problems.DataAnalysis;
 using HeuristicLab.SequentialEngine;
 using Microsoft.VisualStudio.TestTools.UnitTesting;
@@ -22 +26 @@
     }
   }
+
+  [TestClass]
+  public class PenalizedRegressionSplineTests {
+    [TestMethod]
+    public void TestPenalizedRegressionSplinesCrossValidation() {
+      var xs = HeuristicLab.Common.SequenceGenerator.GenerateSteps(-3.5, 3.5, 0.02, includeEnd: true).ToArray();
+      var ys = xs.Select(xi => alglib.normaldistr.normaldistribution(xi)); // 1.0 / (Math.Sqrt(2 * Math.PI) * Math.Exp(-0.5 * xi * xi))).ToArray();
+
+      alglib.hqrndstate state;
+      alglib.hqrndseed(1234, 5678, out state);
+      var ys_noise = ys.Select(yi => yi + alglib.hqrndnormal(state) * 0.01).ToArray();
+
+      double bestRho = -15;
+      double best_loo_rmse = double.PositiveInfinity;
+      var clock = new Stopwatch();
+      clock.Start();
+      int iters = 0;
+      for (int rho = -15; rho <= 15; rho++) {
+        double loo_rmse;
+        double avgTrainRmse;
+        Splines.CalculatePenalizedRegressionSpline(xs, ys_noise, rho, "y", new string[] { "x" }, out avgTrainRmse, out loo_rmse);
+        iters++;
+        Console.WriteLine("{0} {1} {2}", rho, avgTrainRmse, loo_rmse);
+        if (loo_rmse < best_loo_rmse) {
+          best_loo_rmse = loo_rmse;
+          bestRho = rho;
+        }
+      }
+      clock.Stop();
+      Console.WriteLine("Best rho {0}, RMSE (LOO): {1}, ms/run {2}", bestRho, best_loo_rmse, clock.ElapsedMilliseconds / (double)iters);
+    }
+
+    [TestMethod]
+    public void TestReinschSplineCrossValidation() {
+      var xs = HeuristicLab.Common.SequenceGenerator.GenerateSteps(-3.5, 3.5, 0.02, includeEnd: true).ToList();
+      var ys = xs.Select(xi => alglib.normaldistr.normaldistribution(xi)); // 1.0 / (Math.Sqrt(2 * Math.PI) * Math.Exp(-0.5 * xi * xi))).ToArray();
+
+      alglib.hqrndstate state;
+      alglib.hqrndseed(1234, 5678, out state);
+      var ys_noise = ys.Select(yi => yi + alglib.hqrndnormal(state) * 0.01).ToList();
+
+      var tol = ys_noise.StandardDeviation();
+
+      var ds = new Dataset(new string[] { "x", "y" }, new[] { xs, ys_noise });
+      var rows = Enumerable.Range(0, xs.Count);
+
+      var bestTol = double.PositiveInfinity;
+      var bestLooRmse = double.PositiveInfinity;
+      var clock = new Stopwatch();
+      clock.Start();
+      var iters = 0;
+      while (tol > 0.0001) {
+        double loo_rmse;
+        double avgTrainRmse;
+        var model = Splines.CalculateSmoothingSplineReinsch(
+          xs.ToArray(), ys_noise.ToArray(), tol, "y", new string[] { "x" }, out avgTrainRmse, out loo_rmse);
+        var y_pred = model.GetEstimatedValues(ds, rows);
+
+        Console.WriteLine("{0} {1} {2}", tol, avgTrainRmse, loo_rmse);
+        if (loo_rmse < bestLooRmse) {
+          bestLooRmse = loo_rmse;
+          bestTol = tol;
+        }
+        tol *= 0.8;
+        iters++;
+      }
+      clock.Stop();
+      Console.WriteLine("Best tolerance {0}, RMSE (LOO): {1}, ms/run {2}", bestTol, bestLooRmse, clock.ElapsedMilliseconds / (double)iters);
+    }
+
+    [TestMethod]
+    public void TestReinschSplineAutomaticTolerance() {
+      var xs = HeuristicLab.Common.SequenceGenerator.GenerateSteps(-3.5, 3.5, 0.02, includeEnd: true).ToList();
+      var ys = xs.Select(xi => alglib.normaldistr.normaldistribution(xi)); // 1.0 / (Math.Sqrt(2 * Math.PI) * Math.Exp(-0.5 * xi * xi))).ToArray();
+
+      alglib.hqrndstate state;
+      alglib.hqrndseed(1234, 5678, out state);
+      var ys_noise = ys.Select(yi => yi + alglib.hqrndnormal(state) * 0.01).ToList();
+
+      double optTol;
+      double looRMSE;
+      Splines.CalculateSmoothingSplineReinsch(xs.ToArray(), ys_noise.ToArray(), new string[] { "x" }, "y", out optTol, out looRMSE);
+
+      Console.WriteLine("Best tolerance {0}, RMSE (LOO): {1}", optTol, looRMSE);
+
+    }
+  }
 }