source: branches/MathNetNumerics-Exploration-2789/HeuristicLab.Algorithms.DataAnalysis.Experimental/GAM.cs @ 15450

Last change on this file since 15450 was 15450, checked in by gkronber, 4 years ago

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

File size: 15.0 KB
Line 
1#region License Information
2/* HeuristicLab
3 * Copyright (C) 2002-2016 Heuristic and Evolutionary Algorithms Laboratory (HEAL)
4 *
5 * This file is part of HeuristicLab.
6 *
7 * HeuristicLab is free software: you can redistribute it and/or modify
8 * it under the terms of the GNU General Public License as published by
9 * the Free Software Foundation, either version 3 of the License, or
10 * (at your option) any later version.
11 *
12 * HeuristicLab is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
15 * GNU General Public License for more details.
16 *
17 * You should have received a copy of the GNU General Public License
18 * along with HeuristicLab. If not, see <http://www.gnu.org/licenses/>.
19 */
20#endregion
21
22using System;
23using System.Collections.Concurrent;
24using System.Collections.Generic;
25using System.Linq;
26using System.Threading;
27using System.Threading.Tasks;
28using HeuristicLab.Analysis;
29using HeuristicLab.Common;
30using HeuristicLab.Core;
31using HeuristicLab.Data;
32using HeuristicLab.Encodings.SymbolicExpressionTreeEncoding;
33using HeuristicLab.Optimization;
34using HeuristicLab.Parameters;
35using HeuristicLab.Persistence.Default.CompositeSerializers.Storable;
36using HeuristicLab.Problems.DataAnalysis;
37using HeuristicLab.Problems.DataAnalysis.Symbolic;
38using HeuristicLab.Problems.DataAnalysis.Symbolic.Regression;
39
40namespace HeuristicLab.Algorithms.DataAnalysis.Experimental {
41  // UNFINISHED
42  [Item("Generalized Additive Modelling", "GAM")]
43  [Creatable(CreatableAttribute.Categories.DataAnalysisRegression, Priority = 102)]
44  [StorableClass]
45  public sealed class GAM : FixedDataAnalysisAlgorithm<IRegressionProblem> {
46
47    private const string LambdaParameterName = "Lambda";
48    private const string MaxIterationsParameterName = "Max iterations";
49    private const string MaxInteractionsParameterName = "Max interactions";
50
51    public IFixedValueParameter<DoubleValue> LambdaParameter {
52      get { return (IFixedValueParameter<DoubleValue>)Parameters[LambdaParameterName]; }
53    }
54    public IFixedValueParameter<IntValue> MaxIterationsParameter {
55      get { return (IFixedValueParameter<IntValue>)Parameters[MaxIterationsParameterName]; }
56    }
57    public IFixedValueParameter<IntValue> MaxInteractionsParameter {
58      get { return (IFixedValueParameter<IntValue>)Parameters[MaxInteractionsParameterName]; }
59    }
60
61    public double Lambda {
62      get { return LambdaParameter.Value.Value; }
63      set { LambdaParameter.Value.Value = value; }
64    }
65    public int MaxIterations {
66      get { return MaxIterationsParameter.Value.Value; }
67      set { MaxIterationsParameter.Value.Value = value; }
68    }
69    public int MaxInteractions {
70      get { return MaxInteractionsParameter.Value.Value; }
71      set { MaxInteractionsParameter.Value.Value = value; }
72    }
73
74    [StorableConstructor]
75    private GAM(bool deserializing) : base(deserializing) { }
76    [StorableHook(HookType.AfterDeserialization)]
77    private void AfterDeserialization() {
78    }
79
80    private GAM(GAM original, Cloner cloner)
81      : base(original, cloner) {
82    }
83    public override IDeepCloneable Clone(Cloner cloner) {
84      return new GAM(this, cloner);
85    }
86
87    public GAM()
88      : base() {
89      Problem = new RegressionProblem();
90      Parameters.Add(new FixedValueParameter<DoubleValue>(LambdaParameterName, "Regularization for smoothing splines", new DoubleValue(1.0)));
91      Parameters.Add(new FixedValueParameter<IntValue>(MaxIterationsParameterName, "", new IntValue(100)));
92      Parameters.Add(new FixedValueParameter<IntValue>(MaxInteractionsParameterName, "", new IntValue(1)));
93    }
94
95
96    protected override void Run(CancellationToken cancellationToken) {
97      double lambda = Lambda;
98      int maxIters = MaxIterations ;
99      int maxInteractions = MaxInteractions;
100      if (maxInteractions < 1 || maxInteractions > 5) throw new ArgumentException("Max interactions is outside the valid range [1 .. 5]");
101
102      // calculates a GAM model using a linear representation + independent non-linear functions of each variable
103      // using backfitting algorithm (see The Elements of Statistical Learning page 298)
104
105      var problemData = Problem.ProblemData;
106      var y = problemData.TargetVariableTrainingValues.ToArray();
107      var avgY = y.Average();
108      var inputVars = Problem.ProblemData.AllowedInputVariables.ToArray();
109      var nTerms = 0; // inputVars.Length; // LR
110      for (int i = 1; i <= maxInteractions; i++) {
111        nTerms += inputVars.Combinations(i).Count();
112      }
113
114      IRegressionModel[] f = new IRegressionModel[nTerms];
115      for (int i = 0; i < f.Length; i++) {
116        f[i] = new ConstantModel(0.0, problemData.TargetVariable);
117      }
118
119      var rmseTable = new DataTable("RMSE");
120      var rmseRow = new DataRow("RMSE (train)");
121      var rmseRowTest = new DataRow("RMSE (test)");
122      rmseTable.Rows.Add(rmseRow);
123      rmseTable.Rows.Add(rmseRowTest);
124
125      Results.Add(new Result("RMSE", rmseTable));
126      rmseRow.Values.Add(CalculateResiduals(problemData, f, -1, avgY, problemData.TrainingIndices).StandardDeviation()); // -1 index to use all predictors
127      rmseRowTest.Values.Add(CalculateResiduals(problemData, f, -1, avgY, problemData.TestIndices).StandardDeviation());
128
129      // for analytics
130      double[] rss = new double[f.Length];
131      string[] terms = new string[f.Length];
132      Results.Add(new Result("RSS Values", typeof(DoubleMatrix)));
133
134      var combinations = new List<string[]>();
135      for(int i=1;i<=maxInteractions;i++)
136        combinations.AddRange(HeuristicLab.Common.EnumerableExtensions.Combinations(inputVars, i).Select(c => c.ToArray()));
137      // combinations.Add(new string[] { "X1", "X2" });
138      // combinations.Add(new string[] { "X3", "X4" });
139      // combinations.Add(new string[] { "X5", "X6" });
140      // combinations.Add(new string[] { "X1", "X7", "X9" });
141      // combinations.Add(new string[] { "X3", "X6", "X10" });
142
143
144
145      // until convergence
146      int iters = 0;
147      var t = new double[y.Length];
148      while (iters++ < maxIters) {
149        int j = 0;
150        //foreach (var inputVar in inputVars) {
151        //  var res = CalculateResiduals(problemData, f, j, avgY, problemData.TrainingIndices);
152        //  rss[j] = res.Variance();
153        //  terms[j] = inputVar;
154        //  f[j] = RegressLR(problemData, inputVar, res);
155        //  j++;
156        //}
157
158
159
160        foreach (var element in combinations) {
161          var res = CalculateResiduals(problemData, f, j, avgY, problemData.TrainingIndices);
162          rss[j] = res.Variance();
163          terms[j] = string.Format("f({0})", string.Join(",", element));
164          f[j] = RegressSpline(problemData, element.ToArray(), res, lambda);
165          j++;
166        }
167
168        rmseRow.Values.Add(CalculateResiduals(problemData, f, -1, avgY, problemData.TrainingIndices).StandardDeviation()); // -1 index to use all predictors
169        rmseRowTest.Values.Add(CalculateResiduals(problemData, f, -1, avgY, problemData.TestIndices).StandardDeviation());
170
171        // calculate table with residual contributions of each term
172        var rssTable = new DoubleMatrix(rss.Length, 1, new string[] { "RSS" }, terms);
173        for (int i = 0; i < rss.Length; i++) rssTable[i, 0] = rss[i];
174        Results["RSS Values"].Value = rssTable;
175
176        if (cancellationToken.IsCancellationRequested) break;
177      }
178
179      var model = new RegressionEnsembleModel(f.Concat(new[] { new ConstantModel(avgY, problemData.TargetVariable) }));
180      model.AverageModelEstimates = false;
181      var solution = model.CreateRegressionSolution((IRegressionProblemData)problemData.Clone());
182      Results.Add(new Result("Ensemble solution", solution));
183    }
184
185    private double[] CalculateResiduals(IRegressionProblemData problemData, IRegressionModel[] f, int j, double avgY, IEnumerable<int> rows) {
186      var y = problemData.Dataset.GetDoubleValues(problemData.TargetVariable, rows);
187      double[] t = y.Select(yi => yi - avgY).ToArray();
188      // collect other predictions
189      for (int k = 0; k < f.Length; k++) {
190        if (k != j) {
191          var pred = f[k].GetEstimatedValues(problemData.Dataset, rows).ToArray();
192          // determine target for this smoother
193          for (int i = 0; i < t.Length; i++) {
194            t[i] -= pred[i];
195          }
196        }
197      }
198      return t;
199    }
200
201    private IRegressionModel RegressLR(IRegressionProblemData problemData, string inputVar, double[] target) {
202      // Umständlich!
203      var ds = ((Dataset)problemData.Dataset).ToModifiable();
204      ds.ReplaceVariable(problemData.TargetVariable, target.Concat(Enumerable.Repeat(0.0, ds.Rows - target.Length)).ToList<double>());
205      var pd = new RegressionProblemData(ds, new string[] { inputVar }, problemData.TargetVariable);
206      pd.TrainingPartition.Start = problemData.TrainingPartition.Start;
207      pd.TrainingPartition.End = problemData.TrainingPartition.End;
208      pd.TestPartition.Start = problemData.TestPartition.Start;
209      pd.TestPartition.End = problemData.TestPartition.End;
210      double rmsError, cvRmsError;
211      return LinearRegression.CreateLinearRegressionSolution(pd, out rmsError, out cvRmsError).Model;
212    }
213
214    // private IRegressionModel RegressSpline(IRegressionProblemData problemData, string inputVar, double[] target, double lambda) {
215    //   if (problemData.Dataset.VariableHasType<double>(inputVar)) {
216    //     // Umständlich!
217    //     return Splines.CalculatePenalizedRegressionSpline(
218    //       problemData.Dataset.GetDoubleValues(inputVar, problemData.TrainingIndices).ToArray(),
219    //       (double[])target.Clone(), lambda,
220    //       problemData.TargetVariable, new string[] { inputVar }
221    //       );
222    //   } else return new ConstantModel(target.Average(), problemData.TargetVariable);
223    // }
224    private IRegressionModel RegressSpline(IRegressionProblemData problemData, string[] inputVars, double[] target, double lambda) {
225      if (inputVars.All(problemData.Dataset.VariableHasType<double>)) {
226        var product = problemData.Dataset.GetDoubleValues(inputVars.First(), problemData.TrainingIndices).ToArray();
227        for (int i = 1; i < inputVars.Length; i++) {
228          product = product.Zip(problemData.Dataset.GetDoubleValues(inputVars[i], problemData.TrainingIndices), (pi, vi) => pi * vi).ToArray();
229        }
230        // CubicSplineGCV.CubGcvReport report;
231        // return CubicSplineGCV.CalculateCubicSpline(
232        //   product,
233        //   (double[])target.Clone(),
234        //   problemData.TargetVariable, inputVars, out report
235        //   );
236
237        double optTolerance; double cvRMSE;
238        // find tolerance
239        // var ensemble = Splines.CalculateSmoothingSplineReinsch(product, (double[])target.Clone(), inputVars, problemData.TargetVariable, out optTolerance, out cvRMSE);
240        // // train on whole data
241        // return Splines.CalculateSmoothingSplineReinsch(product, (double[])target.Clone(), inputVars, optTolerance, product.Length - 1, problemData.TargetVariable);
242
243
244        // find tolerance
245        var bestLambda = -5.0;
246        double bestCVRMSE = double.PositiveInfinity;
247        double avgTrainRMSE = double.PositiveInfinity;
248        for (double curLambda = -5.0; curLambda <= 6.0; curLambda += 1.0) {
249          var ensemble = Splines.CalculatePenalizedRegressionSpline(product, (double[])target.Clone(), curLambda, problemData.TargetVariable,  inputVars, out avgTrainRMSE, out cvRMSE);
250          Console.Write("{0} {1} {2}", curLambda, avgTrainRMSE, cvRMSE);
251          if (bestCVRMSE > cvRMSE) {
252            Console.Write(" *");
253            bestCVRMSE = cvRMSE;
254            bestLambda = curLambda;
255          }
256          Console.WriteLine();
257        }
258        // train on whole data
259       return Splines.CalculatePenalizedRegressionSpline(product, (double[])target.Clone(), bestLambda, problemData.TargetVariable, inputVars, out avgTrainRMSE, out cvRMSE);
260
261      } else return new ConstantModel(target.Average(), problemData.TargetVariable);
262    }
263
264    private IRegressionModel RegressRF(IRegressionProblemData problemData, string inputVar, double[] target, double lambda) {
265      if (problemData.Dataset.VariableHasType<double>(inputVar)) {
266        // Umständlich!
267        var ds = ((Dataset)problemData.Dataset).ToModifiable();
268        ds.ReplaceVariable(problemData.TargetVariable, target.Concat(Enumerable.Repeat(0.0, ds.Rows - target.Length)).ToList<double>());
269        var pd = new RegressionProblemData(ds, new string[] { inputVar }, problemData.TargetVariable);
270        pd.TrainingPartition.Start = problemData.TrainingPartition.Start;
271        pd.TrainingPartition.End = problemData.TrainingPartition.End;
272        pd.TestPartition.Start = problemData.TestPartition.Start;
273        pd.TestPartition.End = problemData.TestPartition.End;
274        double rmsError, oobRmsError;
275        double avgRelError, oobAvgRelError;
276        return RandomForestRegression.CreateRandomForestRegressionModel(pd, 100, 0.5, 0.5, 1234, out rmsError, out avgRelError, out oobRmsError, out oobAvgRelError);
277      } else return new ConstantModel(target.Average(), problemData.TargetVariable);
278    }
279  }
280
281
282  // UNFINISHED
283  public class RBFModel : NamedItem, IRegressionModel {
284    private alglib.rbfmodel model;
285
286    public string TargetVariable { get; set; }
287
288    public IEnumerable<string> VariablesUsedForPrediction { get; private set; }
289    private ITransformation<double>[] scaling;
290
291    public event EventHandler TargetVariableChanged;
292
293    public RBFModel(RBFModel orig, Cloner cloner) : base(orig, cloner) {
294      this.TargetVariable = orig.TargetVariable;
295      this.VariablesUsedForPrediction = orig.VariablesUsedForPrediction.ToArray();
296      this.model = (alglib.rbfmodel)orig.model.make_copy();
297      this.scaling = orig.scaling.Select(s => cloner.Clone(s)).ToArray();
298    }
299    public RBFModel(alglib.rbfmodel model, string targetVar, string[] inputs, IEnumerable<ITransformation<double>> scaling) : base("RBFModel", "RBFModel") {
300      this.model = model;
301      this.TargetVariable = targetVar;
302      this.VariablesUsedForPrediction = inputs;
303      this.scaling = scaling.ToArray();
304    }
305
306    public override IDeepCloneable Clone(Cloner cloner) {
307      return new RBFModel(this, cloner);
308    }
309
310    public IRegressionSolution CreateRegressionSolution(IRegressionProblemData problemData) {
311      return new RegressionSolution(this, (IRegressionProblemData)problemData.Clone());
312    }
313
314    public IEnumerable<double> GetEstimatedValues(IDataset dataset, IEnumerable<int> rows) {
315      double[] x = new double[VariablesUsedForPrediction.Count()];
316      double[] y;
317      foreach (var r in rows) {
318        int c = 0;
319        foreach (var v in VariablesUsedForPrediction) {
320          x[c] = scaling[c].Apply(dataset.GetDoubleValue(v, r).ToEnumerable()).First(); // OUCH!
321          c++;
322        }
323        alglib.rbfcalc(model, x, out y);
324        yield return y[0];
325      }
326    }
327  }
328}
Note: See TracBrowser for help on using the repository browser.