source: branches/2898_GeneralizedAdditiveModels/HeuristicLab.Algorithms.DataAnalysis/3.4/GAM/GeneralizedAdditiveModelAlgorithm.cs @ 15775

Last change on this file since 15775 was 15775, checked in by gkronber, 2 years ago

#2898 added simple implementation of GAM based on univariate penalized regression splines with the same penalization factor for each term

File size: 10.4 KB
Line 
1#region License Information
2/* HeuristicLab
3 * Copyright (C) 2002-2018 Heuristic and Evolutionary Algorithms Laboratory (HEAL)
4 * and the BEACON Center for the Study of Evolution in Action.
5 *
6 * This file is part of HeuristicLab.
7 *
8 * HeuristicLab is free software: you can redistribute it and/or modify
9 * it under the terms of the GNU General Public License as published by
10 * the Free Software Foundation, either version 3 of the License, or
11 * (at your option) any later version.
12 *
13 * HeuristicLab is distributed in the hope that it will be useful,
14 * but WITHOUT ANY WARRANTY; without even the implied warranty of
15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
16 * GNU General Public License for more details.
17 *
18 * You should have received a copy of the GNU General Public License
19 * along with HeuristicLab. If not, see <http://www.gnu.org/licenses/>.
20 */
21#endregion
22
23using System;
24using System.Collections.Generic;
25using System.Linq;
26using System.Threading;
27using HeuristicLab.Analysis;
28using HeuristicLab.Common;
29using HeuristicLab.Core;
30using HeuristicLab.Data;
31using HeuristicLab.Optimization;
32using HeuristicLab.Parameters;
33using HeuristicLab.Persistence.Default.CompositeSerializers.Storable;
34using HeuristicLab.Problems.DataAnalysis;
35using HeuristicLab.Random;
36
37namespace HeuristicLab.Algorithms.DataAnalysis {
38  [Item("Generalized Additive Model (GAM)",
39    "Generalized Additive Model Algorithm")]
40  [StorableClass]
41  [Creatable(CreatableAttribute.Categories.DataAnalysisRegression, Priority = 600)]
42  public sealed class GeneralizedAdditiveModelAlgorithm : FixedDataAnalysisAlgorithm<IRegressionProblem> {
43    #region ParameterNames
44
45    private const string IterationsParameterName = "Iterations";
46    private const string LambdaParameterName = "Lambda";
47    private const string SeedParameterName = "Seed";
48    private const string SetSeedRandomlyParameterName = "SetSeedRandomly";
49    private const string CreateSolutionParameterName = "CreateSolution";
50    #endregion
51
52    #region ParameterProperties
53
54    public IFixedValueParameter<IntValue> IterationsParameter {
55      get { return (IFixedValueParameter<IntValue>)Parameters[IterationsParameterName]; }
56    }
57
58    public IFixedValueParameter<DoubleValue> LambdaParameter {
59      get { return (IFixedValueParameter<DoubleValue>)Parameters[LambdaParameterName]; }
60    }
61
62    public IFixedValueParameter<IntValue> SeedParameter {
63      get { return (IFixedValueParameter<IntValue>)Parameters[SeedParameterName]; }
64    }
65
66    public FixedValueParameter<BoolValue> SetSeedRandomlyParameter {
67      get { return (FixedValueParameter<BoolValue>)Parameters[SetSeedRandomlyParameterName]; }
68    }
69
70    public IFixedValueParameter<BoolValue> CreateSolutionParameter {
71      get { return (IFixedValueParameter<BoolValue>)Parameters[CreateSolutionParameterName]; }
72    }
73
74    #endregion
75
76    #region Properties
77
78    public int Iterations {
79      get { return IterationsParameter.Value.Value; }
80      set { IterationsParameter.Value.Value = value; }
81    }
82
83    public double Lambda {
84      get { return LambdaParameter.Value.Value; }
85      set { LambdaParameter.Value.Value = value; }
86    }
87
88    public int Seed {
89      get { return SeedParameter.Value.Value; }
90      set { SeedParameter.Value.Value = value; }
91    }
92
93    public bool SetSeedRandomly {
94      get { return SetSeedRandomlyParameter.Value.Value; }
95      set { SetSeedRandomlyParameter.Value.Value = value; }
96    }
97
98    public bool CreateSolution {
99      get { return CreateSolutionParameter.Value.Value; }
100      set { CreateSolutionParameter.Value.Value = value; }
101    }
102
103    #endregion
104
105    [StorableConstructor]
106    private GeneralizedAdditiveModelAlgorithm(bool deserializing)
107      : base(deserializing) {
108    }
109
110    private GeneralizedAdditiveModelAlgorithm(GeneralizedAdditiveModelAlgorithm original, Cloner cloner)
111      : base(original, cloner) {
112    }
113
114    public override IDeepCloneable Clone(Cloner cloner) {
115      return new GeneralizedAdditiveModelAlgorithm(this, cloner);
116    }
117
118    public GeneralizedAdditiveModelAlgorithm() {
119      Problem = new RegressionProblem(); // default problem
120
121      Parameters.Add(new FixedValueParameter<IntValue>(IterationsParameterName,
122        "Number of iterations. Try a large value and check convergence of the error over iterations. Usually, only a few iterations (e.g. 10) are needed for convergence.", new IntValue(10)));
123      Parameters.Add(new FixedValueParameter<DoubleValue>(LambdaParameterName,
124        "The penalty parameter for the penalized regression splines. Set to a value between -8 (weak smoothing) and 8 (strong smooting). Usually, a value between -4 and 4 should be fine", new DoubleValue(3)));
125      Parameters.Add(new FixedValueParameter<IntValue>(SeedParameterName,
126        "The random seed used to initialize the new pseudo random number generator.", new IntValue(0)));
127      Parameters.Add(new FixedValueParameter<BoolValue>(SetSeedRandomlyParameterName,
128        "True if the random seed should be set to a random value, otherwise false.", new BoolValue(true)));
129      Parameters.Add(new FixedValueParameter<BoolValue>(CreateSolutionParameterName,
130        "Flag that indicates if a solution should be produced at the end of the run", new BoolValue(true)));
131      Parameters[CreateSolutionParameterName].Hidden = true;
132    }
133
134    protected override void Run(CancellationToken cancellationToken) {
135      // Set up the algorithm
136      if (SetSeedRandomly) Seed = new System.Random().Next();
137      var rand = new MersenneTwister((uint)Seed);
138
139      // calculates a GAM model using univariate non-linear functions
140      // using backfitting algorithm (see The Elements of Statistical Learning page 298)
141
142      // init
143      var problemData = Problem.ProblemData;
144      var ds = problemData.Dataset;
145      var trainRows = problemData.TrainingIndices;
146      var testRows = problemData.TestIndices;
147      var avgY = problemData.TargetVariableTrainingValues.Average();
148      var inputVars = problemData.AllowedInputVariables.ToArray();
149
150      int nTerms = inputVars.Length;
151
152      #region init results
153      // Set up the results display
154      var iterations = new IntValue(0);
155      Results.Add(new Result("Iterations", iterations));
156
157      var table = new DataTable("Qualities");
158      var rmseRow = new DataRow("RMSE (train)");
159      var rmseRowTest = new DataRow("RMSE (test)");
160      table.Rows.Add(rmseRow);
161      table.Rows.Add(rmseRowTest);
162      Results.Add(new Result("Qualities", table));
163      var curRMSE = new DoubleValue();
164      var curRMSETest = new DoubleValue();
165      Results.Add(new Result("RMSE (train)", curRMSE));
166      Results.Add(new Result("RMSE (test)", curRMSETest));
167
168      // calculate table with residual contributions of each term
169      var rssTable = new DoubleMatrix(nTerms, 1, new string[] { "RSS" }, inputVars);
170      Results.Add(new Result("RSS Values", rssTable));
171      #endregion
172
173      // start with a set of constant models = 0
174      IRegressionModel[] f = new IRegressionModel[nTerms];
175      for (int i = 0; i < f.Length; i++) {
176        f[i] = new ConstantModel(0.0, problemData.TargetVariable);
177      }
178      // init res which contains the current residual vector
179      double[] res = problemData.TargetVariableTrainingValues.Select(yi => yi - avgY).ToArray();
180      double[] resTest = problemData.TargetVariableTestValues.Select(yi => yi - avgY).ToArray();
181
182      curRMSE.Value = res.StandardDeviation();
183      curRMSETest.Value = resTest.StandardDeviation();
184      rmseRow.Values.Add(res.StandardDeviation()); 
185      rmseRowTest.Values.Add(resTest.StandardDeviation());
186
187
188      double lambda = Lambda;
189      var idx = Enumerable.Range(0, nTerms).ToArray();
190
191      // Loop until iteration limit reached or canceled.
192      for (int i = 0; i < Iterations && !cancellationToken.IsCancellationRequested; i++) {
193        // shuffle order of terms in each iteration to remove bias on earlier terms
194        idx.ShuffleInPlace(rand);
195        foreach (var inputIdx in idx) {
196          var inputVar = inputVars[inputIdx];
197          // first remove the effect of the previous model for the inputIdx (by adding the output of the current model to the residual)
198          AddInPlace(res, f[inputIdx].GetEstimatedValues(ds, trainRows));
199          AddInPlace(resTest, f[inputIdx].GetEstimatedValues(ds, testRows));
200
201          rssTable[inputIdx, 0] = res.Variance();
202          f[inputIdx] = RegressSpline(problemData, inputVar, res, lambda);
203
204          SubtractInPlace(res, f[inputIdx].GetEstimatedValues(ds, trainRows));
205          SubtractInPlace(resTest, f[inputIdx].GetEstimatedValues(ds, testRows));
206        }
207
208        curRMSE.Value = res.StandardDeviation(); 
209        curRMSETest.Value = resTest.StandardDeviation();
210        rmseRow.Values.Add(curRMSE.Value);
211        rmseRowTest.Values.Add(curRMSETest.Value);
212        iterations.Value = i;
213      }
214
215      // produce solution
216      if (CreateSolution) {
217        var model = new RegressionEnsembleModel(f.Concat(new[] { new ConstantModel(avgY, problemData.TargetVariable) }));
218        model.AverageModelEstimates = false;
219        var solution = model.CreateRegressionSolution((IRegressionProblemData)problemData.Clone());
220        Results.Add(new Result("Ensemble solution", solution));
221      }
222    }
223
224    private IRegressionModel RegressSpline(IRegressionProblemData problemData, string inputVar, double[] target, double lambda) {
225      var x = problemData.Dataset.GetDoubleValues(inputVar, problemData.TrainingIndices).ToArray();
226      var y = (double[])target.Clone();
227      int info;
228      alglib.spline1dinterpolant s;
229      alglib.spline1dfitreport rep;
230      int numKnots = (int)Math.Min(50, 3 * Math.Sqrt(x.Length)); // heuristic for number of knots  (forgot the source, but it is probably the R documentation or Elements of Statistical Learning)
231
232      alglib.spline1dfitpenalized(x, y, numKnots, lambda, out info, out s, out rep);
233
234      return new Spline1dModel(s.innerobj, problemData.TargetVariable, inputVar);
235    }
236
237
238    private static void AddInPlace(double[] a, IEnumerable<double> enumerable) {
239      int i = 0;
240      foreach (var elem in enumerable) {
241        a[i] += elem;
242        i++;
243      }
244    }
245
246    private static void SubtractInPlace(double[] a, IEnumerable<double> enumerable) {
247      int i = 0;
248      foreach (var elem in enumerable) {
249        a[i] -= elem;
250        i++;
251      }
252    }
253  }
254}
Note: See TracBrowser for help on using the repository browser.