Changeset 15450

Timestamp:

11/06/17 13:12:41 (7 years ago)

Author:

gkronber

Message:

#2789 more tests with CV and automatic determination of smoothing parameter

Location:

branches/MathNetNumerics-Exploration-2789

Files:

• HeuristicLab.Algorithms.DataAnalysis.Experimental/GAM.cs (modified) (4 diffs)
• HeuristicLab.Algorithms.DataAnalysis.Experimental/Splines.cs (modified) (7 diffs)
• Main/Main.csproj (modified) (1 diff)
• Main/Program.cs (modified) (1 diff)
• Test/Test.csproj (modified) (2 diffs)
• Test/UnitTest1.cs (modified) (1 diff)

branches/MathNetNumerics-Exploration-2789/HeuristicLab.Algorithms.DataAnalysis.Experimental/GAM.cs

@@ -111 +111 @@
         nTerms += inputVars.Combinations(i).Count();
       }
+
       IRegressionModel[] f = new IRegressionModel[nTerms];
       for (int i = 0; i < f.Length; i++) {
@@ -130 +131 @@
       string[] terms = new string[f.Length];
       Results.Add(new Result("RSS Values", typeof(DoubleMatrix)));
+
+      var combinations = new List<string[]>();
+      for(int i=1;i<=maxInteractions;i++)
+        combinations.AddRange(HeuristicLab.Common.EnumerableExtensions.Combinations(inputVars, i).Select(c => c.ToArray()));
+      // combinations.Add(new string[] { "X1", "X2" });
+      // combinations.Add(new string[] { "X3", "X4" });
+      // combinations.Add(new string[] { "X5", "X6" });
+      // combinations.Add(new string[] { "X1", "X7", "X9" });
+      // combinations.Add(new string[] { "X3", "X6", "X10" });
+
+
 
       // until convergence
@@ -144 +156 @@
         //}
 
-        for (int interaction = 1; interaction <= maxInteractions; interaction++) {
-          var selectedVars = HeuristicLab.Common.EnumerableExtensions.Combinations(inputVars, interaction);
-
-          foreach (var element in selectedVars) {
-            var res = CalculateResiduals(problemData, f, j, avgY, problemData.TrainingIndices);
-            rss[j] = res.Variance();
-            terms[j] = string.Format("f({0})", string.Join(",", element));
-            f[j] = RegressSpline(problemData, element.ToArray(), res, lambda);
-            j++;
-          }
+
+
+        foreach (var element in combinations) {
+          var res = CalculateResiduals(problemData, f, j, avgY, problemData.TrainingIndices);
+          rss[j] = res.Variance();
+          terms[j] = string.Format("f({0})", string.Join(",", element));
+          f[j] = RegressSpline(problemData, element.ToArray(), res, lambda);
+          j++;
         }
 
@@ -218 +228 @@
           product = product.Zip(problemData.Dataset.GetDoubleValues(inputVars[i], problemData.TrainingIndices), (pi, vi) => pi * vi).ToArray();
         }
-        CubicSplineGCV.CubGcvReport report;
-        return CubicSplineGCV.CalculateCubicSpline(
-          product,
-          (double[])target.Clone(), 
-          problemData.TargetVariable, inputVars, out report
-          );
+        // CubicSplineGCV.CubGcvReport report;
+        // return CubicSplineGCV.CalculateCubicSpline(
+        //   product,
+        //   (double[])target.Clone(), 
+        //   problemData.TargetVariable, inputVars, out report
+        //   );
+
+        double optTolerance; double cvRMSE;
+        // find tolerance
+        // var ensemble = Splines.CalculateSmoothingSplineReinsch(product, (double[])target.Clone(), inputVars, problemData.TargetVariable, out optTolerance, out cvRMSE);
+        // // train on whole data
+        // return Splines.CalculateSmoothingSplineReinsch(product, (double[])target.Clone(), inputVars, optTolerance, product.Length - 1, problemData.TargetVariable); 
+
+
+        // find tolerance
+        var bestLambda = -5.0;
+        double bestCVRMSE = double.PositiveInfinity;
+        double avgTrainRMSE = double.PositiveInfinity;
+        for (double curLambda = -5.0; curLambda <= 6.0; curLambda += 1.0) {
+          var ensemble = Splines.CalculatePenalizedRegressionSpline(product, (double[])target.Clone(), curLambda, problemData.TargetVariable,  inputVars, out avgTrainRMSE, out cvRMSE);
+          Console.Write("{0} {1} {2}", curLambda, avgTrainRMSE, cvRMSE);
+          if (bestCVRMSE > cvRMSE) {
+            Console.Write(" *");
+            bestCVRMSE = cvRMSE;
+            bestLambda = curLambda;
+          }
+          Console.WriteLine();
+        }
+        // train on whole data
+       return Splines.CalculatePenalizedRegressionSpline(product, (double[])target.Clone(), bestLambda, problemData.TargetVariable, inputVars, out avgTrainRMSE, out cvRMSE);
+
       } else return new ConstantModel(target.Average(), problemData.TargetVariable);
     }

branches/MathNetNumerics-Exploration-2789/HeuristicLab.Algorithms.DataAnalysis.Experimental/Splines.cs

@@ -22 +22 @@
 using System;
 using System.Collections.Generic;
+using System.Diagnostics;
 using System.Linq;
 using System.Threading;
@@ -156 +157 @@
           case "CUBGCV": {
               CubicSplineGCV.CubGcvReport report;
-              var model = 
-                CubicSplineGCV.CalculateCubicSpline(x, y, 
+              var model =
+                CubicSplineGCV.CalculateCubicSpline(x, y,
                 Problem.ProblemData.TargetVariable, inputVars, out report);
               var targetVar = Problem.ProblemData.TargetVariable;
@@ -165 +166 @@
               Results.Add(new Result("Estimated error variance", new DoubleValue(report.estimatedErrorVariance)));
               Results.Add(new Result("Estimated RSS degrees of freedom", new DoubleValue(report.estimatedRSSDegreesOfFreedom)));
-              Results.Add(new Result("Estimated treu mean squared error at data points", new DoubleValue(report.estimatedTrueMeanSquaredErrorAtDataPoints)));
+              Results.Add(new Result("Estimated true mean squared error at data points", new DoubleValue(report.estimatedTrueMeanSquaredErrorAtDataPoints)));
               Results.Add(new Result("Mean square of DF", new DoubleValue(report.meanSquareOfDf)));
               Results.Add(new Result("Mean square residual", new DoubleValue(report.meanSquareResiudal)));
@@ -204 +205 @@
     public static IRegressionEnsembleModel CalculatePenalizedRegressionSpline(double[] x, double[] y, double lambda, string targetVar, string[] inputVars,
       out double avgTrainRMSE,
-      out double loocvRMSE) {
+      out double cvRMSE) {
       int info;
       alglib.spline1dinterpolant interpolant;
@@ -211 +212 @@
       n = (int)Math.Max(50, 3 * Math.Sqrt(n));
 
-      double[] x_loo = new double[x.Length - 1];
-      double[] y_loo = new double[y.Length - 1];
+      int folds = 10;
+      int foldSize = x.Length <= 30 ? 1 : x.Length / folds;
+      if (foldSize == 1) folds = x.Length;
+
+      double[] x_loo = new double[(folds - 1) * foldSize];
+      double[] y_loo = new double[(folds - 1) * foldSize];
 
       var models = new List<IRegressionModel>();
-      var sumTrainRmse = 0.0;
-      var sse = 0.0;
-
-      for (int i = 0; i < x.Length; i++) {
-        Array.Copy(x, 0, x_loo, 0, i);
-        Array.Copy(x, i + 1, x_loo, i, x.Length - (i + 1));
-        Array.Copy(y, 0, y_loo, 0, i);
-        Array.Copy(y, i + 1, y_loo, i, y.Length - (i + 1));
+      var sumTrainSE = 0.0;
+      int nTrain = 0;
+      var sumTestSE = 0.0;
+      int nTest = 0;
+
+
+      for (int i = 0; i < folds; i++) {
+        if (i < folds - 1) {
+          //Console.WriteLine("Copy from {0} to {1} n = {2}", (i + 1) * foldSize, i * foldSize, (folds - i - 1) * foldSize);
+          Array.Copy(x, (i + 1) * foldSize, x_loo, i * foldSize, (folds - i - 1) * foldSize);
+          Array.Copy(y, (i + 1) * foldSize, y_loo, i * foldSize, (folds - i - 1) * foldSize);
+        }
         alglib.spline1dfitpenalized(x_loo, y_loo, n, lambda, out info, out interpolant, out rep);
-        var y_pred = alglib.spline1dcalc(interpolant, x[i]);
-        double res = y[i] - y_pred;
-        sse += res * res;
-        sumTrainRmse += rep.rmserror;
-
-        var model = new Spline1dModel(interpolant, targetVar, inputVars);
-        models.Add(model);
-      }
-
-      loocvRMSE = Math.Sqrt(sse / x.Length);
-      avgTrainRMSE = sumTrainRmse / x.Length;
+        Debug.Assert(y_loo.All(yi => yi != 0.0));
+        Debug.Assert(x_loo.All(xi => xi != 0.0));
+
+        // training error
+        for (int j = 0; j < x_loo.Length; j++) {
+          var y_pred = alglib.spline1dcalc(interpolant, x[i]);
+          double res = y[j] - y_pred;
+          sumTrainSE += res * res;
+          nTrain++;
+        }
+        // test
+        for (int j = i * foldSize; j < (i + 1) * foldSize; j++) {
+          var y_pred = alglib.spline1dcalc(interpolant, x[i]);
+          double res = y[j] - y_pred;
+          sumTestSE += res * res;
+          nTest++;
+        }
+
+        models.Add(new Spline1dModel(interpolant, targetVar, inputVars));
+
+        if (i < folds - 1) {
+          //Console.WriteLine("Copy from {0} to {1} n = {2}", i * foldSize, i * foldSize, foldSize);
+          // add current fold for next iteration
+          Array.Copy(x, i * foldSize, x_loo, i * foldSize, foldSize);
+          Array.Copy(y, i * foldSize, y_loo, i * foldSize, foldSize);
+        }
+      }
+
+      cvRMSE = Math.Sqrt(sumTestSE / nTest);
+      avgTrainRMSE = Math.Sqrt(sumTrainSE / nTrain);
 
       return new RegressionEnsembleModel(models);
@@ -352 +380 @@
     public static IRegressionEnsembleModel CalculateSmoothingSplineReinsch(double[] x, double[] y, double tol, string targetVar, string[] inputVars,
       out double avgTrainRMSE,
-      out double loocvRMSE) {
+      out double cvRMSE) {
       int info;
 
-      double[] x_loo = new double[x.Length - 1];
-      double[] y_loo = new double[y.Length - 1];
+      int folds = 10;
+      int foldSize = x.Length <= 30 ? 1 : x.Length / folds;
+      if (foldSize == 1) folds = x.Length;
+
+      double[] x_loo = new double[(folds - 1) * foldSize];
+      double[] y_loo = new double[(folds - 1) * foldSize];
 
       var models = new List<IRegressionModel>();
-      var sumTrainRmse = 0.0;
-      var sse = 0.0;
-
-      for (int i = 0; i < x.Length; i++) {
-        Array.Copy(x, i + 1, x_loo, i, x.Length - (i + 1));
-        Array.Copy(y, i + 1, y_loo, i, y.Length - (i + 1));
-
-        var model = CalculateSmoothingSplineReinsch(x_loo, y_loo, inputVars, tol, x_loo.Length + 1, targetVar);
-
-        var y_pred = model.GetEstimatedValue(x[i]);
-        double res = y[i] - y_pred;
-        sse += res * res;
+      var sumTrainSE = 0.0;
+      int nTrain = 0;
+      var sumTestSE = 0.0;
+      int nTest = 0;
+
+
+      for (int i = 0; i < folds; i++) {
+        if (i < folds - 1) {
+          Array.Copy(x, (i + 1) * foldSize, x_loo, i * foldSize, x.Length - ((i + 1) * foldSize) - 1);
+          Array.Copy(y, (i + 1) * foldSize, y_loo, i * foldSize, y.Length - ((i + 1) * foldSize) - 1);
+        }
+        var model = CalculateSmoothingSplineReinsch(x_loo, y_loo, inputVars, tol, ((folds - 1) * foldSize) - 1, targetVar);
+
+        // training error
+        for(int j = 0;j<x_loo.Length;j++) {
+          var y_pred = model.GetEstimatedValue(x[j]);
+          double res = y[j] - y_pred;
+          sumTrainSE += res * res;
+          nTrain++;
+        }
+        // test
+        for (int j = i * foldSize; j < (i + 1) * foldSize; j++) {
+          var y_pred = model.GetEstimatedValue(x[j]);
+          double res = y[j] - y_pred;
+          sumTestSE += res * res;
+          nTest++;
+        }
 
         models.Add(model);
 
-        // add current row for next iteration
-        if (i < x_loo.Length) {
-          x_loo[i] = x[i];
-          y_loo[i] = y[i];
-        }
-      }
-
-      loocvRMSE = Math.Sqrt(sse / x.Length);
-      avgTrainRMSE = sumTrainRmse / x.Length;
+        if (i < folds - 1) {
+          // add current fold for next iteration
+          Array.Copy(x, i * foldSize, x_loo, i * foldSize, foldSize);
+          Array.Copy(y, i * foldSize, y_loo, i * foldSize, foldSize);
+        }
+      }
+
+      cvRMSE = Math.Sqrt(sumTestSE / nTest);
+      avgTrainRMSE = Math.Sqrt(sumTrainSE / nTrain);
 
       return new RegressionEnsembleModel(models);
@@ -388 +435 @@
 
 
-    // automatically determines tolerance using loo cv
+    // automatically determines tolerance using loo cv, SLOW!
     public static IRegressionModel CalculateSmoothingSplineReinsch(double[] x, double[] y, string[] inputVars, string targetVar,
-      out double optTolerance, out double looRMSE) {
+      out double optTolerance, out double cvRMSE) {
       var maxTolerance = y.StandardDeviation();
       var minTolerance = maxTolerance * 1E-6;
       optTolerance = maxTolerance;
       var tol = maxTolerance;
-      looRMSE = double.PositiveInfinity;
+      cvRMSE = double.PositiveInfinity;
       IRegressionModel bestModel = new ConstantModel(y.Average(), targetVar);
       while (tol > minTolerance) {
-        double loo_rmse;
+        double curCVRMSE;
         double avgTrainRmse;
         var model = Splines.CalculateSmoothingSplineReinsch(
-          x, y, tol, targetVar, inputVars, out avgTrainRmse, out loo_rmse);
-
-        if (loo_rmse < looRMSE) {
-          looRMSE = loo_rmse;
+          x, y, tol, targetVar, inputVars, out avgTrainRmse, out curCVRMSE);
+
+        if (curCVRMSE < cvRMSE) {
+          cvRMSE = curCVRMSE;
           optTolerance = tol;
           bestModel = model;

branches/MathNetNumerics-Exploration-2789/Main/Main.csproj

@@ -38 +38 @@
       <HintPath>..\..\..\trunk\sources\bin\ALGLIB-3.7.0.dll</HintPath>
     </Reference>
+    <Reference Include="HeuristicLab.Algorithms.DataAnalysis-3.4, Version=3.4.0.0, Culture=neutral, PublicKeyToken=ba48961d6f65dcec, processorArchitecture=MSIL">
+      <SpecificVersion>False</SpecificVersion>
+      <HintPath>..\..\..\trunk\sources\bin\HeuristicLab.Algorithms.DataAnalysis-3.4.dll</HintPath>
+    </Reference>
+    <Reference Include="HeuristicLab.Collections-3.3, Version=3.3.0.0, Culture=neutral, PublicKeyToken=ba48961d6f65dcec" />
     <Reference Include="HeuristicLab.Common-3.3">
       <HintPath>..\..\..\trunk\sources\bin\HeuristicLab.Common-3.3.dll</HintPath>
     </Reference>
     <Reference Include="HeuristicLab.Core-3.3, Version=3.3.0.0, Culture=neutral, PublicKeyToken=ba48961d6f65dcec" />
+    <Reference Include="HeuristicLab.Optimization-3.3, Version=3.3.0.0, Culture=neutral, PublicKeyToken=ba48961d6f65dcec" />
     <Reference Include="HeuristicLab.Problems.DataAnalysis-3.4, Version=3.4.0.0, Culture=neutral, PublicKeyToken=ba48961d6f65dcec, processorArchitecture=MSIL">
       <SpecificVersion>False</SpecificVersion>
       <HintPath>..\..\..\trunk\sources\bin\HeuristicLab.Problems.DataAnalysis-3.4.dll</HintPath>
+    </Reference>
+    <Reference Include="HeuristicLab.Problems.Instances-3.3, Version=3.3.0.0, Culture=neutral, PublicKeyToken=ba48961d6f65dcec, processorArchitecture=MSIL">
+      <SpecificVersion>False</SpecificVersion>
+      <HintPath>..\..\..\trunk\sources\bin\HeuristicLab.Problems.Instances-3.3.dll</HintPath>
+    </Reference>
+    <Reference Include="HeuristicLab.Problems.Instances.DataAnalysis-3.3">
+      <HintPath>..\..\..\trunk\sources\bin\HeuristicLab.Problems.Instances.DataAnalysis-3.3.dll</HintPath>
     </Reference>
     <Reference Include="System" />

branches/MathNetNumerics-Exploration-2789/Main/Program.cs

@@ -12 +12 @@
   class Program {
     static void Main(string[] args) {
-      var xs = HeuristicLab.Common.SequenceGenerator.GenerateSteps(-3.5, 3.5, 0.1, includeEnd: true).ToList();
-      var ys = xs.Select(xi => 1.0 / Math.Sqrt(2 * Math.PI) * Math.Exp(-0.5 * xi * xi)).ToArray(); // 1.0 / (Math.Sqrt(2 * Math.PI) * Math.Exp(-0.5 * xi * xi))).ToArray();
+      var provider = new HeuristicLab.Problems.Instances.DataAnalysis.VariousInstanceProvider();
+      var problemData = provider.LoadData(provider.GetDataDescriptors().First(dd => dd.Name.Contains("Poly")));
+      // var provider = new HeuristicLab.Problems.Instances.DataAnalysis.RegressionRealWorldInstanceProvider();
+      // var problemData = provider.LoadData(provider.GetDataDescriptors().First(dd => dd.Name.Contains("Chem")));
 
-      int n = xs.Count();
-      alglib.hqrndstate state;
-      alglib.hqrndseed(1234, 5678, out state);
-      var ys_noise = ys.Select(yi => yi + alglib.hqrndnormal(state) * 0.1).ToList();
+      var gam = new GAM();
+      gam.MaxIterations = 10;
+      gam.MaxInteractions = 3;
+      gam.Problem.ProblemData = problemData;
+      gam.Start();
 
-      CubicSplineGCV.CubGcvReport report;
-      var model = CubicSplineGCV.CalculateCubicSpline(
-        xs.ToArray(), ys_noise.ToArray(), "y", new string[] { "x" }, out report);
+      var solution = (IRegressionSolution)gam.Results["Ensemble solution"].Value;
 
-      Console.WriteLine("Smoothing Parameter (= RHO/(RHO + 1) {0}", report.smoothingParameter);
-      Console.WriteLine("Estimated DOF of RSS {0}", report.estimatedRSSDegreesOfFreedom);
-      Console.WriteLine("GCV {0}", report.generalizedCrossValidation);
-      Console.WriteLine("Mean squared residual {0}", report.meanSquareResiudal);
-      Console.WriteLine("Estimate of true MSE at data points {0}", report.estimatedTrueMeanSquaredErrorAtDataPoints);
-      Console.WriteLine("Estimate of error variance {0}", report.estimatedErrorVariance);
-      Console.WriteLine("Mean square value of DF(I) {0}", report.meanSquareOfDf);
-
-      OnlineCalculatorError error;
-      var ys_smoothed = xs.Select(xi => model.GetEstimatedValue(xi)).ToArray();
-      var mse = OnlineMeanSquaredErrorCalculator.Calculate(ys, ys_smoothed, out error);
-      Console.WriteLine("MSE(ys, ys_smooth) = {0}", mse);
-      mse = OnlineMeanSquaredErrorCalculator.Calculate(ys, ys_noise, out error);
-      Console.WriteLine("MSE(ys, ys_noise) = {0}", mse);
-
-      Thread.CurrentThread.CurrentCulture = CultureInfo.InvariantCulture;
-
-      for (int i = 0; i < n; i++) {
-        Console.WriteLine("{0}\t{1}\t{2}\t{3}\t{4}\t{5}",
-          xs[i], ys[i], ys_smoothed[i], ys_noise[i], ys_smoothed[i] + 1.96 * report.se[i], ys_smoothed[i] - 1.96 * report.se[i]);
-      }
+      Console.WriteLine("RMSE (train) {0}", solution.TrainingRootMeanSquaredError);
+      Console.WriteLine("RMSE (test) {0}", solution.TestRootMeanSquaredError);
     }
   }

branches/MathNetNumerics-Exploration-2789/Test/Test.csproj

@@ -62 +62 @@
       <HintPath>..\..\..\trunk\sources\bin\HeuristicLab.Algorithms.GeneticAlgorithm-3.3.dll</HintPath>
     </Reference>
+    <Reference Include="HeuristicLab.Collections-3.3, Version=3.3.0.0, Culture=neutral, PublicKeyToken=ba48961d6f65dcec" />
     <Reference Include="HeuristicLab.Common-3.3">
       <HintPath>..\..\..\trunk\sources\bin\HeuristicLab.Common-3.3.dll</HintPath>
@@ -90 +91 @@
     <Reference Include="HeuristicLab.Problems.Instances-3.3">
       <HintPath>..\..\..\trunk\sources\bin\HeuristicLab.Problems.Instances-3.3.dll</HintPath>
+    </Reference>
+    <Reference Include="HeuristicLab.Problems.Instances.DataAnalysis-3.3">
+      <HintPath>..\..\..\trunk\sources\bin\HeuristicLab.Problems.Instances.DataAnalysis-3.3.dll</HintPath>
     </Reference>
     <Reference Include="HeuristicLab.SequentialEngine-3.3">

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

@@ -110 +110 @@
 
       Console.WriteLine("Best tolerance {0}, RMSE (LOO): {1}", optTol, looRMSE);
+    }
 
+    [TestMethod]
+    public void TestGAM() {
+      var provider = new HeuristicLab.Problems.Instances.DataAnalysis.VariousInstanceProvider();
+      var problemData = provider.LoadData(provider.GetDataDescriptors().First(dd => dd.Name.Contains("Poly")));
+      // var provider = new HeuristicLab.Problems.Instances.DataAnalysis.RegressionRealWorldInstanceProvider();
+      // var problemData = provider.LoadData(provider.GetDataDescriptors().First(dd => dd.Name.Contains("Chem")));
+
+      var gam = new GAM();
+      gam.MaxIterations = 10;
+      gam.MaxInteractions = 3;
+      gam.Problem.ProblemData = problemData;
+      gam.Start();
+
+      var solution = (IRegressionSolution)gam.Results["Ensemble solution"].Value;
+
+      Console.WriteLine("RMSE (train) {0}", solution.TrainingRootMeanSquaredError);
+      Console.WriteLine("RMSE (test) {0}", solution.TestRootMeanSquaredError);
     }
-      
+
   }
 }