source: branches/M5Regression/HeuristicLab.Algorithms.DataAnalysis/3.4/M5Regression/LeafTypes/M5Leaf.cs @ 15967

Last change on this file since 15967 was 15967, checked in by bwerth, 12 months ago

#2847 added logistic dampening and some minor changes

File size: 9.3 KB
Line 
1#region License Information
2/* HeuristicLab
3 * Copyright (C) 2002-2017 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.Generic;
24using System.Linq;
25using System.Threading;
26using HeuristicLab.Common;
27using HeuristicLab.Core;
28using HeuristicLab.Persistence.Default.CompositeSerializers.Storable;
29using HeuristicLab.Problems.DataAnalysis;
30
31namespace HeuristicLab.Algorithms.DataAnalysis {
32  [StorableClass]
33  [Item("M5Leaf", "A leaf type that uses linear models as leaf models. This is the standard for M5' regression")]
34  public class M5Leaf : LeafBase {
35    #region Constructors & Cloning
36    [StorableConstructor]
37    private M5Leaf(bool deserializing) : base(deserializing) { }
38    private M5Leaf(M5Leaf original, Cloner cloner) : base(original, cloner) { }
39    public M5Leaf() { }
40    public override IDeepCloneable Clone(Cloner cloner) {
41      return new M5Leaf(this, cloner);
42    }
43    #endregion
44
45    #region IModelType
46    public override bool ProvidesConfidence {
47      get { return false; }
48    }
49    public override IRegressionModel Build(IRegressionProblemData pd, IRandom random, CancellationToken cancellationToken, out int noParameters) {
50      if (pd.Dataset.Rows == 0) throw new ArgumentException("The number of training instances is too small to create an M5 leaf model");
51
52      if (pd.Dataset.Rows == 1)
53        return new ConstantLeaf().Build(pd, random, cancellationToken, out noParameters);
54
55      var means = pd.AllowedInputVariables.ToDictionary(n => n, n => pd.Dataset.GetDoubleValues(n, pd.TrainingIndices).Average());
56      var variances = pd.AllowedInputVariables.ToDictionary(n => n, n => pd.Dataset.GetDoubleValues(n, pd.TrainingIndices).Variance());
57      var used = pd.AllowedInputVariables.Where(v => !variances[v].IsAlmost(0.0)).ToList();
58
59      var classMean = pd.TargetVariableTrainingValues.Average();
60      var classVar = pd.TargetVariableTrainingValues.Variance();
61
62      var model = FindBestModel(variances, means, classMean, classVar, pd, used);
63      noParameters = 1 + model.Coefficients.Count;
64      return model;
65    }
66    public override int MinLeafSize(IRegressionProblemData pd) {
67      return 1;
68    }
69    #endregion
70
71    private static PreconstructedLinearModel FindBestModel(Dictionary<string, double> variances, Dictionary<string, double> means, double cMean, double cVar, IRegressionProblemData pd, IList<string> variables) {
72      Dictionary<string, double> coeffs;
73      double intercept;
74      do {
75        coeffs = DoRegression(pd, variables, variances, means, cMean, 1.0e-8, out intercept);
76        variables = DeselectColinear(variances, coeffs, cVar, pd, variables);
77      }
78      while (coeffs.Count != variables.Count);
79      var numAtts = variables.Count;
80      var numInst = pd.TrainingIndices.Count();
81      var fullMse = CalculateSE(coeffs, intercept, pd, variables);
82      var akaike = 1.0 * (numInst - numAtts) + 2 * numAtts;
83
84      var improved = true;
85      var currentNumAttributes = numAtts;
86
87      while (improved && currentNumAttributes > 1) {
88        improved = false;
89        currentNumAttributes--;
90        // Find attribute with smallest SC
91        var candidate = variables.ToDictionary(v => v, v => Math.Abs(coeffs[v] * Math.Sqrt(variances[v] / cVar)))
92          .OrderBy(x => x.Value).Select(x => x.Key).First();
93
94        var currVariables = variables.Where(v => !v.Equals(candidate)).ToList();
95        var currentIntercept = 0.0;
96        var currentCoeffs = DoRegression(pd, currVariables, variances, means, cMean, 1.0e-8, out currentIntercept);
97        var currentMse = CalculateSE(currentCoeffs, currentIntercept, pd, currVariables);
98        var currentAkaike = currentMse / fullMse * (numInst - numAtts) + 2 * currentNumAttributes;
99
100        if (!(currentAkaike < akaike)) continue;
101        improved = true;
102        akaike = currentAkaike;
103        coeffs = currentCoeffs;
104        intercept = currentIntercept;
105        variables = currVariables;
106      }
107
108      var pd2 = new RegressionProblemData(pd.Dataset, variables, pd.TargetVariable);
109      pd2.TestPartition.End = pd.TestPartition.End;
110      pd2.TestPartition.Start = pd.TestPartition.Start;
111      pd2.TrainingPartition.End = pd.TrainingPartition.End;
112      pd2.TrainingPartition.Start = pd.TrainingPartition.Start;
113
114      return new PreconstructedLinearModel(coeffs, intercept, pd.TargetVariable);
115    }
116
117    private static Dictionary<string, double> DoRegression(IRegressionProblemData pd, IList<string> variables, Dictionary<string, double> variances, Dictionary<string, double> means, double cmean, double ridge, out double intercept) {
118      //if (pd.TrainingIndices.Count() > variables.Count) {
119      //  var pd2 = new RegressionProblemData(pd.Dataset, variables, pd.TargetVariable);
120      //  pd2.TestPartition.End = pd.TestPartition.End;
121      //  pd2.TestPartition.Start = pd.TestPartition.Start;
122      //  pd2.TrainingPartition.End = pd.TrainingPartition.End;
123      //  pd2.TrainingPartition.Start = pd.TrainingPartition.Start;
124      //
125      //  double x1, x2;
126      //  var lm = PreconstructedLinearModel.CreateLinearModel(pd2, out x1, out x2);
127      //  intercept = lm.Intercept;
128      //  return lm.Coefficients;
129
130      var n = variables.Count;
131      var m = pd.TrainingIndices.Count();
132
133      var inTr = new double[n, m];
134      for (var i = 0; i < n; i++) {
135        var v = variables[i];
136        var vdata = pd.Dataset.GetDoubleValues(v, pd.TrainingIndices).ToArray();
137        var sd = Math.Sqrt(variances[v]);
138        var mean = means[v];
139        for (var j = 0; j < m; j++) {
140          inTr[i, j] = (vdata[j] - mean) / sd;
141        }
142      }
143
144      var y = new double[m, 1];
145      var ydata = pd.TargetVariableTrainingValues.ToArray();
146      for (var i = 0; i < m; i++)
147        y[i, 0] = ydata[i]; //no scaling for targets;
148
149
150      var aTy = new double[n, 1];
151      alglib.rmatrixgemm(n, 1, m, 1, inTr, 0, 0, 0, y, 0, 0, 0, 0, ref aTy, 0, 0); //aTy = inTr * y;
152      var aTa = new double[n, n];
153      alglib.rmatrixgemm(n, n, m, 1, inTr, 0, 0, 0, inTr, 0, 0, 1, 0, ref aTa, 0, 0); //aTa = inTr * t(inTr) +aTa //
154
155      var aTaDecomp = new double[n, n];
156      bool success;
157      var tries = 0;
158      double[] coefficients = null;
159      do {
160        for (var i = 0; i < n; i++) aTa[i, i] += ridge; // add ridge to diagonal to enforce singularity
161        try {
162          //solve "aTa * coefficients = aTy" for coefficients;
163          Array.Copy(aTa, 0, aTaDecomp, 0, aTa.Length);
164          alglib.spdmatrixcholesky(ref aTaDecomp, n, true);
165          int info;
166          alglib.densesolverreport report;
167          alglib.spdmatrixcholeskysolve(aTaDecomp, n, true, ydata, out info, out report, out coefficients);
168
169          if (info != 1) throw new Exception();
170          success = true;
171        }
172        catch (Exception) {
173          for (var i = 0; i < n; i++) aTa[i, i] -= ridge;
174          ridge *= 10; // increase ridge;
175          success = false;
176        }
177        finally {
178          tries++;
179        }
180      }
181      while (!success && tries < 100);
182      if (coefficients == null) throw new ArgumentException("No linear model could be built");
183
184      intercept = cmean;
185      var res = new Dictionary<string, double>();
186      for (var i = 0; i < n; i++) {
187        var v = variables[i];
188        res.Add(v, coefficients[i] /= Math.Sqrt(variances[v]));
189        intercept -= coefficients[i] * means[v];
190      }
191
192      return res;
193    }
194
195    private static IList<string> DeselectColinear(Dictionary<string, double> variances, Dictionary<string, double> coeffs, double cVar, IRegressionProblemData pd, IList<string> variables) {
196      var candidates = variables.ToDictionary(v => v, v => Math.Abs(coeffs[v] * Math.Sqrt(variances[v] / cVar))).Where(x => x.Value > 1.5).OrderBy(x => -x.Value).ToList();
197      if (candidates.Count == 0) return variables;
198      var c = candidates.First().Key;
199      return variables.Where(v => !v.Equals(c)).ToList();
200    }
201
202    private static double CalculateSE(Dictionary<string, double> coefficients, double intercept, IRegressionProblemData pd, IList<string> variables) {
203      return pd.TrainingIndices.Select(i => RegressionPrediction(i, pd, variables, coefficients, intercept) - pd.Dataset.GetDoubleValue(pd.TargetVariable, i)).Select(error => error * error).Sum();
204    }
205    private static double RegressionPrediction(int i, IRegressionProblemData pd, IList<string> variables, Dictionary<string, double> coefficients, double intercept) {
206      return intercept + variables.Sum(v => pd.Dataset.GetDoubleValue(v, i) * coefficients[v]);
207    }
208  }
209}
Note: See TracBrowser for help on using the repository browser.