Changeset 15365

Timestamp:

09/17/17 18:06:27 (7 years ago)

Author:

gkronber

Message:

#2789 experiments with spline smoothing and additive models

Location:

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

Files:

• GAM.cs (modified) (3 diffs)
• Splines.cs (modified) (15 diffs)

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

@@ -109 +109 @@
 
         rmseRow.Values.Add(CalculateResiduals(problemData, f, -1, avgY, problemData.TrainingIndices).StandardDeviation()); // -1 index to use all predictors
-        rmseRowTest.Values.Add(CalculateResiduals(problemData, f, -1, avgY, problemData.TestIndices).StandardDeviation()); 
+        rmseRowTest.Values.Add(CalculateResiduals(problemData, f, -1, avgY, problemData.TestIndices).StandardDeviation());
+
+        if (cancellationToken.IsCancellationRequested) break;
       }
 
@@ -135 +137 @@
     private IRegressionModel RegressLR(IRegressionProblemData problemData, string inputVar, double[] target) {
       // Umständlich!
-      var ds = new Dataset(new string[] { inputVar, "y" }, new[] { problemData.Dataset.GetDoubleValues(inputVar, problemData.TrainingIndices).ToList(), target.ToList() });
-      var pd = new RegressionProblemData(ds, new string[] { inputVar }, "y");
+      var ds = ((Dataset)problemData.Dataset).ToModifiable();
+      ds.ReplaceVariable(problemData.TargetVariable, target.Concat(Enumerable.Repeat(0.0, ds.Rows - target.Length)).ToList<double>());
+      var pd = new RegressionProblemData(ds, new string[] { inputVar }, problemData.TargetVariable);
+      pd.TrainingPartition.Start = problemData.TrainingPartition.Start;
+      pd.TrainingPartition.End = problemData.TrainingPartition.End;
+      pd.TestPartition.Start = problemData.TestPartition.Start;
+      pd.TestPartition.End = problemData.TestPartition.End;
       double rmsError, cvRmsError;
       return LinearRegression.CreateLinearRegressionSolution(pd, out rmsError, out cvRmsError).Model;
@@ -142 +149 @@
 
     private IRegressionModel RegressSpline(IRegressionProblemData problemData, string inputVar, double[] target, double lambda) {
-      // Umständlich!
-      return Splines.CalculateSmoothingSplineReinsch(
-        problemData.Dataset.GetDoubleValues(inputVar, problemData.TrainingIndices).ToArray(),
-        (double[])target.Clone(),
-        new string[] { inputVar },
-        s: lambda,
-        targetVar: problemData.TargetVariable);
-    }
-
+      if (problemData.Dataset.VariableHasType<double>(inputVar)) {
+        // Umständlich!
+        return Splines.CalculateSmoothingSplineReinschWithAutomaticTolerance(
+          problemData.Dataset.GetDoubleValues(inputVar, problemData.TrainingIndices).ToArray(),
+          (double[])target.Clone(),
+          new string[] { inputVar },
+          targetVar: problemData.TargetVariable);
+      } else return new ConstantModel(target.Average(), problemData.TargetVariable);
+    }
+    private IRegressionModel RegressRF(IRegressionProblemData problemData, string inputVar, double[] target, double lambda) {
+      if (problemData.Dataset.VariableHasType<double>(inputVar)) {
+        // Umständlich!
+        var ds = ((Dataset)problemData.Dataset).ToModifiable();
+        ds.ReplaceVariable(problemData.TargetVariable, target.Concat(Enumerable.Repeat(0.0, ds.Rows - target.Length)).ToList<double>());
+        var pd = new RegressionProblemData(ds, new string[] { inputVar }, problemData.TargetVariable);
+        pd.TrainingPartition.Start = problemData.TrainingPartition.Start;
+        pd.TrainingPartition.End = problemData.TrainingPartition.End;
+        pd.TestPartition.Start = problemData.TestPartition.Start;
+        pd.TestPartition.End = problemData.TestPartition.End;
+        double rmsError, oobRmsError;
+        double avgRelError, oobAvgRelError; 
+        return RandomForestRegression.CreateRandomForestRegressionModel(pd, 100, 0.5, 0.5, 1234, out rmsError, out avgRelError, out oobRmsError, out oobAvgRelError);
+      } else return new ConstantModel(target.Average(), problemData.TargetVariable);
+    }
   }
 

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

@@ -36 +36 @@
 using HeuristicLab.Problems.DataAnalysis.Symbolic;
 using HeuristicLab.Problems.DataAnalysis.Symbolic.Regression;
+using HeuristicLab.Random;
 using MathNet.Numerics.Interpolation;
 using MathNet.Numerics.LinearAlgebra.Double;
@@ -70 +71 @@
           "Smoothing Spline (Mine)",
           "Smoothing Spline (Reinsch)",
+          "Smoothing Spline (Reinsch with automatic tolerance determination using CV)",
+          "B-Spline Smoothing"
       }.Select(s => new StringValue(s)));
 
       Parameters.Add(new ConstrainedValueParameter<StringValue>("Type", "The type of spline (as supported by alglib)", validTypes, validTypes.First()));
       Parameters.Add(new ValueParameter<DoubleValue>("Lambda", "Regularization parameter for smoothing splines (0..+inf)", new DoubleValue(100)));
+      Parameters.Add(new ValueParameter<DoubleValue>("StdDev (noise)", "Known error in y values. Used only be Reinsch Smoothing Splines", new DoubleValue(0.1)));
       Problem = new RegressionProblem();
     }
@@ -145 +149 @@
               break;
             }
+          case "Smoothing Spline (Reinsch with automatic tolerance determination using CV)": {
+              var model = CalculateSmoothingSplineReinschWithAutomaticTolerance(x, y, inputVars, Problem.ProblemData.TargetVariable);
+              Results.Add(new Result("Solution", new RegressionSolution(model, (IRegressionProblemData)Problem.ProblemData.Clone())));
+              break;
+            }
+          case "B-Spline Smoothing": {
+              BSplineSmoothing(x, y, inputVars);
+              break;
+            }
           default: throw new NotSupportedException();
         }
@@ -151 +164 @@
     }
 
-    // array with non-zero lower index
-    private class MyArray<T> {
-      private T[] arr;
-      private int lowerBound;
-
-      public T this[int key] {
-        get {
-          return arr[key - lowerBound];
-        }
-        set {
-          arr[key - lowerBound] = value;
-        }
-      }
-
-      public MyArray(int lowerBound, int numElements) {
-        this.lowerBound = lowerBound;
-        arr = new T[numElements];
-      }
-      public MyArray(int lowerBound, T[] source) : this(lowerBound, source.Length) {
-        Array.Copy(source, arr, source.Length);
-      }
-
-      public T[] ToArray() {
-        var res = new T[arr.Length];
-        Array.Copy(arr, res, res.Length);
-        return res;
-      }
+    public void BSplineSmoothing(double[] x, double[] y, string[] inputVars) {
+      Array.Sort(x, y);
+
+      CubicSpline2(x, y, inputVars);
+
+      //  // see Elements of Statistical Learning Appendix: Computations for Splines
+      //  int M = 4; // order of spline (cubic)
+      //  int K = x.Length;
+      //  int n = x.Length;
+      // 
+      //  double[] t = new double[K + 2 * M];
+      //  for (int i = 0; i < M; i++) t[i] = x[0];
+      //  for (int i = K; i < K + 2 * M; i++) t[i] = x[K - 1];
+      // 
+      //  double[,] B = new double[n, n + 4]; // basis function matrix
+      //  double[,] W = new double[n + 4, n + 4]; // penalty matrix (W[j,k] = 
+      //  for (int i = 0; i < M; i++) B[i] = new double[K + 2 * M];
+      //  for (int j = 0; j < K; j++) {
+      //    for (int i = 0; i < K + 2M - 1; i++) {
+      //      if (t[i] <= x[j] && x[j] < t[i + 1])
+      //        B[1][i] = 1.0;
+      //    }
+      //  }
+
+
+    }
+
+    private struct SplineData {
+      public double a, b, c, d, x, y;
+    };
+
+
+    /*
+     * The procedure Quincunx, takes as arguments
+       the vectors u, v and w which are respectively the diagonal, the first supradiagonal
+       and the second supradiagonal of the banded matrix. The vector on
+       the LHS of the equation (80) is placed in q which contains the solution on the
+       completion of the procedure.
+    */
+    private void Quincunx(int n, MyArray<double> u, MyArray<double> v, MyArray<double> w, MyArray<double> q) {
+      // { factorisation}
+      u[-1] = 0;
+      u[0] = 0;
+      for (int j = 1; j <= n - 1; j++) {
+        u[j] = u[j] - u[j - 2] * Math.Pow(w[j - 2], 2) - u[j - 1] * Math.Pow(v[j - 1], 2);
+        v[j] = (v[j] - u[j - 1] * v[j - 1] * w[j - 1]) / u[j];
+        w[j] = w[j] / u[j];
+      }
+      // { forward substitution}
+      for (int j = 1; j <= n - 1; j++) {
+        q[j] = q[j] - v[j - 1] * q[j - 1] - w[j - 2] * q[j - 2];
+      }
+      for (int j = 1; j < n - 1; j++) {
+        q[j] = q[j] / u[j];
+      }
+      // { back substitution}
+      q[n + 1] = 0;
+      q[n] = 0;
+      for (int j = n - 1; j >= 1; j--) {
+        q[j] = q[j] - v[j] * q[j + 1] - w[j] * q[j + 2];
+      }
+    }
+
+    public void CubicSpline2(double[] x, double[] y, string[] inputVars) {
+      // see Pollock: Smoothing Splines
+      int n = x.Length;
+      double lambda = 1.0;
+      double[] sigma = Enumerable.Repeat(1.0, n).ToArray();
+      SplineData[] S = new SplineData[x.Length + 1];
+      for (int i = 0; i < x.Length; i++) {
+        S[i].x = x[i];
+        S[i].y = y[i];
+      }
+      S[x.Length].x = S[x.Length - 1].x;
+      S[x.Length].y = S[x.Length - 1].y;
+
+      // var
+      double[] h = new double[n],
+        r = new double[n],
+        f = new double[n],
+        p = new double[n];
+      var q = new MyArray<double>(-1, n + 3);
+      var u = new MyArray<double>(-1, n + 3);
+      var v = new MyArray<double>(-1, n + 3);
+      var w = new MyArray<double>(-1, n + 3);
+      double mu;
+
+      mu = 2 * (1 - lambda) / (3 * lambda);
+      h[0] = S[1].x - S[0].x;
+      r[0] = 3 / h[0];
+      for (int i = 1; i < n - 1; i++) {
+        h[i] = S[i + 1].x - S[i].x;
+        r[i] = 3 / h[i];
+        f[i] = -(r[i - 1] + r[i]);
+        p[i] = 2 * (S[i + 1].x - S[i - 1].x);
+        q[i] = 3 * (S[i + 1].y - S[i].y) / h[i] - 3 * (S[i].y - S[i - 1].y) / h[i - 1];
+      }
+      for (int i = 1; i < n; i++) {
+        u[i] = Math.Pow(r[i - 1], 2) * sigma[i - 1] + Math.Pow(f[i], 2) * sigma[i] + Math.Pow(r[i], 2) * sigma[i + 1];    // TODO Sigma 1..n
+        u[i] = mu * u[i] + p[i];
+        v[i] = f[i] * r[i] * sigma[i] + r[i] * f[i + 1] * sigma[i + 1];
+        v[i] = mu * v[i] + h[i];
+        w[i] = mu * r[i] * r[i + 1] * sigma[i + 1];
+      }
+      Quincunx(n, u, v, w, q);
+      // { Spline P arameters}
+      S[0].d = S[0].y - mu * r[0] * q[1] * sigma[0];
+      S[1].d = S[1].y - mu * (f[1] * q[1] + r[1] * q[2]) * sigma[0];
+      S[0].a = q[1] / (3 * h[0]);
+      S[0].b = 0;
+      S[0].c = (S[1].d - S[0].d) / h[0] - q[1] * h[0] / 3;
+      r[0] = 0;
+      for (int j = 1; j < n; j++) {
+        S[j].a = (q[j + 1] - q[j]) / (3 * h[j]);
+        S[j].b = q[j];
+        S[j].c = (q[j] + q[j - 1]) * h[j - 1] + S[j - 1].c;
+        S[j].d = r[j - 1] * q[j - 1] + f[j] * q[j] + r[j] * q[j + 1];
+        S[j].d = y[j] - mu * S[j].d * sigma[j];
+      }
+
+      //   { SmoothingSpline}
     }
 
     public void CalculateSmoothingSplineReinsch(double[] xOrig, double[] yOrig, string[] inputVars) {
-      double s = ((IValueParameter<DoubleValue>)Parameters["Lambda"]).Value.Value; // controls extent of smoothing
-
-      var model = CalculateSmoothingSplineReinsch(xOrig, yOrig, inputVars, s, Problem.ProblemData.TargetVariable);
+      double s = xOrig.Length + 1;
+      double stdTol = ((IValueParameter<DoubleValue>)Parameters["StdDev (noise)"]).Value.Value; // controls extent of smoothing
+
+      var model = CalculateSmoothingSplineReinsch(xOrig, yOrig, inputVars, stdTol, s, Problem.ProblemData.TargetVariable);
 
       Results.Add(new Result("Solution", new RegressionSolution(model, (IRegressionProblemData)Problem.ProblemData.Clone())));
@@ -190 +298 @@
     }
 
-    public static IRegressionModel CalculateSmoothingSplineReinsch(double[] xOrig, double[] yOrig, string[] inputVars, double s, string targetVar) {
+
+
+    public static IRegressionModel CalculateSmoothingSplineReinschWithAutomaticTolerance(double[] xOrig, double[] yOrig, string[] inputVars, string targetVar) {
+      int n = xOrig.Length;
+      var rows = Enumerable.Range(0, n);
+      var rand = new FastRandom(1234);
+      var trainRows = rows.Shuffle(rand).Take((int)(n * 0.66)).ToArray();
+      var testRows = rows.Except(trainRows).ToArray();
+
+      var trainX = trainRows.Select(i => xOrig[i]).ToArray();
+      var trainY = trainRows.Select(i => yOrig[i]).ToArray();
+      var testX = testRows.Select(i => xOrig[i]).ToArray();
+      var testY = testRows.Select(i => yOrig[i]).ToArray();
+
+      var shuffledDs = new Dataset(inputVars.Concat(new string[] { targetVar }), new[] { trainX.Concat(testX).ToArray(), trainY.Concat(testY).ToArray() });
+      var shuffeledProblemData = new RegressionProblemData(shuffledDs, inputVars, targetVar);
+      shuffeledProblemData.TrainingPartition.Start = 0;
+      shuffeledProblemData.TrainingPartition.End = trainX.Length;
+      shuffeledProblemData.TestPartition.Start = trainX.Length;
+      shuffeledProblemData.TestPartition.End = shuffledDs.Rows + 1;
+
+      double curTol = trainY.StandardDeviation();
+      double prevTestRMSE = double.PositiveInfinity;
+      double testRMSE = double.PositiveInfinity;
+      IRegressionModel prevModel = null;
+      IRegressionModel model = null;
+      do {                                     
+        prevModel = model;
+        prevTestRMSE = testRMSE;
+        model = CalculateSmoothingSplineReinsch(trainX, trainY, inputVars, curTol, s: trainX.Length + 1, targetVar: targetVar);
+
+        var sol = model.CreateRegressionSolution(shuffeledProblemData);
+        var trainRMSE = sol.TrainingRootMeanSquaredError;
+        testRMSE = sol.TestRootMeanSquaredError;
+        curTol = Math.Max(curTol * 0.5, 1e-12 * trainY.StandardDeviation());
+      } while (testRMSE < prevTestRMSE);
+      return prevModel;
+    }
+
+
+    public static IRegressionModel CalculateSmoothingSplineReinsch(double[] xOrig, double[] yOrig, string[] inputVars, double stdTol, double s, string targetVar) {
       var minX = xOrig.Min();
       var maxX = xOrig.Max();
       var range = maxX - minX;
 
-      SortAndBin(xOrig, yOrig, out xOrig, out yOrig, scaling: true);
+      double[] w = Enumerable.Repeat(stdTol, xOrig.Length).ToArray();
+      SortAndBin((double[])xOrig.Clone(), (double[])yOrig.Clone(), w, out xOrig, out yOrig, out w, scaling: false);
 
       // See Smoothing by Spline Functions, Reinsch, 1967
@@ -202 +351 @@
       var x = new MyArray<double>(1, xOrig);
       var y = new MyArray<double>(1, yOrig);
-      var inv_dy = new MyArray<double>(1, Enumerable.Repeat(1.0, n).ToArray()); // equal weights
+      var inv_dy = new MyArray<double>(1, w);
 
       int n1 = 1;
-      int n2 = n - 1;
+      int n2 = n;
 
       // results
@@ -283 +432 @@
       }
 
-      return new ReinschSmoothingSplineModel(a.ToArray(), b.ToArray(), c.ToArray(), d.ToArray(), xOrig, targetVar, inputVars, minX, 1.0 / range);
+      return new ReinschSmoothingSplineModel(a, b, c, d, x, targetVar, inputVars);
     }
 
     private void CalculateSmoothingSpline(double[] x, double[] y, string[] inputVars) {
       // see Smoothing and Non-Parametric Regression, Germán Rodríguez, 2001 2.3.1
-      SortAndBin(x, y, out x, out y, scaling: false);
+      double[] w = Enumerable.Repeat(1.0, x.Length).ToArray();  // weights necessary for sortAndBin but are ignored below (TODO)
+      SortAndBin(x, y, w, out x, out y, out w, scaling: false);
       int n = x.Length;
 
@@ -314 +464 @@
       }
 
+      // WInvD = W.Cholesky().Solve(delta);
       var solvResult = W.TrySolveIterative(delta, WInvD, new MlkBiCgStab());
 
@@ -339 +490 @@
 
       Results.Add(new Result("Solution", new RegressionSolution(new SmoothingSplineModel(s, x, Problem.ProblemData.TargetVariable, inputVars),
- (IRegressionProblemData)Problem.ProblemData.Clone())));
-    }
-
-    private static void SortAndBin(double[] x, double[] y, out double[] x2, out double[] y2, bool scaling = false) {
-      // sort y by x
-      Array.Sort(x, y);
+  (IRegressionProblemData)Problem.ProblemData.Clone())));
+    }
+
+    private static void SortAndBin(double[] x, double[] y, double[] w, out double[] x2, out double[] y2, out double[] w2, bool scaling = false) {
+      var sortedIdx = Enumerable.Range(0, x.Length).ToArray();
+      // sort by x
+      Array.Sort(x, sortedIdx);
 
       var xl = new List<double>();
       var yl = new List<double>();
+      var wl = new List<double>();
 
       int n = x.Length;
@@ -360 +513 @@
           double sumX = 0.0;
           double sumY = 0.0;
+          double sumW = 0.0;
           while (j < n && x[j] - x[i] <= binSize) {
             k++;
             sumX += x[j];
-            sumY += y[j];
+            sumY += y[sortedIdx[j]];
+            sumW += w[sortedIdx[j]];
             j++;
           }
           var avgX = sumX / k;
           if (scaling) avgX = (avgX - x[0]) / range;
-          xl.Add(avgX);     // scale 
+          xl.Add(avgX);
           yl.Add(sumY / k);
+          wl.Add(sumW);
           i += k;
         }
@@ -376 +532 @@
       x2 = xl.ToArray();
       y2 = yl.ToArray();
+      w2 = wl.ToArray();
     }
 
@@ -390 +547 @@
 
 
+  // array with non-zero lower index
+  internal class MyArray<T> {
+    private T[] arr;
+    private int lowerBound;
+
+    public int Length { get { return arr.Length; } }
+
+    public T this[int key] {
+      get {
+        return arr[key - lowerBound];
+      }
+      set {
+        arr[key - lowerBound] = value;
+      }
+    }
+
+    public MyArray(int lowerBound, int numElements) {
+      this.lowerBound = lowerBound;
+      arr = new T[numElements];
+    }
+    public MyArray(int lowerBound, T[] source) : this(lowerBound, source.Length) {
+      Array.Copy(source, arr, source.Length);
+    }
+
+    public T[] ToArray() {
+      var res = new T[arr.Length];
+      Array.Copy(arr, res, res.Length);
+      return res;
+    }
+  }
+
 
   // UNFINISHED
-  public class ReinschSmoothingSplineModel : NamedItem, IRegressionModel {
-    private double[] a;
-    private double[] b;
-    private double[] c;
-    private double[] d;
-    private double[] x;
+  internal class ReinschSmoothingSplineModel : NamedItem, IRegressionModel {
+    private MyArray<double> a;
+    private MyArray<double> b;
+    private MyArray<double> c;
+    private MyArray<double> d;
+    private MyArray<double> x;
     private double offset;
     private double scale;
@@ -418 +606 @@
       this.offset = orig.offset;
     }
-    public ReinschSmoothingSplineModel(double[] a, double[] b, double[] c, double[] d, double[] x, string targetVar, string[] inputs, double offset, double scale) : base("SplineModel", "SplineModel") {
+    public ReinschSmoothingSplineModel(
+      MyArray<double> a,
+      MyArray<double> b,
+      MyArray<double> c,
+      MyArray<double> d,
+      MyArray<double> x,
+      string targetVar, string[] inputs, double offset = 0, double scale = 1) : base("SplineModel", "SplineModel") {
       this.a = a;
       this.b = b;
@@ -428 +622 @@
       this.scale = scale;
       this.offset = offset;
+
+      // extrapolate for xx > x[n2]
+      b[b.Length] = b[b.Length - 1];
+      d[b.Length] = d[d.Length - 1];
     }
 
@@ -441 +639 @@
       int n = x.Length;
       foreach (var xx in dataset.GetDoubleValues(VariablesUsedForPrediction.First(), rows).Select(xi => (xi - offset) * scale)) {
-        if (xx < x[0]) yield return a[0];
-        else if (xx >= x[n - 1]) yield return a[n - 2];
-        else {
+        if (xx <= x[1]) {
+          double h = xx - x[1];
+          yield return a[1] + h * (b[1] + h * (c[1] + h * d[1]));
+        } else if (xx >= x[n]) {
+          double h = xx - x[n];
+          yield return a[n] + h * (b[n] + h * (c[n] + h * d[n]));
+        } else {
           // binary search
-          int lower = 0;
-          int upper = n - 1;
+          int lower = 1;
+          int upper = n;
           while (true) {
             if (upper < lower) throw new InvalidProgramException();