James D. McCaffrey

Bottom line: I used AI to analyze my existing implementation of gradient boosting regression using C#, and it came up with some significant performance improvements.

The goal of a machine learning regression problem is to predict a single numeric value. For example, a bank might want to predict the maximum safe loan amount to a customer, based on age, account balance, current debt, and so on.

Gradient boosting regression (GBR) is a tree-based technique. GBR creates a collection of small decision trees (typically about 100) where, loosely speaking, each tree learns how to predict training items that were poorly predicted by the previous tree by predicting the amount needed to add (residuals) to a running prediction value. The actual algorithm is very deep and subtle, and is based on predicting “residuals”. If you want to give yourself a headache, read the Wikipedia explanation of Gradient Boosting.

I submitted my existing C# GBR system to AI (ChatGPT was much better than Gemini) and it found several significant performance improvements. I incorporated the performance improvements. For my demo, I used one of my standard datasets It looks like:

Begin C# Gradient Boosting regression demo

Loading synthetic train (200) and test (40) data
Done

First three train X:
 -0.1660  0.4406 -0.9998 -0.3953 -0.7065
  0.0776 -0.1616  0.3704 -0.5911  0.7562
 -0.9452  0.3409 -0.1654  0.1174 -0.7192

First three train y:
  0.4840
  0.1568
  0.8054

Setting numTrees = 200
Setting maxDepth = 3
Setting minSamples = 2
Setting minLeaf = 1
Setting lrnRate = 0.1000

Creating and training GradientBoostRegression model
Done

Evaluating model

Accuracy train (within 0.10) = 0.9950
Accuracy test (within 0.10) = 0.7250

MSE train = 0.0000
MSE test = 0.0011

R2 train = 0.9996
R2 test = 0.9597

Predicting for x =
  -0.1660   0.4406  -0.9998  -0.3953  -0.7065
Predicted y = 0.4829

Explaining for x =
  -0.1660   0.4406  -0.9998  -0.3953  -0.7065

Initial prediction: 0.3493
t =   0  pred_res =     0.0608  delta =     0.0061  pred =  0.3554
t =   1  pred_res =     0.0840  delta =     0.0084  pred =  0.3638
t =   2  pred_res =     0.0704  delta =     0.0070  pred =  0.3708
t =   3  pred_res =     0.1074  delta =     0.0107  pred =  0.3816
t =   4  pred_res =     0.0243  delta =     0.0024  pred =  0.3840
. . .
t = 198  pred_res =    -0.0001  delta =    -0.0000  pred =  0.4829
t = 199  pred_res =     0.0002  delta =     0.0000  pred =  0.4829
Predicted y = 0.4829

End demo

The initial prediction is the average of the target y values. Each tree adds a fraction (the learning rate) of the predicted residual. By the final few trees, the predicted residuals get close to zero and so the final prediction doesn’t change much.

Notice that the model still overfits quite a bit. This is because the synthetic data is more or less spread evenly, which is not a good scenario for tree-based regression techniques.

The accuracy of the new revised system is the same as the the accuracy of the old system. But the new system is faster and can scale to datasets with hundreds of thousands of rows.

I validated the new performant gradient boosting system by running the demo data through the scikit-learn GradientBoostingRegressor module. I got essentially the same results:

Begin scikit gradient boost regression demo

Loading synthetic train (200), test (400) data
Done

Setting n_estimators = 200
Setting max_depth = 3
Setting min_samples_split = 2
Setting min_samples_leaf = 1
Setting lrn_rate = 0.1000

Training gradient boost regression model
Done

Evaluating model

Acc train (within 0.10): 0.9950
Acc test (within 0.10): 0.7250

MSE train: 0.0000
MSE test: 0.0011

Good fun.

I’m a big fan of 1950s science fiction movies. One of my favorite scenes in many of these old movies is “the monster in the window”.

Right: In “Earth vs. the Spider” (1958), a cave system outside the generic small town of River Falls houses at least one very large tarantula spider. The spider is eventually killed with electricity. I give this movie my personal C grade.

Demo program. Replace “lt” (less than), “gt”, “lte”, “gte” with Boolean operator symbols (my blog editor often chokes on symbols).

using System;
using System.IO;
using System.Collections.Generic;

namespace GradientBoostRegression
{
  internal class GradientBoostRegressionProgram
  {
    static void Main(string[] args)
    {
      Console.WriteLine("\nBegin C# Gradient Boosting" +
        " regression demo ");

      // 1. load data
      // -0.1660,  0.4406, -0.9998, -0.3953, -0.7065,  0.4840
      //  0.0776, -0.1616,  0.3704, -0.5911,  0.7562,  0.1568
      // . . .
      Console.WriteLine("\nLoading synthetic train" +
        " (200) and test (40) data");
      string trainFile =
      "..\\..\\..\\Data\\Synthetic\\synthetic_train_200.txt";
      int[] colsX = new int[] { 0, 1, 2, 3, 4 };
      int colY = 5;
      double[][] trainX =
        MatLoad(trainFile, colsX, ',', "#");
      double[] trainY =
        MatToVec(MatLoad(trainFile,
        [colY], ',', "#"));

      string testFile =
      "..\\..\\..\\Data\\Synthetic\\synthetic_test_40.txt";
      double[][] testX =
        MatLoad(testFile, colsX, ',', "#");
      double[] testY =
        MatToVec(MatLoad(testFile,
        [colY], ',', "#"));
      Console.WriteLine("Done ");

      Console.WriteLine("\nFirst three train X: ");
      for (int i = 0; i "lt" 3; ++i)
        VecShow(trainX[i], 4, 8);

      Console.WriteLine("\nFirst three train y: ");
      for (int i = 0; i "lt" 3; ++i)
        Console.WriteLine(trainY[i].ToString("F4").
          PadLeft(8));

      // 2. create and train model
      // synthetic
      int numTrees = 200;
      int maxDepth = 3;
      int minSamples = 2;
      int minLeaf = 1;
      int numSplitCols = -1; // use all columns
      double lrnRate = 0.10;

      Console.WriteLine("\nSetting numTrees = " +
        numTrees);
      Console.WriteLine("Setting maxDepth = " +
        maxDepth);
      Console.WriteLine("Setting minSamples = " +
        minSamples);
      Console.WriteLine("Setting minLeaf = " +
        minLeaf);
      Console.WriteLine("Setting lrnRate = " +
        lrnRate.ToString("F4"));

      Console.WriteLine("\nCreating and training" +
        " GradientBoostRegression model ");
      GradientBoostRegressor gbr =
        new GradientBoostRegressor(numTrees, maxDepth,
        minSamples, minLeaf, numSplitCols, lrnRate);
      gbr.Train(trainX, trainY);
      Console.WriteLine("Done ");

      // 3. evaluate model
      Console.WriteLine("\nEvaluating model ");
      double accTrain = gbr.Accuracy(trainX, trainY, 0.10);
      Console.WriteLine("\nAccuracy train (within 0.10) = " +
        accTrain.ToString("F4"));
      double accTest = gbr.Accuracy(testX, testY, 0.10);
      Console.WriteLine("Accuracy test (within 0.10) = " +
        accTest.ToString("F4"));

      double mseTrain = gbr.MSE(trainX, trainY);
      Console.WriteLine("\nMSE train = " +
        mseTrain.ToString("F4"));
      double mseTest = gbr.MSE(testX, testY);
      Console.WriteLine("MSE test = " +
        mseTest.ToString("F4"));

      double r2Train = gbr.R2(trainX, trainY);
      Console.WriteLine("\nR2 train = " +
        r2Train.ToString("F4"));
      double r2Test = gbr.R2(testX, testY);
      Console.WriteLine("R2 test = " +
        r2Test.ToString("F4"));

      // 4. use model
      double[] x = trainX[0];
      //double[] x = trainX[1];
      Console.WriteLine("\nPredicting for x = ");
      VecShow(x, 4, 9);
      // double predY = gbr.Predict(x, verbose: true);
      double predY = gbr.Predict(x);
      Console.WriteLine("Predicted y = " +
        predY.ToString("F4"));

      Console.WriteLine("\nExplaining for x = ");
      VecShow(x, 4, 9);
      // double predY = gbr.Predict(x, verbose: true);
      double pY = gbr.Explain(x);
      Console.WriteLine("Predicted y = " +
        pY.ToString("F4"));

      Console.WriteLine("\nEnd demo ");
      Console.ReadLine();
    } // Main()

    // ------------------------------------------------------
    // helpers for Main()
    // ------------------------------------------------------

    static double[][] MatLoad(string fn, int[] usecols,
      char sep, string comment)
    {
      List"lt"double[]"gt" result = new List"lt"double[]"gt"();
      string line = "";
      FileStream ifs = new FileStream(fn, FileMode.Open);
      StreamReader sr = new StreamReader(ifs);
      while ((line = sr.ReadLine()) != null)
      {
        if (line.StartsWith(comment) == true)
          continue;
        string[] tokens = line.Split(sep);
        List"lt"double"gt" lst = new List"lt"double"gt"();
        for (int j = 0; j "lt" usecols.Length; ++j)
          lst.Add(double.Parse(tokens[usecols[j]]));
        double[] row = lst.ToArray();
        result.Add(row);
      }
      sr.Close(); ifs.Close();
      return result.ToArray();
    }

    static double[] MatToVec(double[][] mat)
    {
      int nRows = mat.Length;
      int nCols = mat[0].Length;
      double[] result = new double[nRows * nCols];
      int k = 0;
      for (int i = 0; i "lt" nRows; ++i)
        for (int j = 0; j "lt" nCols; ++j)
          result[k++] = mat[i][j];
      return result;
    }

    static void VecShow(double[] vec, int dec, int wid)
    {
      for (int i = 0; i "lt" vec.Length; ++i)
        Console.Write(vec[i].ToString("F" + dec).
          PadLeft(wid));
      Console.WriteLine("");
    }

  } // class Program

  // ========================================================

  public class GradientBoostRegressor
  {
    public double lrnRate;
    public int nTrees;
    public int maxDepth;
    public int minSamples;
    public int minLeaf;
    public int numSplitCols;
    public List"lt"DecisionTreeRegressor"gt" trees;
    public double pred0;  // initial prediction
    private Random rnd;

    public GradientBoostRegressor(int nTrees, int maxDepth,
      int minSamples, int minLeaf, int numSplitCols,
      double lrnRate, int seed = 0)
    {
      this.nTrees = nTrees;
      this.maxDepth = maxDepth;
      this.minSamples = minSamples;
      this.minLeaf = minLeaf;
      this.numSplitCols = numSplitCols;
      this.trees = new List"lt"DecisionTreeRegressor"gt"();
      this.lrnRate = lrnRate;
      this.rnd = new Random(seed);
    }

    // ------------------------------------------------------

    public void Train(double[][] trainX, double[] trainY)
    {
      int n = trainX.Length;
      this.pred0 = Mean(trainY);

      double[] preds = new double[n]; //each data item
      for (int i = 0; i "lt" n; ++i)
        preds[i] = this.pred0;

      double[] residuals = new double[n]; // for curr tree
      for (int t = 0; t "lt" this.nTrees; ++t) // each tree
      {
        // compute residuals from prev tree predictions
        for (int i = 0; i "lt" n; ++i)
          residuals[i] = trainY[i] - preds[i];

        int treeSeed = this.rnd.Next();
        DecisionTreeRegressor dtr =
          new DecisionTreeRegressor(
            this.maxDepth,
            this.minSamples,
            this.minLeaf,
            this.numSplitCols,
            false, // don't save rows info
            treeSeed);
        dtr.Train(trainX, residuals); // predict residuals

        // predict the residuals
        for (int i = 0; i "lt" n; ++i)
        {
          double predResidual = dtr.Predict(trainX[i]);
          preds[i] += this.lrnRate * predResidual;
        }
        this.trees.Add(dtr);
      }
    } // Train

    // ------------------------------------------------------

    public double Predict(double[] x)
    {
      double result = this.pred0;
      for (int t = 0; t "lt" this.nTrees; ++t)
      {
        double predResidual = this.trees[t].Predict(x);
        double delta = this.lrnRate * predResidual;
        result += delta;
      }
      return result;
    }

    // ------------------------------------------------------

    public double Explain(double[] x)
    {
      double result = this.pred0;

      Console.WriteLine("\nInitial prediction: " +
      result.ToString("F4"));

      for (int t = 0; t "lt" this.nTrees; ++t)
      {
        double predResidual = this.trees[t].Predict(x);
        double delta = this.lrnRate * predResidual;
        result += delta;

        if (t "gte" 0 && t "lte" 4 ||
          t "gte" this.nTrees - 2 &&
          t "lte" this.nTrees - 1)
        {
          Console.Write("t = " + t.ToString().PadLeft(3) +
            "  pred_res = " + predResidual.ToString("F4").
            PadLeft(10));
          Console.Write("  delta = " + delta.ToString("F4").
            PadLeft(10));
          Console.WriteLine("  pred = " +
            result.ToString("F4").PadLeft(7));
        }
        if (t == 5)
          Console.WriteLine(". . . ");
      }
      return result;
    }

    // ------------------------------------------------------

    public double Accuracy(double[][] dataX, double[] dataY,
      double pctClose)
    {
      int numCorrect = 0; int numWrong = 0;
      for (int i = 0; i "lt" dataX.Length; ++i)
      {
        double actualY = dataY[i];
        double predY = Predict(dataX[i]);
        if (Math.Abs(predY - actualY) "lt"
          (pctClose * Math.Abs(actualY)))
          ++numCorrect;
        else
          ++numWrong;
      }
      return (numCorrect * 1.0) / (numWrong + numCorrect);
    }

    // ------------------------------------------------------

    public double MSE(double[][] dataX,
      double[] dataY)
    {
      int n = dataX.Length;
      double sum = 0.0;
      for (int i = 0; i "lt" n; ++i)
      {
        double actualY = dataY[i];
        double predY = this.Predict(dataX[i]);
        sum += (actualY - predY) * (actualY - predY);
      }
      return sum / n;
    }

    // ------------------------------------------------------

    public double R2(double[][] dataX, double[] dataY)
    {
      // coefficient of determination
      int n = dataX.Length;
      double sum = 0.0;

      for (int i = 0; i "lt" n; ++i)
        sum += dataY[i];
      double meanActual = sum / n;

      double ssRes = 0.0;
      double ssTot = 0.0;
      for (int i = 0; i "lt" n; ++i)
      {
        double predY = this.Predict(dataX[i]);
        ssRes += (dataY[i] - predY) * (dataY[i] - predY);
        ssTot += (dataY[i] - meanActual) *
          (dataY[i] - meanActual);
      }
      if (Math.Abs(ssTot) "lt" 1.0e-12)
        return 0.0;
      else
        return 1.0 - (ssRes / ssTot);
    }

    // ------------------------------------------------------

    private static double Mean(double[] data)
    {
      int n = data.Length;
      double sum = 0.0;
      for (int i = 0; i "lt" n; ++i)
        sum += data[i];
      return sum / n;
    }

  } // class GradientBoostRegression

  // ========================================================
  // ========================================================

  public class DecisionTreeRegressor
  {
    // ChatGPT: "I don't see any fundamental tree-building or
    // indexing bugs. If the inputs are guaranteed to be
    // finite doubles and the parameters are validated 
    // elsewhere, the core regression algorithm is correct."

    public int maxDepth;
    public int minSamples;
    public int minLeaf;
    public int numSplitCols;
    public List"lt"Node"gt" tree;
    public Random rnd;
    public bool saveRows;
    public double[][] trainX;
    public double[] trainY;

    // ............................................

    public class Node
    {
      public int id;
      public int colIdx = -1;
      public double thresh = 0.0;
      public int left = -1;
      public int right = -1;
      public double value = 0.0;
      public bool isLeaf = false;
      public List"lt"int"gt" rows;
    }

    // ............................................

    public DecisionTreeRegressor(int maxDepth = 3,
      int minSamples = 2, int minLeaf = 1,
      int numSplitCols = -1, bool saveRows = false,
      int seed = 0)
    {
      this.maxDepth = maxDepth;
      this.minSamples = minSamples;
      this.minLeaf = minLeaf;
      this.numSplitCols = numSplitCols;
      this.saveRows = saveRows;
      this.tree = new List"lt"Node"gt"();

      int numNodes = (int)Math.Pow(2, (maxDepth + 1)) - 1;
      for (int i = 0; i "lt" numNodes; ++i)
        this.tree.Add(null);
      this.rnd = new Random(seed);
    }

    // ------------------------------------------------------
    // ------------------------------------------------------
    // public: ctor(), Train(), Predict()
    // private: BestSplit(), TreeTargetMean()
    // ------------------------------------------------------

    public void Train(double[][] trainX, double[] trainY)
    {
      this.trainX = trainX;
      this.trainY = trainY;

      int maxParentID = (int)Math.Pow(2, this.maxDepth) - 2;

      List"lt"int"gt" allRows =
        new List"lt"int"gt"(this.trainX.Length);
      for (int i = 0; i "lt" this.trainX.Length; ++i)
        allRows.Add(i);

      // construct root node
      Node root = new Node()
      {
        id = 0,
        value = this.TreeTargetMean(allRows),
        rows = allRows
      };
      this.tree[0] = root;

      for (int i = 0; i "lt" this.tree.Count; ++i)
      {
        Node currNode = this.tree[i];
        if (currNode == null) continue;

        int leftID = currNode.id * 2 + 1;
        int rightID = currNode.id * 2 + 2;

        // 1. check before expensive split calcs
        if (currNode.id "gt" maxParentID ||
            currNode.rows.Count "lt" this.minSamples ||
            rightID "gte" this.tree.Count)
        {
          currNode.isLeaf = true;
          continue;
        }

        // 2. find best split
        double[] splitInfo = this.BestSplit(currNode.rows);
        int colIdx = (int)splitInfo[0];
        double splitVal = splitInfo[1];

        if (colIdx == -1) // unable to split
        {
          currNode.isLeaf = true;
          continue;
        }

        // 3. apply split properties
        currNode.colIdx = colIdx;  // aka split col
        currNode.thresh = splitVal;

        List"lt"int"gt" leftIdxs = new List"lt"int"gt"();
        List"lt"int"gt" rightIdxs = new List"lt"int"gt"();

        for (int k = 0; k "lt" currNode.rows.Count; ++k)
        {
          int r = currNode.rows[k];
          if (this.trainX[r][colIdx] "lte" splitVal)
            leftIdxs.Add(r);
          else
            rightIdxs.Add(r);
        }

        // 4. create and attach children
        currNode.left = leftID;
        this.tree[leftID] = new Node()
        {
          id = leftID,
          rows = leftIdxs,
          value = this.TreeTargetMean(leftIdxs)
        };

        currNode.right = rightID;
        this.tree[rightID] = new Node()
        {
          id = rightID,
          rows = rightIdxs,
          value = this.TreeTargetMean(rightIdxs)
        };
      } // main loop

      // delete row info (for ensembles)
      if (this.saveRows == false)
      {
        for (int i = 0; i "lt" this.tree.Count; ++i)
          if (this.tree[i] != null)
            this.tree[i].rows = null;
      }
    } // Train()

    // ------------------------------------------------------

    public double Predict(double[] x)
    {
      int p = 0;
      double lastValidValue = 0.0;
      while (p != -1 && p "lt" this.tree.Count)
      {
        Node currNode = this.tree[p];
        if (currNode == null) break;

        lastValidValue = currNode.value;
        if (currNode.isLeaf || 
          currNode.colIdx == -1) break;

        if (x[currNode.colIdx] "lte" currNode.thresh)
          p = currNode.left;
        else
          p = currNode.right;
      }
      return lastValidValue;
    }

    // ------------------------------------------------------

    private double[] BestSplit(List"lt"int"gt" rows)
    {
      int bestColIdx = -1;
      double bestThresh = 0.0;
      double bestGain = -1.0;
      int nRows = rows.Count;
      int nCols = this.trainX[0].Length;

      int[] colIndices = new int[nCols];
      for (int k = 0; k "lt" nCols; ++k)
        colIndices[k] = k;

      // feature sampling
      for (int i = 0; i "lt" nCols; ++i)
      {
        int ri = rnd.Next(i, nCols);
        int tmp = colIndices[i];
        colIndices[i] = colIndices[ri];
        colIndices[ri] = tmp;
      }

      int numColsToUse;
      if (this.numSplitCols == -1)
        numColsToUse = nCols;
      else
        numColsToUse = Math.Min(this.numSplitCols, nCols);

      int[] sortedRows = new int[nRows];
      double[] featureKeys = new double[nRows];

      double totalSum = 0.0;
      for (int i = 0; i "lt" nRows; ++i)
        totalSum += this.trainY[rows[i]];

      for (int j = 0; j "lt" numColsToUse; ++j)
      {
        int colIdx = colIndices[j];
        for (int i = 0; i "lt" nRows; ++i)
        {
          int r = rows[i];
          sortedRows[i] = r;
          featureKeys[i] = this.trainX[r][colIdx];
        }

        Array.Sort(featureKeys, sortedRows);
        double leftSum = 0.0;

        for (int i = 0; i "lt" nRows - 1; ++i)
        {
          leftSum += this.trainY[sortedRows[i]];

          int leftCount = i + 1;
          int rightCount = nRows - leftCount;

          if (featureKeys[i] == featureKeys[i + 1])
            continue; // duplicate values
          if (leftCount "lt" this.minLeaf ||
            rightCount "lt" this.minLeaf)
            continue;

          double rightSum = totalSum - leftSum;

          // maximizing Gain = (G_L^2 / N_L) + (G_R^2 / N_R)
          // equivalent to maximizing explicit weighted
          // variance reduction
          double gain = (leftSum * leftSum / leftCount) +
            (rightSum * rightSum / rightCount);

          if (gain "gt" bestGain)
          {
            bestGain = gain;
            bestColIdx = colIdx;
            // bestThresh = 
            //  (featureKeys[i] + featureKeys[i + 1]) / 2.0;
            bestThresh = featureKeys[i] +
              (featureKeys[i + 1] - featureKeys[i]) / 2.0;
          }
        }
      }

      // using a tuple would create version dependency
      double[] result = new double[2];
      result[0] = 1.0 * bestColIdx;
      result[1] = bestThresh;
      return result;
    }

    private double TreeTargetMean(List"lt"int"gt" rows)
    {
      if (rows == null || rows.Count == 0)
        return 0.0;
      double sum = 0.0;
      for (int i = 0; i "lt" rows.Count; ++i)
        sum += this.trainY[rows[i]];
      return sum / rows.Count;
    }

  } // DecisionTreeRegressor

  // ========================================================

} // ns

Training data:

# synthetic_train_200.txt
#
-0.1660,  0.4406, -0.9998, -0.3953, -0.7065,  0.4840
 0.0776, -0.1616,  0.3704, -0.5911,  0.7562,  0.1568
-0.9452,  0.3409, -0.1654,  0.1174, -0.7192,  0.8054
 0.9365, -0.3732,  0.3846,  0.7528,  0.7892,  0.1345
-0.8299, -0.9219, -0.6603,  0.7563, -0.8033,  0.7955
 0.0663,  0.3838, -0.3690,  0.3730,  0.6693,  0.3206
-0.9634,  0.5003,  0.9777,  0.4963, -0.4391,  0.7377
-0.1042,  0.8172, -0.4128, -0.4244, -0.7399,  0.4801
-0.9613,  0.3577, -0.5767, -0.4689, -0.0169,  0.6861
-0.7065,  0.1786,  0.3995, -0.7953, -0.1719,  0.5569
 0.3888, -0.1716, -0.9001,  0.0718,  0.3276,  0.2500
 0.1731,  0.8068, -0.7251, -0.7214,  0.6148,  0.3297
-0.2046, -0.6693,  0.8550, -0.3045,  0.5016,  0.2129
 0.2473,  0.5019, -0.3022, -0.4601,  0.7918,  0.2613
-0.1438,  0.9297,  0.3269,  0.2434, -0.7705,  0.5171
 0.1568, -0.1837, -0.5259,  0.8068,  0.1474,  0.3307
-0.9943,  0.2343, -0.3467,  0.0541,  0.7719,  0.5581
 0.2467, -0.9684,  0.8589,  0.3818,  0.9946,  0.1092
-0.6553, -0.7257,  0.8652,  0.3936, -0.8680,  0.7018
 0.8460,  0.4230, -0.7515, -0.9602, -0.9476,  0.1996
-0.9434, -0.5076,  0.7201,  0.0777,  0.1056,  0.5664
 0.9392,  0.1221, -0.9627,  0.6013, -0.5341,  0.1533
 0.6142, -0.2243,  0.7271,  0.4942,  0.1125,  0.1661
 0.4260,  0.1194, -0.9749, -0.8561,  0.9346,  0.2230
 0.1362, -0.5934, -0.4953,  0.4877, -0.6091,  0.3810
 0.6937, -0.5203, -0.0125,  0.2399,  0.6580,  0.1460
-0.6864, -0.9628, -0.8600, -0.0273,  0.2127,  0.5387
 0.9772,  0.1595, -0.2397,  0.1019,  0.4907,  0.1611
 0.3385, -0.4702, -0.8673, -0.2598,  0.2594,  0.2270
-0.8669, -0.4794,  0.6095, -0.6131,  0.2789,  0.4700
 0.0493,  0.8496, -0.4734, -0.8681,  0.4701,  0.3516
 0.8639, -0.9721, -0.5313,  0.2336,  0.8980,  0.1412
 0.9004,  0.1133,  0.8312,  0.2831, -0.2200,  0.1782
 0.0991,  0.8524,  0.8375, -0.2102,  0.9265,  0.2150
-0.6521, -0.7473, -0.7298,  0.0113, -0.9570,  0.7422
 0.6190, -0.3105,  0.8802,  0.1640,  0.7577,  0.1056
 0.6895,  0.8108, -0.0802,  0.0927,  0.5972,  0.2214
 0.1982, -0.9689,  0.1870, -0.1326,  0.6147,  0.1310
-0.3695,  0.7858,  0.1557, -0.6320,  0.5759,  0.3773
-0.1596,  0.3581,  0.8372, -0.9992,  0.9535,  0.2071
-0.2468,  0.9476,  0.2094,  0.6577,  0.1494,  0.4132
 0.1737,  0.5000,  0.7166,  0.5102,  0.3961,  0.2611
 0.7290, -0.3546,  0.3416, -0.0983, -0.2358,  0.1332
-0.3652,  0.2438, -0.1395,  0.9476,  0.3556,  0.4170
-0.6029, -0.1466, -0.3133,  0.5953,  0.7600,  0.4334
-0.4596, -0.4953,  0.7098,  0.0554,  0.6043,  0.2775
 0.1450,  0.4663,  0.0380,  0.5418,  0.1377,  0.2931
-0.8636, -0.2442, -0.8407,  0.9656, -0.6368,  0.7429
 0.6237,  0.7499,  0.3768,  0.1390, -0.6781,  0.2185
-0.5499,  0.1850, -0.3755,  0.8326,  0.8193,  0.4399
-0.4858, -0.7782, -0.6141, -0.0008,  0.4572,  0.4197
 0.7033, -0.1683,  0.2334, -0.5327, -0.7961,  0.1776
 0.0317, -0.0457, -0.6947,  0.2436,  0.0880,  0.3345
 0.5031, -0.5559,  0.0387,  0.5706, -0.9553,  0.3107
-0.3513,  0.7458,  0.6894,  0.0769,  0.7332,  0.3170
 0.2205,  0.5992, -0.9309,  0.5405,  0.4635,  0.3532
-0.4806, -0.4859,  0.2646, -0.3094,  0.5932,  0.3202
 0.9809, -0.3995, -0.7140,  0.8026,  0.0831,  0.1600
 0.9495,  0.2732,  0.9878,  0.0921,  0.0529,  0.1289
-0.9476, -0.6792,  0.4913, -0.9392, -0.2669,  0.5966
 0.7247,  0.3854,  0.3819, -0.6227, -0.1162,  0.1550
-0.5922, -0.5045, -0.4757,  0.5003, -0.0860,  0.5863
-0.8861,  0.0170, -0.5761,  0.5972, -0.4053,  0.7301
 0.6877, -0.2380,  0.4997,  0.0223,  0.0819,  0.1404
 0.9189,  0.6079, -0.9354,  0.4188, -0.0700,  0.1907
-0.1428, -0.7820,  0.2676,  0.6059,  0.3936,  0.2790
 0.5324, -0.3151,  0.6917, -0.1425,  0.6480,  0.1071
-0.8432, -0.9633, -0.8666, -0.0828, -0.7733,  0.7784
-0.9444,  0.5097, -0.2103,  0.4939, -0.0952,  0.6787
-0.0520,  0.6063, -0.1952,  0.8094, -0.9259,  0.4836
 0.5477, -0.7487,  0.2370, -0.9793,  0.0773,  0.1241
 0.2450,  0.8116,  0.9799,  0.4222,  0.4636,  0.2355
 0.8186, -0.1983, -0.5003, -0.6531, -0.7611,  0.1511
-0.4714,  0.6382, -0.3788,  0.9648, -0.4667,  0.5950
 0.0673, -0.3711,  0.8215, -0.2669, -0.1328,  0.2677
-0.9381,  0.4338,  0.7820, -0.9454,  0.0441,  0.5518
-0.3480,  0.7190,  0.1170,  0.3805, -0.0943,  0.4724
-0.9813,  0.1535, -0.3771,  0.0345,  0.8328,  0.5438
-0.1471, -0.5052, -0.2574,  0.8637,  0.8737,  0.3042
-0.5454, -0.3712, -0.6505,  0.2142, -0.1728,  0.5783
 0.6327, -0.6297,  0.4038, -0.5193,  0.1484,  0.1153
-0.5424,  0.3282, -0.0055,  0.0380, -0.6506,  0.6613
 0.1414,  0.9935,  0.6337,  0.1887,  0.9520,  0.2540
-0.9351, -0.8128, -0.8693, -0.0965, -0.2491,  0.7353
 0.9507, -0.6640,  0.9456,  0.5349,  0.6485,  0.1059
-0.0462, -0.9737, -0.2940, -0.0159,  0.4602,  0.2606
-0.0627, -0.0852, -0.7247, -0.9782,  0.5166,  0.2977
 0.0478,  0.5098, -0.0723, -0.7504, -0.3750,  0.3335
 0.0090,  0.3477,  0.5403, -0.7393, -0.9542,  0.4415
-0.9748,  0.3449,  0.3736, -0.1015,  0.8296,  0.4358
 0.2887, -0.9895, -0.0311,  0.7186,  0.6608,  0.2057
 0.1570, -0.4518,  0.1211,  0.3435, -0.2951,  0.3244
 0.7117, -0.6099,  0.4946, -0.4208,  0.5476,  0.1096
-0.2929, -0.5726,  0.5346, -0.3827,  0.4665,  0.2465
 0.4889, -0.5572, -0.5718, -0.6021, -0.7150,  0.2163
-0.7782,  0.3491,  0.5996, -0.8389, -0.5366,  0.6516
-0.5847,  0.8347,  0.4226,  0.1078, -0.3910,  0.6134
 0.8469,  0.4121, -0.0439, -0.7476,  0.9521,  0.1571
-0.6803, -0.5948, -0.1376, -0.1916, -0.7065,  0.7156
 0.2878,  0.5086, -0.5785,  0.2019,  0.4979,  0.2980
 0.2764,  0.1943, -0.4090,  0.4632,  0.8906,  0.2960
-0.8877,  0.6705, -0.6155, -0.2098, -0.3998,  0.7107
-0.8398,  0.8093, -0.2597,  0.0614, -0.0118,  0.6502
-0.8476,  0.0158, -0.4769, -0.2859, -0.7839,  0.7715
 0.5751, -0.7868,  0.9714, -0.6457,  0.1448,  0.1175
 0.4802, -0.7001,  0.1022, -0.5668,  0.5184,  0.1090
 0.4458, -0.6469,  0.7239, -0.9604,  0.7205,  0.0779
 0.5175,  0.4339,  0.9747, -0.4438, -0.9924,  0.2879
 0.8678,  0.7158,  0.4577,  0.0334,  0.4139,  0.1678
 0.5406,  0.5012,  0.2264, -0.1963,  0.3946,  0.2088
-0.9938,  0.5498,  0.7928, -0.5214, -0.7585,  0.7687
 0.7661,  0.0863, -0.4266, -0.7233, -0.4197,  0.1466
 0.2277, -0.3517, -0.0853, -0.1118,  0.6563,  0.1767
 0.3499, -0.5570, -0.0655, -0.3705,  0.2537,  0.1632
 0.7547, -0.1046,  0.5689, -0.0861,  0.3125,  0.1257
 0.8186,  0.2110,  0.5335,  0.0094, -0.0039,  0.1391
 0.6858, -0.8644,  0.1465,  0.8855,  0.0357,  0.1845
-0.4967,  0.4015,  0.0805,  0.8977,  0.2487,  0.4663
 0.6760, -0.9841,  0.9787, -0.8446, -0.3557,  0.1509
-0.1203, -0.4885,  0.6054, -0.0443, -0.7313,  0.4854
 0.8557,  0.7919, -0.0169,  0.7134, -0.1628,  0.2002
 0.0115, -0.6209,  0.9300, -0.4116, -0.7931,  0.4052
-0.7114, -0.9718,  0.4319,  0.1290,  0.5892,  0.3661
 0.3915,  0.5557, -0.1870,  0.2955, -0.6404,  0.2954
-0.3564, -0.6548, -0.1827, -0.5172, -0.1862,  0.4622
 0.2392, -0.4959,  0.5857, -0.1341, -0.2850,  0.2470
-0.3394,  0.3947, -0.4627,  0.6166, -0.4094,  0.5325
 0.7107,  0.7768, -0.6312,  0.1707,  0.7964,  0.2757
-0.1078,  0.8437, -0.4420,  0.2177,  0.3649,  0.4028
-0.3139,  0.5595, -0.6505, -0.3161, -0.7108,  0.5546
 0.4335,  0.3986,  0.3770, -0.4932,  0.3847,  0.1810
-0.2562, -0.2894, -0.8847,  0.2633,  0.4146,  0.4036
 0.2272,  0.2966, -0.6601, -0.7011,  0.0284,  0.2778
-0.0743, -0.1421, -0.0054, -0.6770, -0.3151,  0.3597
-0.4762,  0.6891,  0.6007, -0.1467,  0.2140,  0.4266
-0.4061,  0.7193,  0.3432,  0.2669, -0.7505,  0.6147
-0.0588,  0.9731,  0.8966,  0.2902, -0.6966,  0.4955
-0.0627, -0.1439,  0.1985,  0.6999,  0.5022,  0.3077
 0.1587,  0.8494, -0.8705,  0.9827, -0.8940,  0.4263
-0.7850,  0.2473, -0.9040, -0.4308, -0.8779,  0.7199
 0.4070,  0.3369, -0.2428, -0.6236,  0.4940,  0.2215
-0.0242,  0.0513, -0.9430,  0.2885, -0.2987,  0.3947
-0.5416, -0.1322, -0.2351, -0.0604,  0.9590,  0.3683
 0.1055,  0.7783, -0.2901, -0.5090,  0.8220,  0.2984
-0.9129,  0.9015,  0.1128, -0.2473,  0.9901,  0.4776
-0.9378,  0.1424, -0.6391,  0.2619,  0.9618,  0.5368
 0.7498, -0.0963,  0.4169,  0.5549, -0.0103,  0.1614
-0.2612, -0.7156,  0.4538, -0.0460, -0.1022,  0.3717
 0.7720,  0.0552, -0.1818, -0.4622, -0.8560,  0.1685
-0.4177,  0.0070,  0.9319, -0.7812,  0.3461,  0.3052
-0.0001,  0.5542, -0.7128, -0.8336, -0.2016,  0.3803
 0.5356, -0.4194, -0.5662, -0.9666, -0.2027,  0.1776
-0.2378,  0.3187, -0.8582, -0.6948, -0.9668,  0.5474
-0.1947, -0.3579,  0.1158,  0.9869,  0.6690,  0.2992
 0.3992,  0.8365, -0.9205, -0.8593, -0.0520,  0.3154
-0.0209,  0.0793,  0.7905, -0.1067,  0.7541,  0.1864
-0.4928, -0.4524, -0.3433,  0.0951, -0.5597,  0.6261
-0.8118,  0.7404, -0.5263, -0.2280,  0.1431,  0.6349
 0.0516, -0.8480,  0.7483,  0.9023,  0.6250,  0.1959
-0.3212,  0.1093,  0.9488, -0.3766,  0.3376,  0.2735
-0.3481,  0.5490, -0.3484,  0.7797,  0.5034,  0.4379
-0.5785, -0.9170, -0.3563, -0.9258,  0.3877,  0.4121
 0.3407, -0.1391,  0.5356,  0.0720, -0.9203,  0.3458
-0.3287, -0.8954,  0.2102,  0.0241,  0.2349,  0.3247
-0.1353,  0.6954, -0.0919, -0.9692,  0.7461,  0.3338
 0.9036, -0.8982, -0.5299, -0.8733, -0.1567,  0.1187
 0.7277, -0.8368, -0.0538, -0.7489,  0.5458,  0.0830
 0.9049,  0.8878,  0.2279,  0.9470, -0.3103,  0.2194
 0.7957, -0.1308, -0.5284,  0.8817,  0.3684,  0.2172
 0.4647, -0.4931,  0.2010,  0.6292, -0.8918,  0.3371
-0.7390,  0.6849,  0.2367,  0.0626, -0.5034,  0.7039
-0.1567, -0.8711,  0.7940, -0.5932,  0.6525,  0.1710
 0.7635, -0.0265,  0.1969,  0.0545,  0.2496,  0.1445
 0.7675,  0.1354, -0.7698, -0.5460,  0.1920,  0.1728
-0.5211, -0.7372, -0.6763,  0.6897,  0.2044,  0.5217
 0.1913,  0.1980,  0.2314, -0.8816,  0.5006,  0.1998
 0.8964,  0.0694, -0.6149,  0.5059, -0.9854,  0.1825
 0.1767,  0.7104,  0.2093,  0.6452,  0.7590,  0.2832
-0.3580, -0.7541,  0.4426, -0.1193, -0.7465,  0.5657
-0.5996,  0.5766, -0.9758, -0.3933, -0.9572,  0.6800
 0.9950,  0.1641, -0.4132,  0.8579,  0.0142,  0.2003
-0.4717, -0.3894, -0.2567, -0.5111,  0.1691,  0.4266
 0.3917, -0.8561,  0.9422,  0.5061,  0.6123,  0.1212
-0.0366, -0.1087,  0.3449, -0.1025,  0.4086,  0.2475
 0.3633,  0.3943,  0.2372, -0.6980,  0.5216,  0.1925
-0.5325, -0.6466, -0.2178, -0.3589,  0.6310,  0.3568
 0.2271,  0.5200, -0.1447, -0.8011, -0.7699,  0.3128
 0.6415,  0.1993,  0.3777, -0.0178, -0.8237,  0.2181
-0.5298, -0.0768, -0.6028, -0.9490,  0.4588,  0.4356
 0.6870, -0.1431,  0.7294,  0.3141,  0.1621,  0.1632
-0.5985,  0.0591,  0.7889, -0.3900,  0.7419,  0.2945
 0.3661,  0.7984, -0.8486,  0.7572, -0.6183,  0.3449
 0.6995,  0.3342, -0.3113, -0.6972,  0.2707,  0.1712
 0.2565,  0.9126,  0.1798, -0.6043, -0.1413,  0.2893
-0.3265,  0.9839, -0.2395,  0.9854,  0.0376,  0.4770
 0.2690, -0.1722,  0.9818,  0.8599, -0.7015,  0.3954
-0.2102, -0.0768,  0.1219,  0.5607, -0.0256,  0.3949
 0.8216, -0.9555,  0.6422, -0.6231,  0.3715,  0.0801
-0.2896,  0.9484, -0.7545, -0.6249,  0.7789,  0.4370
-0.9985, -0.5448, -0.7092, -0.5931,  0.7926,  0.5402

Test data:

# synthetic_test_40.txt
#
 0.7462,  0.4006, -0.0590,  0.6543, -0.0083,  0.1935
 0.8495, -0.2260, -0.0142, -0.4911,  0.7699,  0.1078
-0.2335, -0.4049,  0.4352, -0.6183, -0.7636,  0.5088
 0.1810, -0.5142,  0.2465,  0.2767, -0.3449,  0.3136
-0.8650,  0.7611, -0.0801,  0.5277, -0.4922,  0.7140
-0.2358, -0.7466, -0.5115, -0.8413, -0.3943,  0.4533
 0.4834,  0.2300,  0.3448, -0.9832,  0.3568,  0.1360
-0.6502, -0.6300,  0.6885,  0.9652,  0.8275,  0.3046
-0.3053,  0.5604,  0.0929,  0.6329, -0.0325,  0.4756
-0.7995,  0.0740, -0.2680,  0.2086,  0.9176,  0.4565
-0.2144, -0.2141,  0.5813,  0.2902, -0.2122,  0.4119
-0.7278, -0.0987, -0.3312, -0.5641,  0.8515,  0.4438
 0.3793,  0.1976,  0.4933,  0.0839,  0.4011,  0.1905
-0.8568,  0.9573, -0.5272,  0.3212, -0.8207,  0.7415
-0.5785,  0.0056, -0.7901, -0.2223,  0.0760,  0.5551
 0.0735, -0.2188,  0.3925,  0.3570,  0.3746,  0.2191
 0.1230, -0.2838,  0.2262,  0.8715,  0.1938,  0.2878
 0.4792, -0.9248,  0.5295,  0.0366, -0.9894,  0.3149
-0.4456,  0.0697,  0.5359, -0.8938,  0.0981,  0.3879
 0.8629, -0.8505, -0.4464,  0.8385,  0.5300,  0.1769
 0.1995,  0.6659,  0.7921,  0.9454,  0.9970,  0.2330
-0.0249, -0.3066, -0.2927, -0.4923,  0.8220,  0.2437
 0.4513, -0.9481, -0.0770, -0.4374, -0.9421,  0.2879
-0.3405,  0.5931, -0.3507, -0.3842,  0.8562,  0.3987
 0.9538,  0.0471,  0.9039,  0.7760,  0.0361,  0.1706
-0.0887,  0.2104,  0.9808,  0.5478, -0.3314,  0.4128
-0.8220, -0.6302,  0.0537, -0.1658,  0.6013,  0.4306
-0.4123, -0.2880,  0.9074, -0.0461, -0.4435,  0.5144
 0.0060,  0.2867, -0.7775,  0.5161,  0.7039,  0.3599
-0.7968, -0.5484,  0.9426, -0.4308,  0.8148,  0.2979
 0.7811,  0.8450, -0.6877,  0.7594,  0.2640,  0.2362
-0.6802, -0.1113, -0.8325, -0.6694, -0.6056,  0.6544
 0.3821,  0.1476,  0.7466, -0.5107,  0.2592,  0.1648
 0.7265,  0.9683, -0.9803, -0.4943, -0.5523,  0.2454
-0.9049, -0.9797, -0.0196, -0.9090, -0.4433,  0.6447
-0.4607,  0.1811, -0.2389,  0.4050, -0.0078,  0.5229
 0.2664, -0.2932, -0.4259, -0.7336,  0.8742,  0.1834
-0.4507,  0.1029, -0.6294, -0.1158, -0.6294,  0.6081
 0.8948, -0.0124,  0.9278,  0.2899, -0.0314,  0.1534
-0.1323, -0.8813, -0.0146, -0.0697,  0.6135,  0.2386

Bottom line: I decided to implement random forest regression using Extra (“Extremely Randomized”) Tree learners. For my demo dataset, the technique worked slightly better, but not significantly better, than the standard random forest architecture that uses regular decision trees as the learners.

As is often the case with complex problems, explaining what the problem is, is more difficult than explaining the solution. So bear with me.

The goal of a machine learning regression problem is to predict a single numeric value. For example, a bank might want to predict the maximum safe loan amount to a customer, based on age, account balance, current debt, and so on.

There are many regression techniques. The techniques fall into two main categories: 1.) classical math-based techniques (linear regression, nearest neighbors regression, quadratic regression, kernel ridge regression, neural network regression, and others) and 2.) tree-based techniques (random forest regression, Extra Trees (“extremely randomized”) regression, AdaBoost regression, Gradient Boost regression, and others).

The tree-based random forest regression technique uses a collection of simple decision trees — they’re called the learners or the estimators. Each tree is trained using a different subset of the source training data. The rows are selected with replacement so duplicates are common. The columns are selected without replacement so there are no duplicate columns in a subset.

The final prediction is average of the predictions of the trees. The idea is that a single tree regressor will overfit the training data, but a collection of weak tree regressors will average out to a good result.

Although a random forest regression almost always uses standard decision trees as the learners, in theory, any kind of simple regression technique can be used. Paradoxically, the base learners need to be weak instead of powerful, but that’s a long and complicated story.

Extra Tree regressors are even weaker than standard decision trees. When constructing a standard decision tree, for each tree node, an optimal split column and split value are found. When constructing an Extra Tree tree, for each node, a random split value is selected from each column, and then the best of those is used. So Extra Tree trees train much faster than regular decision trees, have a kind of built-in regularization that, in theory, might help the overall random forest prediction model.

For my demo, I used one of my standard datasets It looks like:

-0.1660,  0.4406, -0.9998, -0.3953, -0.7065,  0.4840
 0.0776, -0.1616,  0.3704, -0.5911,  0.7562,  0.1568
-0.9452,  0.3409, -0.1654,  0.1174, -0.7192,  0.8054
. . .

The data is synthetic. The first five values on each line are the predictors. The last value on each line is the target to predict. There are 200 training items and 40 test items.

The key parts of the output of my demo are:

Begin C# Random Forest with Extra Trees regression demo

Loading synthetic train (200) and test (40) data
Done

First three train X:
 -0.1660  0.4406 -0.9998 -0.3953 -0.7065
  0.0776 -0.1616  0.3704 -0.5911  0.7562
 -0.9452  0.3409 -0.1654  0.1174 -0.7192

First three train y:
  0.4840
  0.1568
  0.8054

Setting nTrees = 400
Setting maxDepth = 7
Setting minSamples = 2
Setting minLeaf = 1
Setting numSplitCols = 4
Setting nRows = 180

Creating and training RandomForestRegressor model
Done

Accuracy train (within 0.10) = 0.8100
Accuracy test (within 0.10) = 0.7000

MSE train = 0.0004
MSE test = 0.0011

Predicting for x =
  -0.1660   0.4406  -0.9998  -0.3953  -0.7065
Predicted y = 0.4860

End demo

The random forest regressor with Extra Tree learners scored 81.00% accuracy on the training data (162 out of 200 correct), and 70.00% accuracy on the test data (28 out of 40 correct). A prediction is scored correct if it’s within 10% of the true target value. The corresponding accuracies when using standard decision trees are 83.50% (train) and 65.00% (test). Because the test accuracy is more important, the random forest with Extra Tree learners scored slightly better than the standard random forest with decision tree learners. But I do not rate the difference as significant.

A very interesting experiment.

Random forest regression code is a wrapper around a collection of weak learners. I’m a big fan of science fiction movies of the 1950s. Many of these movies had creatures that were actors wrapped in costumes of some kind. Even though the costumes were usually weak and not realistic, I still liked many of these movies. Here are two that have a aquatic theme.

Right: In “The Phantom from 10,000 Leagues” (1955), a rogue scientist creates a sea creature using radiation (naturally — this was the 1050s). The creature kills several people. The scientist eventually realizes that this experiment is not going well and he blows up his lab, the creature, and himself. This is a very low-budget effort but the actors do a decent job with limited plot and sets. My grade = C (but my grading is very lenient).

Demo program. Very long, quite complex. Replace “lt” (less than), “gt”, “lte”, “gte” with Boolean operator symbols. (My blog editor chokes on symbols).

using System;
using System.Collections.Generic;
using System.IO;

namespace RandomForestWithExtraTrees
{
  internal class RandomForestWithExtraTreesProgram
  {
    static void Main(string[] args)
    {
      Console.WriteLine("\nBegin C# Random Forest" +
        " with Extra Trees regression demo ");

      // 1. load data
      Console.WriteLine("\nLoading synthetic train (200)" +
        " and test (40) data");
      string trainFile =
        "..\\..\\..\\Data\\synthetic_train_200.txt";
      int[] colsX = new int[] { 0, 1, 2, 3, 4 };
      double[][] trainX =
        MatLoad(trainFile, colsX, ',', "#");
      double[] trainY =
        MatToVec(MatLoad(trainFile,
        new int[] { 5 }, ',', "#"));

      string testFile =
        "..\\..\\..\\Data\\synthetic_test_40.txt";
      double[][] testX =
        MatLoad(testFile, colsX, ',', "#");
      double[] testY =
        MatToVec(MatLoad(testFile,
        new int[] { 5 }, ',', "#"));
      Console.WriteLine("Done ");
      
      Console.WriteLine("\nFirst three train X: ");
      for (int i = 0; i "lt" 3; ++i)
        VecShow(trainX[i], 4, 8);

      Console.WriteLine("\nFirst three train y: ");
      for (int i = 0; i "lt" 3; ++i)
        Console.WriteLine(trainY[i].ToString("F4").
          PadLeft(8));

      //Random Forest with decision tree learners
      //Setting nTrees = 100
      //Setting maxDepth = 6
      //Setting minSamples = 2
      //Setting minLeaf = 1
      //Setting numSplitCols = 4
      //Setting nRows = 180
      //Accuracy train(within 0.10) = 0.8350
      //Accuracy test(within 0.10) = 0.6500
      //MSE train = 0.0004
      //MSE test = 0.0015

      int nTrees = 400;
      int maxDepth = 7;
      int minSamples = 2;
      int minLeaf = 1;
      // Standard heuristic for Random Forest Regression: 
      // p / 3 (where p is total feature columns = 5).
      //int numSplitCols = (int)(colsX.Length / 3.0); 
      //if (numSplitCols "lt" 1) numSplitCols = 1; // safety
      int numSplitCols = 4;
      
      int nRows = 180;  // train data for each tree

      Console.WriteLine("\nSetting nTrees = " + nTrees);
      Console.WriteLine("Setting maxDepth = " + maxDepth);
      Console.WriteLine("Setting minSamples = " +
        minSamples);
      Console.WriteLine("Setting minLeaf = " + minLeaf);
      Console.WriteLine("Setting numSplitCols = " +
        numSplitCols);
      Console.WriteLine("Setting nRows = " + nRows);

      Console.WriteLine("\nCreating and training" +
        " RandomForestRegressor model ");
      RandomForestRegressor rfr =
        new RandomForestRegressor(nTrees, nRows, maxDepth,
        minSamples, minLeaf, numSplitCols, seed: 0);
      rfr.Train(trainX, trainY);
      Console.WriteLine("Done ");

      // 3. evaluate model
      double accTrain = rfr.Accuracy(trainX, trainY, 0.10);
      Console.WriteLine("\nAccuracy train (within 0.10) = " +
        accTrain.ToString("F4"));
      double accTest = rfr.Accuracy(testX, testY, 0.10);
      Console.WriteLine("Accuracy test (within 0.10) = " +
        accTest.ToString("F4"));

      double mseTrain = rfr.MSE(trainX, trainY);
      Console.WriteLine("\nMSE train = " +
        mseTrain.ToString("F4"));
      double mseTest = rfr.MSE(testX, testY);
      Console.WriteLine("MSE test = " +
        mseTest.ToString("F4"));

      // 4. use model
      double[] x = trainX[0];
      Console.WriteLine("\nPredicting for x = ");
      VecShow(x, 4, 9);
      double predY = rfr.Predict(x);
      Console.WriteLine("Predicted y = " +
        predY.ToString("F4"));

      Console.WriteLine("\nEnd demo ");
      Console.ReadLine();
    } // Main()

    // ------------------------------------------------------
    // helpers for Main()
    // ------------------------------------------------------

    static double[][] MatLoad(string fn, int[] usecols,
      char sep, string comment)
    {
      List"lt"double[]"gt" result = 
        new List"lt"double[]"gt"();
      string line = "";
      FileStream ifs = new FileStream(fn, FileMode.Open);
      StreamReader sr = new StreamReader(ifs);
      while ((line = sr.ReadLine()) != null)
      {
        if (line.StartsWith(comment) == true)
          continue;
        string[] tokens = line.Split(sep);
        List"lt"double"gt" lst = 
          new List"lt"double"gt"();
        for (int j = 0; j "lt" usecols.Length; ++j)
          lst.Add(double.Parse(tokens[usecols[j]]));
        double[] row = lst.ToArray();
        result.Add(row);
      }
      sr.Close(); ifs.Close();
      return result.ToArray();
    }

    static double[] MatToVec(double[][] mat)
    {
      int nRows = mat.Length;
      int nCols = mat[0].Length;
      double[] result = new double[nRows * nCols];
      int k = 0;
      for (int i = 0; i "lt" nRows; ++i)
        for (int j = 0; j "lt" nCols; ++j)
          result[k++] = mat[i][j];
      return result;
    }

    static void VecShow(double[] vec, int dec, int wid)
    {
      for (int i = 0; i "lt" vec.Length; ++i)
        Console.Write(vec[i].ToString("F" + dec).
          PadLeft(wid));
      Console.WriteLine("");
    }

  } // class Program

  // ========================================================

  public class RandomForestRegressor // w/ ExtraTree
  {
    public int nTrees;
    public int maxDepth;
    public int minSamples;
    public int minLeaf;
    public int numSplitCols;
    public int nRows;
    public List"lt"ExtraTreeRegressor"gt" trees;
    public Random rnd;

    public RandomForestRegressor(int nTrees, int nRows,
      int maxDepth, int minSamples, int minLeaf,
      int numSplitCols, int seed = 0)
    {
      this.nTrees = nTrees;
      this.nRows = nRows;  // num train rows, each tree
      this.maxDepth = maxDepth;
      this.minSamples = minSamples;
      this.minLeaf = minLeaf;
      this.numSplitCols = numSplitCols;
      this.trees = new List"lt"ExtraTreeRegressor"gt"();
      this.rnd = new Random(seed);
    }

    public void Train(double[][] trainX, double[] trainY)
    {
      int totalRows = trainX.Length;
      int nCols = trainX[0].Length;

      // reusable buffers for all trees
      double[][] subsetX = new double[this.nRows][];
      for (int i = 0; i "lt" this.nRows; ++i)
        subsetX[i] = new double[nCols];
      
      double[] subsetY = new double[this.nRows];

      // train each tree using buffers
      for (int t = 0; t "lt" this.nTrees; ++t)
      {
        // populate the buffers
        for (int r = 0; r "lt" this.nRows; ++r)
        {
          int rndIndex = this.rnd.Next(0, totalRows);

          // Array.Copy() uses low-level native memory moves
          // much faster than a nested C# for-loop
          Array.Copy(trainX[rndIndex], subsetX[r], nCols);
          subsetY[r] = trainY[rndIndex];
        }

        ExtraTreeRegressor etr = new ExtraTreeRegressor(
          this.maxDepth,
          this.minSamples,
          this.minLeaf,
          this.numSplitCols, // Passed to sample columns
          false,
          this.rnd.Next(0, 1_000_000)
        );

        etr.Train(subsetX, subsetY); // safe to reuse
        this.trees.Add(etr);
      }
    }

    // ------------------------------------------------------

    public double Predict(double[] x)
    {
      double sum = 0.0;
      for (int t = 0; t "lt" this.nTrees; ++t)
        sum += this.trees[t].Predict(x);
      return sum / this.nTrees;
    }

    // ------------------------------------------------------

    public double Accuracy(double[][] dataX, double[] dataY,
      double pctClose)
    {
      int numCorrect = 0; int numWrong = 0;
      for (int i = 0; i "lt" dataX.Length; ++i)
      {
        double actualY = dataY[i];
        double predY = Predict(dataX[i]);
        if (Math.Abs(predY - actualY) "lt"
          (pctClose * Math.Abs(actualY)))
          ++numCorrect;
        else
          ++numWrong;
      }
      return (numCorrect * 1.0) / (numWrong + numCorrect);
    }

    // ------------------------------------------------------

    public double MSE(double[][] dataX, double[] dataY)
    {
      int n = dataX.Length;
      double sum = 0.0;
      for (int i = 0; i "lt" n; ++i)
      {
        double actualY = dataY[i];
        double predY = this.Predict(dataX[i]);
        sum += (actualY - predY) * (actualY - predY);
      }
      return sum / n;
    }

  } // class RandomForestRegressor

  // ========================================================

  
  public class ExtraTreeRegressor
  {
    // same as DecisionTreeRegressor except BestSplit()
    public int maxDepth;
    public int minSamples;  // aka min_samples_split
    public int minLeaf;  // min number of values in a leaf
    public int numSplitCols;
    public List"lt"Node"gt" tree;
    public Random rnd;  // order in which cols are searched
    public bool saveRows;  // keep rows in Nodes after train

    public double[][] trainX;  // store data by ref
    public double[] trainY;    // more convenient

    // ------------------------------------------------------

    public class Node
    {
      public int id;
      public int colIdx;      // aka feature index
      public double thresh;   // aka split value
      public int left;        // index into List
      public int right;
      public double value;    // aka predicted y
      public bool isLeaf;
      public List"lt"int"gt" rows;  // rows in train data

      public Node()
      {
        this.id = -1;
        this.colIdx = -1;
        this.thresh = 0.0;
        this.left = -1;
        this.right = -1;
        this.value = 0.0;  // aka pred y
        this.isLeaf = false;
        this.rows = null;
      }
    } // class Node

    // ------------------------------------------------------

    public ExtraTreeRegressor(int maxDepth = 2,
      int minSamples = 2, int minLeaf = 1,
      int numSplitCols = -1, bool saveRows = false,
      int seed = 0)
    {
      // if maxDepth = n, at most 2^(n+1) - 1 nodes
      this.maxDepth = maxDepth;
      this.minSamples = minSamples;
      this.minLeaf = minLeaf;
      this.numSplitCols = numSplitCols;  // -1 = all
      this.saveRows = saveRows;

      this.tree = new List"lt"Node"gt"();

      // create full tree List with null nodes
      int numNodes = (int)Math.Pow(2, (maxDepth + 1)) - 1;
      for (int i = 0; i "lt" numNodes; ++i)
        this.tree.Add(null);  // empty nodes

      this.rnd = new Random(seed);
    }

    // ------------------------------------------------------
    // public: ctor(), Train(), Predict()
    // private helpers: BestSplit(), TreeTargetMean()
    // ------------------------------------------------------

    public void Train(double[][] trainX, double[] trainY)
    {
      this.trainX = trainX; // useful to avoid passing args
      this.trainY = trainY;

      int maxID = (int)Math.Pow(2, (this.maxDepth + 1)) - 2;
      int maxStartID = (int)Math.Pow(2, this.maxDepth) - 1;

      // prepare root node
      List"lt"int"gt" allRows = new List"lt"int"gt"();
      for (int i = 0; i "lt" this.trainX.Length; ++i)
        allRows.Add(i);
      double grandMean = this.TreeTargetMean(allRows);

      // wait to supply root colIdx and thresh in main loop
      Node root = new Node();
      root.id = 0;
      root.left = 1;
      root.right = 2;
      root.value = grandMean;
      root.isLeaf = false; // (already set)
      root.rows = allRows;
      this.tree[0] = root;

      for (int i = 0; i "lt" this.tree.Count; ++i)
      {
        Node currNode = this.tree[i];
        if (currNode == null) continue;

        if (currNode.id "gte" maxStartID ||
          currNode.rows.Count "lt" this.minSamples)
        {
          // unable to attempt split
          currNode.isLeaf = true;
          continue;
        }

        // try to split curr node
        double[] splitInfo = this.BestSplit(currNode.rows);
        int colIdx = (int)splitInfo[0];
        double splitVal = splitInfo[1];  //split value

        if (colIdx == -1)  // bad split
        {
          currNode.isLeaf = true;
          currNode.left = -1;
          currNode.right = -1;
          continue;
        }

        // got successful split info
        // complete the fields for curr node
        currNode.colIdx = colIdx;
        currNode.thresh = splitVal;

        // compute associated rows for the children
        List"lt"int"gt" leftIdxs = 
          new List"lt"int"gt"();
        List"lt"int"gt" rightIdxs = 
          new List"lt"int"gt"();
        for (int k = 0; k "lt" currNode.rows.Count; ++k)
        {
          int r = currNode.rows[k];
          if (this.trainX[r][colIdx] "lte" splitVal)
            leftIdxs.Add(r);
          else
            rightIdxs.Add(r);
        }

        // explicitly assign child structural
        // pointers to the parent node
        int leftID = currNode.id * 2 + 1;
        if (leftID "lte" maxID && 
        leftIdxs.Count "gte" this.minLeaf)
        {
          currNode.left = leftID;
          Node leftNode = new Node();
          leftNode.id = leftID;
          leftNode.rows = leftIdxs;
          leftNode.value =
            this.TreeTargetMean(leftNode.rows);
          this.tree[leftID] = leftNode;
        }
        else
        {
          currNode.left = -1;
        }

        int rightID = currNode.id * 2 + 2;
        if (rightID "lte" maxID &&
          rightIdxs.Count "gte" this.minLeaf)
        {
          currNode.right = rightID;
          Node rightNode = new Node();
          rightNode.id = rightID;
          rightNode.rows = rightIdxs;
          rightNode.value =
            this.TreeTargetMean(rightNode.rows);
          this.tree[rightID] = rightNode;
        }
        else
        {
          currNode.right = -1;
        }

        // if child splits could not be formed
        // adequately, default parent to leaf
        if (currNode.left == -1 && currNode.right == -1)
        {
          currNode.isLeaf = true;
        }
      }

      if (this.saveRows == false)
      {
        for (int i = 0; i "lt" this.tree.Count; ++i)
          if (this.tree[i] != null)
            this.tree[i].rows = null;
      }
    } // Train()

    // ------------------------------------------------------

    public double Predict(double[] x)
    {
      int p = 0;
      double lastValidValue = 0.0;

      while (p != -1 && p "lt" this.tree.Count)
      {
        Node currNode = this.tree[p];
        if (currNode == null) break;

        lastValidValue = currNode.value;
        if (currNode.isLeaf == true) break;

        if (x[currNode.colIdx] "lte" currNode.thresh)
          p = currNode.left;
        else
          p = currNode.right;
      }
      return lastValidValue;
    }

    // ------------------------------------------------------
    // helpers: BestSplit(), TreeTargetMean()
    // ------------------------------------------------------

    private double[] BestSplit(List"lt"int"gt" rows)
    {
      // special algorithm for Extra Trees Regression
      // pick random col (from active cols)
      // pick random split val
      int bestColIdx = -1;
      double bestThresh = 0.0;
      double bestVar = double.MaxValue;

      int nRows = rows.Count;
      int nCols = this.trainX[0].Length;

      if (nRows == 0)
        throw new Exception("Empty data in BestSplit()");

      // 1. Fisher-Yates the columns
      int[] colIndices = new int[nCols];
      for (int k = 0; k "lt" nCols; ++k) colIndices[k] = k;

      for (int i = 0; i "lt" nCols - 1; ++i)
      {
        int ri = rnd.Next(i, nCols);
        int tmp = colIndices[i];
        colIndices[i] = colIndices[ri];
        colIndices[ri] = tmp;
      }

      int nColsToUse;
      if (this.numSplitCols != -1)
        nColsToUse = Math.Min(this.numSplitCols, nCols);
      else
        nColsToUse = nCols;

      int[] activeCols = new int[nColsToUse];
      for (int i = 0; i "lt" nColsToUse; ++i)
        activeCols[i] = colIndices[i];

      // pre-calculate total sums for the rows at this node
      double totalSumY = 0.0;
      double totalSumSqY = 0.0;
      for (int i = 0; i "lt" nRows; ++i)
      {
        int r = rows[i];
        double y = this.trainY[r];
        totalSumY += y;
        totalSumSqY += y * y;
      }

      // evaluate a random thresh in each active column
      for (int j = 0; j "lt" activeCols.Length; ++j)
      {
        int currColIdx = activeCols[j];

        // find min and max val in curr column
        // for the current rows
        double minVal = double.MaxValue;
        double maxVal = double.MinValue;
        for (int i = 0; i "lt" nRows; ++i)
        {
          int r = rows[i];
          double currVal = this.trainX[r][currColIdx];
          if (currVal "lt" minVal) minVal = currVal;
          if (currVal "gt" maxVal) maxVal = currVal;
        }

        if (minVal == maxVal)
          continue; // Cannot split this column

        double ranSplitVal = (maxVal - minVal) * 
          this.rnd.NextDouble() + minVal;

        int leftCount = 0;
        int rightCount = 0;
        double leftSumY = 0.0;
        double leftSumSqY = 0.0;

        for (int i = 0; i "lt" nRows; ++i)
        {
          int r = rows[i];
          double yCurr = this.trainY[r];
          if (this.trainX[r][currColIdx] "lte" ranSplitVal)
          {
            ++leftCount;
            leftSumY += yCurr;
            leftSumSqY += yCurr * yCurr;
          }
          else
          {
            ++rightCount;
          }
        }

        // enforce minLeaf
        if (leftCount "lt" this.minLeaf || 
          rightCount "lt" this.minLeaf)
          continue;

        double rightSumY = totalSumY - leftSumY;
        double rightSumSqY = totalSumSqY - leftSumSqY;

        double tmp1 = 
        (leftSumY / leftCount) * (leftSumY / leftCount);
        double leftVar = (leftSumSqY / leftCount) - tmp1;

        double tmp2 = 
        (rightSumY / rightCount) * (rightSumY / rightCount);
        double rightVar = (rightSumSqY / rightCount) - tmp2;

        if (leftVar "lt" 0.0) leftVar = 0.0;
        if (rightVar "lt" 0.0) rightVar = 0.0;

        double weightedVar = 
          ((leftCount * leftVar) + 
          (rightCount * rightVar)) / nRows;

        if (weightedVar "lt" bestVar)
        {
          bestVar = weightedVar;
          bestColIdx = currColIdx;
          bestThresh = ranSplitVal;
        }
      }

      double[] result = new double[2];
      result[0] = 1.0 * bestColIdx;
      result[1] = bestThresh;
      return result;
    }

    // ------------------------------------------------------

    private double TreeTargetMean(List"lt"int"gt" rows)
    {
      // mean of rows items in trainY
      // for node prediction
      if (rows == null || rows.Count == 0) return 0.0;
      double sum = 0.0;
      for (int i = 0; i "lt" rows.Count; ++i)
      {
        int r = rows[i];
        sum += this.trainY[r];
      }
      return sum / rows.Count;
    }

  } // class ExtraTreeRegressor

  // ========================================================

} // ns

Training data:

# synthetic_train_200.txt
#
-0.1660,  0.4406, -0.9998, -0.3953, -0.7065,  0.4840
 0.0776, -0.1616,  0.3704, -0.5911,  0.7562,  0.1568
-0.9452,  0.3409, -0.1654,  0.1174, -0.7192,  0.8054
 0.9365, -0.3732,  0.3846,  0.7528,  0.7892,  0.1345
-0.8299, -0.9219, -0.6603,  0.7563, -0.8033,  0.7955
 0.0663,  0.3838, -0.3690,  0.3730,  0.6693,  0.3206
-0.9634,  0.5003,  0.9777,  0.4963, -0.4391,  0.7377
-0.1042,  0.8172, -0.4128, -0.4244, -0.7399,  0.4801
-0.9613,  0.3577, -0.5767, -0.4689, -0.0169,  0.6861
-0.7065,  0.1786,  0.3995, -0.7953, -0.1719,  0.5569
 0.3888, -0.1716, -0.9001,  0.0718,  0.3276,  0.2500
 0.1731,  0.8068, -0.7251, -0.7214,  0.6148,  0.3297
-0.2046, -0.6693,  0.8550, -0.3045,  0.5016,  0.2129
 0.2473,  0.5019, -0.3022, -0.4601,  0.7918,  0.2613
-0.1438,  0.9297,  0.3269,  0.2434, -0.7705,  0.5171
 0.1568, -0.1837, -0.5259,  0.8068,  0.1474,  0.3307
-0.9943,  0.2343, -0.3467,  0.0541,  0.7719,  0.5581
 0.2467, -0.9684,  0.8589,  0.3818,  0.9946,  0.1092
-0.6553, -0.7257,  0.8652,  0.3936, -0.8680,  0.7018
 0.8460,  0.4230, -0.7515, -0.9602, -0.9476,  0.1996
-0.9434, -0.5076,  0.7201,  0.0777,  0.1056,  0.5664
 0.9392,  0.1221, -0.9627,  0.6013, -0.5341,  0.1533
 0.6142, -0.2243,  0.7271,  0.4942,  0.1125,  0.1661
 0.4260,  0.1194, -0.9749, -0.8561,  0.9346,  0.2230
 0.1362, -0.5934, -0.4953,  0.4877, -0.6091,  0.3810
 0.6937, -0.5203, -0.0125,  0.2399,  0.6580,  0.1460
-0.6864, -0.9628, -0.8600, -0.0273,  0.2127,  0.5387
 0.9772,  0.1595, -0.2397,  0.1019,  0.4907,  0.1611
 0.3385, -0.4702, -0.8673, -0.2598,  0.2594,  0.2270
-0.8669, -0.4794,  0.6095, -0.6131,  0.2789,  0.4700
 0.0493,  0.8496, -0.4734, -0.8681,  0.4701,  0.3516
 0.8639, -0.9721, -0.5313,  0.2336,  0.8980,  0.1412
 0.9004,  0.1133,  0.8312,  0.2831, -0.2200,  0.1782
 0.0991,  0.8524,  0.8375, -0.2102,  0.9265,  0.2150
-0.6521, -0.7473, -0.7298,  0.0113, -0.9570,  0.7422
 0.6190, -0.3105,  0.8802,  0.1640,  0.7577,  0.1056
 0.6895,  0.8108, -0.0802,  0.0927,  0.5972,  0.2214
 0.1982, -0.9689,  0.1870, -0.1326,  0.6147,  0.1310
-0.3695,  0.7858,  0.1557, -0.6320,  0.5759,  0.3773
-0.1596,  0.3581,  0.8372, -0.9992,  0.9535,  0.2071
-0.2468,  0.9476,  0.2094,  0.6577,  0.1494,  0.4132
 0.1737,  0.5000,  0.7166,  0.5102,  0.3961,  0.2611
 0.7290, -0.3546,  0.3416, -0.0983, -0.2358,  0.1332
-0.3652,  0.2438, -0.1395,  0.9476,  0.3556,  0.4170
-0.6029, -0.1466, -0.3133,  0.5953,  0.7600,  0.4334
-0.4596, -0.4953,  0.7098,  0.0554,  0.6043,  0.2775
 0.1450,  0.4663,  0.0380,  0.5418,  0.1377,  0.2931
-0.8636, -0.2442, -0.8407,  0.9656, -0.6368,  0.7429
 0.6237,  0.7499,  0.3768,  0.1390, -0.6781,  0.2185
-0.5499,  0.1850, -0.3755,  0.8326,  0.8193,  0.4399
-0.4858, -0.7782, -0.6141, -0.0008,  0.4572,  0.4197
 0.7033, -0.1683,  0.2334, -0.5327, -0.7961,  0.1776
 0.0317, -0.0457, -0.6947,  0.2436,  0.0880,  0.3345
 0.5031, -0.5559,  0.0387,  0.5706, -0.9553,  0.3107
-0.3513,  0.7458,  0.6894,  0.0769,  0.7332,  0.3170
 0.2205,  0.5992, -0.9309,  0.5405,  0.4635,  0.3532
-0.4806, -0.4859,  0.2646, -0.3094,  0.5932,  0.3202
 0.9809, -0.3995, -0.7140,  0.8026,  0.0831,  0.1600
 0.9495,  0.2732,  0.9878,  0.0921,  0.0529,  0.1289
-0.9476, -0.6792,  0.4913, -0.9392, -0.2669,  0.5966
 0.7247,  0.3854,  0.3819, -0.6227, -0.1162,  0.1550
-0.5922, -0.5045, -0.4757,  0.5003, -0.0860,  0.5863
-0.8861,  0.0170, -0.5761,  0.5972, -0.4053,  0.7301
 0.6877, -0.2380,  0.4997,  0.0223,  0.0819,  0.1404
 0.9189,  0.6079, -0.9354,  0.4188, -0.0700,  0.1907
-0.1428, -0.7820,  0.2676,  0.6059,  0.3936,  0.2790
 0.5324, -0.3151,  0.6917, -0.1425,  0.6480,  0.1071
-0.8432, -0.9633, -0.8666, -0.0828, -0.7733,  0.7784
-0.9444,  0.5097, -0.2103,  0.4939, -0.0952,  0.6787
-0.0520,  0.6063, -0.1952,  0.8094, -0.9259,  0.4836
 0.5477, -0.7487,  0.2370, -0.9793,  0.0773,  0.1241
 0.2450,  0.8116,  0.9799,  0.4222,  0.4636,  0.2355
 0.8186, -0.1983, -0.5003, -0.6531, -0.7611,  0.1511
-0.4714,  0.6382, -0.3788,  0.9648, -0.4667,  0.5950
 0.0673, -0.3711,  0.8215, -0.2669, -0.1328,  0.2677
-0.9381,  0.4338,  0.7820, -0.9454,  0.0441,  0.5518
-0.3480,  0.7190,  0.1170,  0.3805, -0.0943,  0.4724
-0.9813,  0.1535, -0.3771,  0.0345,  0.8328,  0.5438
-0.1471, -0.5052, -0.2574,  0.8637,  0.8737,  0.3042
-0.5454, -0.3712, -0.6505,  0.2142, -0.1728,  0.5783
 0.6327, -0.6297,  0.4038, -0.5193,  0.1484,  0.1153
-0.5424,  0.3282, -0.0055,  0.0380, -0.6506,  0.6613
 0.1414,  0.9935,  0.6337,  0.1887,  0.9520,  0.2540
-0.9351, -0.8128, -0.8693, -0.0965, -0.2491,  0.7353
 0.9507, -0.6640,  0.9456,  0.5349,  0.6485,  0.1059
-0.0462, -0.9737, -0.2940, -0.0159,  0.4602,  0.2606
-0.0627, -0.0852, -0.7247, -0.9782,  0.5166,  0.2977
 0.0478,  0.5098, -0.0723, -0.7504, -0.3750,  0.3335
 0.0090,  0.3477,  0.5403, -0.7393, -0.9542,  0.4415
-0.9748,  0.3449,  0.3736, -0.1015,  0.8296,  0.4358
 0.2887, -0.9895, -0.0311,  0.7186,  0.6608,  0.2057
 0.1570, -0.4518,  0.1211,  0.3435, -0.2951,  0.3244
 0.7117, -0.6099,  0.4946, -0.4208,  0.5476,  0.1096
-0.2929, -0.5726,  0.5346, -0.3827,  0.4665,  0.2465
 0.4889, -0.5572, -0.5718, -0.6021, -0.7150,  0.2163
-0.7782,  0.3491,  0.5996, -0.8389, -0.5366,  0.6516
-0.5847,  0.8347,  0.4226,  0.1078, -0.3910,  0.6134
 0.8469,  0.4121, -0.0439, -0.7476,  0.9521,  0.1571
-0.6803, -0.5948, -0.1376, -0.1916, -0.7065,  0.7156
 0.2878,  0.5086, -0.5785,  0.2019,  0.4979,  0.2980
 0.2764,  0.1943, -0.4090,  0.4632,  0.8906,  0.2960
-0.8877,  0.6705, -0.6155, -0.2098, -0.3998,  0.7107
-0.8398,  0.8093, -0.2597,  0.0614, -0.0118,  0.6502
-0.8476,  0.0158, -0.4769, -0.2859, -0.7839,  0.7715
 0.5751, -0.7868,  0.9714, -0.6457,  0.1448,  0.1175
 0.4802, -0.7001,  0.1022, -0.5668,  0.5184,  0.1090
 0.4458, -0.6469,  0.7239, -0.9604,  0.7205,  0.0779
 0.5175,  0.4339,  0.9747, -0.4438, -0.9924,  0.2879
 0.8678,  0.7158,  0.4577,  0.0334,  0.4139,  0.1678
 0.5406,  0.5012,  0.2264, -0.1963,  0.3946,  0.2088
-0.9938,  0.5498,  0.7928, -0.5214, -0.7585,  0.7687
 0.7661,  0.0863, -0.4266, -0.7233, -0.4197,  0.1466
 0.2277, -0.3517, -0.0853, -0.1118,  0.6563,  0.1767
 0.3499, -0.5570, -0.0655, -0.3705,  0.2537,  0.1632
 0.7547, -0.1046,  0.5689, -0.0861,  0.3125,  0.1257
 0.8186,  0.2110,  0.5335,  0.0094, -0.0039,  0.1391
 0.6858, -0.8644,  0.1465,  0.8855,  0.0357,  0.1845
-0.4967,  0.4015,  0.0805,  0.8977,  0.2487,  0.4663
 0.6760, -0.9841,  0.9787, -0.8446, -0.3557,  0.1509
-0.1203, -0.4885,  0.6054, -0.0443, -0.7313,  0.4854
 0.8557,  0.7919, -0.0169,  0.7134, -0.1628,  0.2002
 0.0115, -0.6209,  0.9300, -0.4116, -0.7931,  0.4052
-0.7114, -0.9718,  0.4319,  0.1290,  0.5892,  0.3661
 0.3915,  0.5557, -0.1870,  0.2955, -0.6404,  0.2954
-0.3564, -0.6548, -0.1827, -0.5172, -0.1862,  0.4622
 0.2392, -0.4959,  0.5857, -0.1341, -0.2850,  0.2470
-0.3394,  0.3947, -0.4627,  0.6166, -0.4094,  0.5325
 0.7107,  0.7768, -0.6312,  0.1707,  0.7964,  0.2757
-0.1078,  0.8437, -0.4420,  0.2177,  0.3649,  0.4028
-0.3139,  0.5595, -0.6505, -0.3161, -0.7108,  0.5546
 0.4335,  0.3986,  0.3770, -0.4932,  0.3847,  0.1810
-0.2562, -0.2894, -0.8847,  0.2633,  0.4146,  0.4036
 0.2272,  0.2966, -0.6601, -0.7011,  0.0284,  0.2778
-0.0743, -0.1421, -0.0054, -0.6770, -0.3151,  0.3597
-0.4762,  0.6891,  0.6007, -0.1467,  0.2140,  0.4266
-0.4061,  0.7193,  0.3432,  0.2669, -0.7505,  0.6147
-0.0588,  0.9731,  0.8966,  0.2902, -0.6966,  0.4955
-0.0627, -0.1439,  0.1985,  0.6999,  0.5022,  0.3077
 0.1587,  0.8494, -0.8705,  0.9827, -0.8940,  0.4263
-0.7850,  0.2473, -0.9040, -0.4308, -0.8779,  0.7199
 0.4070,  0.3369, -0.2428, -0.6236,  0.4940,  0.2215
-0.0242,  0.0513, -0.9430,  0.2885, -0.2987,  0.3947
-0.5416, -0.1322, -0.2351, -0.0604,  0.9590,  0.3683
 0.1055,  0.7783, -0.2901, -0.5090,  0.8220,  0.2984
-0.9129,  0.9015,  0.1128, -0.2473,  0.9901,  0.4776
-0.9378,  0.1424, -0.6391,  0.2619,  0.9618,  0.5368
 0.7498, -0.0963,  0.4169,  0.5549, -0.0103,  0.1614
-0.2612, -0.7156,  0.4538, -0.0460, -0.1022,  0.3717
 0.7720,  0.0552, -0.1818, -0.4622, -0.8560,  0.1685
-0.4177,  0.0070,  0.9319, -0.7812,  0.3461,  0.3052
-0.0001,  0.5542, -0.7128, -0.8336, -0.2016,  0.3803
 0.5356, -0.4194, -0.5662, -0.9666, -0.2027,  0.1776
-0.2378,  0.3187, -0.8582, -0.6948, -0.9668,  0.5474
-0.1947, -0.3579,  0.1158,  0.9869,  0.6690,  0.2992
 0.3992,  0.8365, -0.9205, -0.8593, -0.0520,  0.3154
-0.0209,  0.0793,  0.7905, -0.1067,  0.7541,  0.1864
-0.4928, -0.4524, -0.3433,  0.0951, -0.5597,  0.6261
-0.8118,  0.7404, -0.5263, -0.2280,  0.1431,  0.6349
 0.0516, -0.8480,  0.7483,  0.9023,  0.6250,  0.1959
-0.3212,  0.1093,  0.9488, -0.3766,  0.3376,  0.2735
-0.3481,  0.5490, -0.3484,  0.7797,  0.5034,  0.4379
-0.5785, -0.9170, -0.3563, -0.9258,  0.3877,  0.4121
 0.3407, -0.1391,  0.5356,  0.0720, -0.9203,  0.3458
-0.3287, -0.8954,  0.2102,  0.0241,  0.2349,  0.3247
-0.1353,  0.6954, -0.0919, -0.9692,  0.7461,  0.3338
 0.9036, -0.8982, -0.5299, -0.8733, -0.1567,  0.1187
 0.7277, -0.8368, -0.0538, -0.7489,  0.5458,  0.0830
 0.9049,  0.8878,  0.2279,  0.9470, -0.3103,  0.2194
 0.7957, -0.1308, -0.5284,  0.8817,  0.3684,  0.2172
 0.4647, -0.4931,  0.2010,  0.6292, -0.8918,  0.3371
-0.7390,  0.6849,  0.2367,  0.0626, -0.5034,  0.7039
-0.1567, -0.8711,  0.7940, -0.5932,  0.6525,  0.1710
 0.7635, -0.0265,  0.1969,  0.0545,  0.2496,  0.1445
 0.7675,  0.1354, -0.7698, -0.5460,  0.1920,  0.1728
-0.5211, -0.7372, -0.6763,  0.6897,  0.2044,  0.5217
 0.1913,  0.1980,  0.2314, -0.8816,  0.5006,  0.1998
 0.8964,  0.0694, -0.6149,  0.5059, -0.9854,  0.1825
 0.1767,  0.7104,  0.2093,  0.6452,  0.7590,  0.2832
-0.3580, -0.7541,  0.4426, -0.1193, -0.7465,  0.5657
-0.5996,  0.5766, -0.9758, -0.3933, -0.9572,  0.6800
 0.9950,  0.1641, -0.4132,  0.8579,  0.0142,  0.2003
-0.4717, -0.3894, -0.2567, -0.5111,  0.1691,  0.4266
 0.3917, -0.8561,  0.9422,  0.5061,  0.6123,  0.1212
-0.0366, -0.1087,  0.3449, -0.1025,  0.4086,  0.2475
 0.3633,  0.3943,  0.2372, -0.6980,  0.5216,  0.1925
-0.5325, -0.6466, -0.2178, -0.3589,  0.6310,  0.3568
 0.2271,  0.5200, -0.1447, -0.8011, -0.7699,  0.3128
 0.6415,  0.1993,  0.3777, -0.0178, -0.8237,  0.2181
-0.5298, -0.0768, -0.6028, -0.9490,  0.4588,  0.4356
 0.6870, -0.1431,  0.7294,  0.3141,  0.1621,  0.1632
-0.5985,  0.0591,  0.7889, -0.3900,  0.7419,  0.2945
 0.3661,  0.7984, -0.8486,  0.7572, -0.6183,  0.3449
 0.6995,  0.3342, -0.3113, -0.6972,  0.2707,  0.1712
 0.2565,  0.9126,  0.1798, -0.6043, -0.1413,  0.2893
-0.3265,  0.9839, -0.2395,  0.9854,  0.0376,  0.4770
 0.2690, -0.1722,  0.9818,  0.8599, -0.7015,  0.3954
-0.2102, -0.0768,  0.1219,  0.5607, -0.0256,  0.3949
 0.8216, -0.9555,  0.6422, -0.6231,  0.3715,  0.0801
-0.2896,  0.9484, -0.7545, -0.6249,  0.7789,  0.4370
-0.9985, -0.5448, -0.7092, -0.5931,  0.7926,  0.5402

Test data:

# synthetic_test_40.txt
#
 0.7462,  0.4006, -0.0590,  0.6543, -0.0083,  0.1935
 0.8495, -0.2260, -0.0142, -0.4911,  0.7699,  0.1078
-0.2335, -0.4049,  0.4352, -0.6183, -0.7636,  0.5088
 0.1810, -0.5142,  0.2465,  0.2767, -0.3449,  0.3136
-0.8650,  0.7611, -0.0801,  0.5277, -0.4922,  0.7140
-0.2358, -0.7466, -0.5115, -0.8413, -0.3943,  0.4533
 0.4834,  0.2300,  0.3448, -0.9832,  0.3568,  0.1360
-0.6502, -0.6300,  0.6885,  0.9652,  0.8275,  0.3046
-0.3053,  0.5604,  0.0929,  0.6329, -0.0325,  0.4756
-0.7995,  0.0740, -0.2680,  0.2086,  0.9176,  0.4565
-0.2144, -0.2141,  0.5813,  0.2902, -0.2122,  0.4119
-0.7278, -0.0987, -0.3312, -0.5641,  0.8515,  0.4438
 0.3793,  0.1976,  0.4933,  0.0839,  0.4011,  0.1905
-0.8568,  0.9573, -0.5272,  0.3212, -0.8207,  0.7415
-0.5785,  0.0056, -0.7901, -0.2223,  0.0760,  0.5551
 0.0735, -0.2188,  0.3925,  0.3570,  0.3746,  0.2191
 0.1230, -0.2838,  0.2262,  0.8715,  0.1938,  0.2878
 0.4792, -0.9248,  0.5295,  0.0366, -0.9894,  0.3149
-0.4456,  0.0697,  0.5359, -0.8938,  0.0981,  0.3879
 0.8629, -0.8505, -0.4464,  0.8385,  0.5300,  0.1769
 0.1995,  0.6659,  0.7921,  0.9454,  0.9970,  0.2330
-0.0249, -0.3066, -0.2927, -0.4923,  0.8220,  0.2437
 0.4513, -0.9481, -0.0770, -0.4374, -0.9421,  0.2879
-0.3405,  0.5931, -0.3507, -0.3842,  0.8562,  0.3987
 0.9538,  0.0471,  0.9039,  0.7760,  0.0361,  0.1706
-0.0887,  0.2104,  0.9808,  0.5478, -0.3314,  0.4128
-0.8220, -0.6302,  0.0537, -0.1658,  0.6013,  0.4306
-0.4123, -0.2880,  0.9074, -0.0461, -0.4435,  0.5144
 0.0060,  0.2867, -0.7775,  0.5161,  0.7039,  0.3599
-0.7968, -0.5484,  0.9426, -0.4308,  0.8148,  0.2979
 0.7811,  0.8450, -0.6877,  0.7594,  0.2640,  0.2362
-0.6802, -0.1113, -0.8325, -0.6694, -0.6056,  0.6544
 0.3821,  0.1476,  0.7466, -0.5107,  0.2592,  0.1648
 0.7265,  0.9683, -0.9803, -0.4943, -0.5523,  0.2454
-0.9049, -0.9797, -0.0196, -0.9090, -0.4433,  0.6447
-0.4607,  0.1811, -0.2389,  0.4050, -0.0078,  0.5229
 0.2664, -0.2932, -0.4259, -0.7336,  0.8742,  0.1834
-0.4507,  0.1029, -0.6294, -0.1158, -0.6294,  0.6081
 0.8948, -0.0124,  0.9278,  0.2899, -0.0314,  0.1534
-0.1323, -0.8813, -0.0146, -0.0697,  0.6135,  0.2386

Bottom line: I compared my from-scratch C# implementation of gradient boosting regression on training data with categorical predictor variables with the scikit-learn GradientBoostingRegression library module and got essentially the same results.

The goal of a machine learning regression problem is to predict a single numeric value. For example, you might want to predict the Income of a person based on their Sex, Age, City, Height, and Psychology (timid, alert, brave). There are many different regression techniques. One of the most common techniques is gradient boosting regression.

Gradient boosting regression (GBR) is a tree-based technique. GBR creates a collection of small decision trees (typically about 100) where each tree learns how to predict training items that were poorly predicted by the previous tree.

I wanted to validate my from-scratch C# implementation of GBR against the Python language scikit-learn library GradientBoostingRegressor module. I generated a set of synthetic training data that looks like:

M, 29, Austin, 68.5, timid, 54600.00
F, 51, Juneau, 60.0, brave, 25700.00
M, 32 ,Boston, 73.5, alert, 74900.00
. . .

The goal is to predict Income (in the last column) from Sex, Age, City, Height (in inches), and Psychology (timid, alert, brave). There are 200 training items and 40 test items.

The Sex variable is encoded as M = 0, F = 1. There are 10 possible City values: Austin, Boston, Casper, Denver, Eugene, Fresno, Gallup, Helena, Irvine, Juneau.

I used sorted-target-mean ordinal encoding, where the city that has the smallest average Income (Boston = $30,600.00) is encoded as 0, the next larger average Income city (Casper = $33,100.00) is encoded as 1, and so on to the largest average income city (Helena = $64,194.74) is 9.

I used one-hot encoding for Psychology: timid = 1,0,0 alert = 0,1,0 brave = 0,0,1.

When using tree-based regression, it’s not necessary to normalize numeric predictor variables, but it doesn’t hurt and can be helpful in some scenarios. I divided all Age values by 100, and all Height values by 100. I normalized the target Income variable by dividing all values by 100,000.

The resulting normalized and encoded data looks like:

0, 0.29, 2, 0.685, 1,0,0, 0.54600
1, 0.51, 8, 0.600, 0,0,1, 0.25700
0, 0.32, 0, 0.735, 0,1,0, 0.74900
. . . 

The key parts of the output of my C# implementation of gradient boosting regression are:

Begin C# Gradient Boosting regression demo
Predict Income from Sex, Age, City Height, Psychology

Loading people train (200) and test (40) data
Done

First three raw train X:
M,29,Austin,68.5,timid,54600.00
F,51,Juneau,60.0,brave,25700.00
M,32,Boston,73.5,alert,74900.00

First three train X:
  0.00000  0.29000  2.00000  0.68500  1.00000  0.00000  0.00000
  1.00000  0.51000  8.00000  0.60000  0.00000  0.00000  1.00000
  0.00000  0.32000  0.00000  0.73500  0.00000  1.00000  0.00000

First three train y:
  0.54600
  0.25700
  0.74900

Setting numTrees = 200
Setting maxDepth = 4
Setting minSamples = 2
Setting minLeaf = 1
Setting lrnRate = 0.0500

Creating and training GradientBoostRegression model
Done

Evaluating model

Accuracy train (within 0.10) = 1.0000
Accuracy test (within 0.10) = 0.7750

MSE train = 0.0000
MSE test = 0.0022

R2 train = 0.9993
R2 test = 0.9344

Predicting for x =
   0.0000   0.2900   2.0000   0.6850   1.0000   0.0000   0.0000
Predicted y = 0.5391

Explaining for x =
   0.0000   0.2900   2.0000   0.6850   1.0000   0.0000   0.0000

Initial prediction: 0.4685
t =   0  pred_res =  -0.0342  delta =  -0.0017  pred =  0.4668
t =   1  pred_res =   0.0079  delta =   0.0004  pred =  0.4672
t =   2  pred_res =  -0.0315  delta =  -0.0016  pred =  0.4656
t =   3  pred_res =  -0.0361  delta =  -0.0018  pred =  0.4638
t =   4  pred_res =  -0.0136  delta =  -0.0007  pred =  0.4632
. . .
t = 198  pred_res =  -0.0000  delta =  -0.0000  pred =  0.5391
t = 199  pred_res =   0.0000  delta =   0.0000  pred =  0.5391
Predicted y = 0.5391

End demo

I ran the data through the scikit GradientBoostingRegression module and got essentially the same results:

Begin scikit gradient boost regression demo
Predict Income from Sex, Age, City, Height, Psych

Setting n_estimators = 200
Setting max_depth = 4
Setting min_samples_split = 2
Setting min_samples_leaf = 1
Setting lrn_rate = 0.0500

Training gradient boost regression model
Done

Evaluating model

Acc train (within 0.10): 1.0000
Acc test (within 0.10): 0.8000

MSE train: 0.0000
MSE test: 0.0022

R2 train: 0.9993
R2 test: 0.9347

Predicting for x =
[0.0000 0.2900 2.0000 0.6850 1.0000 0.0000 0.0000]
Predicted y = 0.5391

End demo

A big advantage of my from-scratch C# implementation of GBR is that it was easy for me to add statements that show exactly how a prediction is made.

Machine learning regression models are often beautiful to my eye. But I’m not a good judge of other kinds of models. I’ve noticed that fashion models often look like they walk with an exaggerated limp. Here are three images that were returned by an Internet search for “ugly fashion model on runway”.

Right: OK, there isn’t much to disagree about “ugly fashion model on runway” here. I don’t want to sound mean, but she doesn’t belong on a fashion runway. An airport runway maybe, but not a fashion runway.

Demo C# program. Very long and complex. Replace “lt” (less than), “gt”, “lte”, “gte” with Boolean operator symbols (my lame blog editor consistently chokes on symbols).

using System;
using System.IO;
using System.Collections.Generic;

namespace GradientBoostRegression
{
  internal class GradientBoostRegressionProgram
  {
    static void Main(string[] args)
    {
      Console.WriteLine("\nBegin C# Gradient Boosting" +
        " regression demo ");

      Console.WriteLine("Predict Income from Sex, Age," +
        " City Height, Psychology ");

      // 1. load data
      Console.WriteLine("\nLoading people train" +
        " (200) and test (40) data");
      string trainFile =
        "..\\..\\..\\Data\\People\\people_train_200.txt";
      int[] colsX = new int[] { 0, 1, 2, 3, 4,5,6 };
      double[][] trainX =
        MatLoad(trainFile, colsX, ',', "#");
      double[] trainY =
        MatToVec(MatLoad(trainFile,
        [7], ',', "#"));

      string testFile =
        "..\\..\\..\\Data\\People\\people_test_40.txt";
      double[][] testX =
        MatLoad(testFile, colsX, ',', "#");
      double[] testY =
        MatToVec(MatLoad(testFile,
        [7], ',', "#"));
      Console.WriteLine("Done ");

      Console.WriteLine("\nFirst three raw train X: ");
      Console.WriteLine("M,29,Austin,68.5,timid,54600.00 ");
      Console.WriteLine("F,51,Juneau,60.0,brave,25700.00 ");
      Console.WriteLine("M,32,Boston,73.5,alert,74900.00 ");

      Console.WriteLine("\nFirst three train X: ");
      for (int i = 0; i "lt" 3; ++i)
        VecShow(trainX[i], 4, 8);

      Console.WriteLine("\nFirst three train y: ");
      for (int i = 0; i "lt" 3; ++i)
        Console.WriteLine(trainY[i].ToString("F4").
          PadLeft(8));

      // 2. create and train model
      int numTrees = 200;
      int maxDepth = 4;
      int minSamples = 2;
      int minLeaf = 1;
      int numSplitCols = -1; // use all columns
      double lrnRate = 0.10;


      Console.WriteLine("\nSetting numTrees = " +
        numTrees);
      Console.WriteLine("Setting maxDepth = " +
        maxDepth);
      Console.WriteLine("Setting minSamples = " +
        minSamples);
      Console.WriteLine("Setting minLeaf = " +
        minLeaf);
      Console.WriteLine("Setting lrnRate = " +
        lrnRate.ToString("F4"));

      Console.WriteLine("\nCreating and training" +
        " GradientBoostRegression model ");
      GradientBoostRegressor gbr =
        new GradientBoostRegressor(numTrees, maxDepth,
        minSamples, minLeaf, numSplitCols, lrnRate);
      gbr.Train(trainX, trainY);
      Console.WriteLine("Done ");

      // 3. evaluate model
      Console.WriteLine("\nEvaluating model ");
      double accTrain = gbr.Accuracy(trainX, trainY, 0.10);
      Console.WriteLine("\nAccuracy train (within 0.10) = " +
        accTrain.ToString("F4"));
      double accTest = gbr.Accuracy(testX, testY, 0.10);
      Console.WriteLine("Accuracy test (within 0.10) = " +
        accTest.ToString("F4"));

      double mseTrain = gbr.MSE(trainX, trainY);
      Console.WriteLine("\nMSE train = " +
        mseTrain.ToString("F4"));
      double mseTest = gbr.MSE(testX, testY);
      Console.WriteLine("MSE test = " +
        mseTest.ToString("F4"));

      double r2Train = gbr.R2(trainX, trainY);
      Console.WriteLine("\nR2 train = " +
        r2Train.ToString("F4"));
      double r2Test = gbr.R2(testX, testY);
      Console.WriteLine("R2 test = " +
        r2Test.ToString("F4"));

      // 4. use model
      double[] x = trainX[0];
      Console.WriteLine("\nPredicting for x = ");
      VecShow(x, 4, 9);
      double predY = gbr.Predict(x, verbose: true);
      Console.WriteLine("Predicted y = " +
        predY.ToString("F4"));

      Console.WriteLine("\nEnd demo ");
      Console.ReadLine();
    } // Main()

    // ------------------------------------------------------
    // helpers for Main()
    // ------------------------------------------------------

    static double[][] MatLoad(string fn, int[] usecols,
      char sep, string comment)
    {
      List"lt"double[]"gt" result = 
        new List"lt"double[]"gt"();
      string line = "";
      FileStream ifs = new FileStream(fn, FileMode.Open);
      StreamReader sr = new StreamReader(ifs);
      while ((line = sr.ReadLine()) != null)
      {
        if (line.StartsWith(comment) == true)
          continue;
        string[] tokens = line.Split(sep);
        List"lt"double"gt" lst = new List"lt"double"gt"();
        for (int j = 0; j "lt" usecols.Length; ++j)
          lst.Add(double.Parse(tokens[usecols[j]]));
        double[] row = lst.ToArray();
        result.Add(row);
      }
      sr.Close(); ifs.Close();
      return result.ToArray();
    }

    static double[] MatToVec(double[][] mat)
    {
      int nRows = mat.Length;
      int nCols = mat[0].Length;
      double[] result = new double[nRows * nCols];
      int k = 0;
      for (int i = 0; i "lt" nRows; ++i)
        for (int j = 0; j "lt" nCols; ++j)
          result[k++] = mat[i][j];
      return result;
    }

    static void VecShow(double[] vec, int dec, int wid)
    {
      for (int i = 0; i "lt" vec.Length; ++i)
        Console.Write(vec[i].ToString("F" + dec).
          PadLeft(wid));
      Console.WriteLine("");
    }

  } // class Program

  // ========================================================

  public class GradientBoostRegressor
  {
    public double lrnRate;
    public int nTrees;
    public int maxDepth;
    public int minSamples;
    public int minLeaf;
    public int numSplitCols;
    public List"lt"DecisionTreeRegressor"gt" trees;
    public double pred0;  // initial prediction
    private Random rnd;

    public GradientBoostRegressor(int nTrees, int maxDepth,
      int minSamples, int minLeaf, int numSplitCols,
      double lrnRate, int seed=0)
    {
      this.nTrees = nTrees;
      this.maxDepth = maxDepth;
      this.minSamples = minSamples;
      this.minLeaf = minLeaf;
      this.numSplitCols = numSplitCols;
      this.trees = new List"lt"DecisionTreeRegressor"gt"();
      this.lrnRate = lrnRate;
      this.rnd = new Random(seed);
    }

    public void Train(double[][] trainX, double[] trainY)
    {
      int n = trainX.Length;
      this.pred0 = Mean(trainY);

      double[] preds = new double[n]; //each data item
      for (int i = 0; i "lt" n; ++i)
        preds[i] = this.pred0;

      double[] residuals = new double[n]; // for curr tree
      for (int t = 0; t "lt" this.nTrees; ++t) // each tree
      {
        // double[] residuals = new double[n]; // curr tree
        for (int i = 0; i "lt" n; ++i)
          residuals[i] = trainY[i] - preds[i];

        DecisionTreeRegressor dtr =
          new DecisionTreeRegressor(this.maxDepth,
          this.minSamples, this.minLeaf, this.numSplitCols,
          false, this.rnd.Next(0,1_000_000));
        dtr.Train(trainX, residuals); // predict residuals

        for (int i = 0; i "lt" n; ++i)
        {
          double predResidual = dtr.Predict(trainX[i]);
          preds[i] += this.lrnRate * predResidual;
        }
        this.trees.Add(dtr);
      }
    } // Train

    public double Predict(double[] x, bool verbose = false)
    {
      double result = this.pred0;
      if (verbose == true)
      {
        Console.WriteLine("\nInitial prediction: " +
        result.ToString("F4"));
      }
      for (int t = 0; t "lt" this.nTrees; ++t)
      {
        double predResidual = this.trees[t].Predict(x);
        double delta = this.lrnRate * predResidual;
        result += delta;

        if (verbose == true)
        {
          if (t "gte" 0 && t "lte" 4 ||
            t "gte" this.nTrees - 2 &&
            t "lte" this.nTrees - 1)
          {
            Console.Write("t = " + t.ToString().PadLeft(3) +
              "  pred_res = " + predResidual.ToString("F4").
              PadLeft(10));
            Console.Write("  delta = " + delta.ToString("F4").
              PadLeft(10));
            Console.WriteLine("  pred = " +
              result.ToString("F4").PadLeft(7));
          }
          if (t == 5)
            Console.WriteLine(". . . ");
        }
      }
      return result;
    }

    public double Accuracy(double[][] dataX, double[] dataY,
      double pctClose)
    {
      int numCorrect = 0; int numWrong = 0;
      for (int i = 0; i "lt" dataX.Length; ++i)
      {
        double actualY = dataY[i];
        double predY = Predict(dataX[i]);
        if (Math.Abs(predY - actualY) "lt"
          (pctClose * Math.Abs(actualY)))
          ++numCorrect;
        else
          ++numWrong;
      }
      return (numCorrect * 1.0) / (numWrong + numCorrect);
    }

    public double MSE(double[][] dataX,
      double[] dataY)
    {
      int n = dataX.Length;
      double sum = 0.0;
      for (int i = 0; i "lt" n; ++i)
      {
        double actualY = dataY[i];
        double predY = this.Predict(dataX[i]);
        sum += (actualY - predY) * (actualY - predY);
      }
      return sum / n;
    }

    public double R2(double[][] dataX, double[] dataY)
    {
      // coefficient of determination
      int n = dataX.Length;
      double sum = 0.0;

      for (int i = 0; i "lt" n; ++i)
        sum += dataY[i];
      double meanActual = sum / n;

      double ssRes = 0.0;
      double ssTot = 0.0;
      for (int i = 0; i "lt" n; ++i)
      {
        double predY = this.Predict(dataX[i]);
        ssRes += (dataY[i] - predY) * (dataY[i] - predY);
        ssTot += (dataY[i] - meanActual) *
          (dataY[i] - meanActual);
      }
      if (Math.Abs(ssTot) "lt" 1.0e-12)
        return 0.0;
      else
        return 1.0 - (ssRes / ssTot);
    }

    private static double Mean(double[] data)
    {
      int n = data.Length;
      double sum = 0.0;
      for (int i = 0; i "lt" n; ++i)
        sum += data[i];
      return sum / n;
    }

  } // class GradientBoostRegressor

  // ========================================================

  public class DecisionTreeRegressor
  {
    // ChatGPT: "I don't see any fundamental tree-building or
    // indexing bugs. If the inputs are guaranteed to be
    // finite doubles and the parameters are validated 
    // elsewhere, the core regression algorithm is correct."

    public int maxDepth;
    public int minSamples;
    public int minLeaf;
    public int numSplitCols;
    public List"lt"Node"gt" tree;
    public Random rnd;
    public bool saveRows;
    public double[][] trainX;
    public double[] trainY;

    // ............................................

    public class Node
    {
      public int id;
      public int colIdx = -1;
      public double thresh = 0.0;
      public int left = -1;
      public int right = -1;
      public double value = 0.0;
      public bool isLeaf = false;
      public List"lt"int"gt" rows;
    }

    // ............................................

    public DecisionTreeRegressor(int maxDepth = 3,
      int minSamples = 2, int minLeaf = 1,
      int numSplitCols = -1, bool saveRows = false,
      int seed = 0)
    {
      this.maxDepth = maxDepth;
      this.minSamples = minSamples;
      this.minLeaf = minLeaf;
      this.numSplitCols = numSplitCols;
      this.saveRows = saveRows;
      this.tree = new List"lt"Node"gt"();

      int numNodes = (int)Math.Pow(2, (maxDepth + 1)) - 1;
      for (int i = 0; i "lt" numNodes; ++i)
        this.tree.Add(null);
      this.rnd = new Random(seed);
    }

    // ------------------------------------------------------
    // ------------------------------------------------------
    // public: ctor(), Train(), Predict()
    // private: BestSplit(), TreeTargetMean()
    // ------------------------------------------------------

    public void Train(double[][] trainX, double[] trainY)
    {
      this.trainX = trainX;
      this.trainY = trainY;

      int maxParentID = (int)Math.Pow(2, this.maxDepth) - 2;

      List"lt"int"gt" allRows =
        new List"lt"int"gt"(this.trainX.Length);
      for (int i = 0; i "lt" this.trainX.Length; ++i)
        allRows.Add(i);

      // construct root node
      Node root = new Node()
      {
        id = 0,
        value = this.TreeTargetMean(allRows),
        rows = allRows
      };
      this.tree[0] = root;

      for (int i = 0; i "lt" this.tree.Count; ++i)
      {
        Node currNode = this.tree[i];
        if (currNode == null) continue;

        int leftID = currNode.id * 2 + 1;
        int rightID = currNode.id * 2 + 2;

        // 1. check before expensive split calcs
        if (currNode.id "gt" maxParentID ||
            currNode.rows.Count "lt" this.minSamples ||
            rightID "gte" this.tree.Count)
        {
          currNode.isLeaf = true;
          continue;
        }

        // 2. find best split
        double[] splitInfo = this.BestSplit(currNode.rows);
        int colIdx = (int)splitInfo[0];
        double splitVal = splitInfo[1];

        if (colIdx == -1) // unable to split
        {
          currNode.isLeaf = true;
          continue;
        }

        // 3. apply split properties
        currNode.colIdx = colIdx;  // aka split col
        currNode.thresh = splitVal;

        List"lt"int"gt" leftIdxs = new List"lt"int"gt"();
        List"lt"int"gt" rightIdxs = new List"lt"int"gt"();

        for (int k = 0; k "lt" currNode.rows.Count; ++k)
        {
          int r = currNode.rows[k];
          if (this.trainX[r][colIdx] "lte" splitVal)
            leftIdxs.Add(r);
          else
            rightIdxs.Add(r);
        }

        // 4. create and attach children
        currNode.left = leftID;
        this.tree[leftID] = new Node()
        {
          id = leftID,
          rows = leftIdxs,
          value = this.TreeTargetMean(leftIdxs)
        };

        currNode.right = rightID;
        this.tree[rightID] = new Node()
        {
          id = rightID,
          rows = rightIdxs,
          value = this.TreeTargetMean(rightIdxs)
        };
      } // main loop

      // delete row info (for ensembles)
      if (this.saveRows == false)
      {
        for (int i = 0; i "lt" this.tree.Count; ++i)
          if (this.tree[i] != null)
            this.tree[i].rows = null;
      }

    } // Train()

    // ------------------------------------------------------

    public double Predict(double[] x)
    {
      int p = 0;
      double lastValidValue = 0.0;
      while (p != -1 && p "lt" this.tree.Count)
      {
        Node currNode = this.tree[p];
        if (currNode == null) break;

        lastValidValue = currNode.value;
        if (currNode.isLeaf || 
          currNode.colIdx == -1) break;

        if (x[currNode.colIdx] "lte" currNode.thresh)
          p = currNode.left;
        else
          p = currNode.right;
      }
      return lastValidValue;
    }

    // ------------------------------------------------------

    private double[] BestSplit(List"lt"int"gt" rows)
    {
      int bestColIdx = -1;
      double bestThresh = 0.0;
      double bestGain = -1.0;
      int nRows = rows.Count;
      int nCols = this.trainX[0].Length;

      int[] colIndices = new int[nCols];
      for (int k = 0; k "lt" nCols; ++k)
        colIndices[k] = k;

      // feature sampling
      for (int i = 0; i "lt" nCols; ++i)
      {
        int ri = rnd.Next(i, nCols);
        int tmp = colIndices[i];
        colIndices[i] = colIndices[ri];
        colIndices[ri] = tmp;
      }

      int numColsToUse;
      if (this.numSplitCols == -1)
        numColsToUse = nCols;
      else
        numColsToUse = Math.Min(this.numSplitCols, nCols);

      int[] sortedRows = new int[nRows];
      double[] featureKeys = new double[nRows];

      double totalSum = 0.0;
      for (int i = 0; i "lt" nRows; ++i)
        totalSum += this.trainY[rows[i]];

      for (int j = 0; j "lt" numColsToUse; ++j)
      {
        int colIdx = colIndices[j];
        for (int i = 0; i "lt" nRows; ++i)
        {
          int r = rows[i];
          sortedRows[i] = r;
          featureKeys[i] = this.trainX[r][colIdx];
        }

        Array.Sort(featureKeys, sortedRows);
        double leftSum = 0.0;

        for (int i = 0; i "lt" nRows - 1; ++i)
        {
          leftSum += this.trainY[sortedRows[i]];

          int leftCount = i + 1;
          int rightCount = nRows - leftCount;

          if (featureKeys[i] == featureKeys[i + 1])
            continue; // duplicate values
          if (leftCount "lt" this.minLeaf ||
            rightCount "lt" this.minLeaf)
            continue;

          double rightSum = totalSum - leftSum;

          // maximizing Gain = (G_L^2 / N_L) + (G_R^2 / N_R)
          // equivalent to maximizing explicit weighted
          // variance reduction
          double gain = (leftSum * leftSum / leftCount) +
            (rightSum * rightSum / rightCount);

          if (gain "gt" bestGain)
          {
            bestGain = gain;
            bestColIdx = colIdx;
            // bestThresh = 
            //  (featureKeys[i] + featureKeys[i + 1]) / 2.0;
            bestThresh = featureKeys[i] +
              (featureKeys[i + 1] - featureKeys[i]) / 2.0;
          }
        }
      }

      // using a tuple would create version dependency
      double[] result = new double[2];
      result[0] = 1.0 * bestColIdx;
      result[1] = bestThresh;
      return result;
    }

    private double TreeTargetMean(List"lt"int"gt" rows)
    {
      if (rows == null || rows.Count == 0)
        return 0.0;
      double sum = 0.0;
      for (int i = 0; i "lt" rows.Count; ++i)
        sum += this.trainY[rows[i]];
      return sum / rows.Count;
    }

  } // DecisionTreeRegressor

  // ========================================================

  

Demo scikit program.

# scikit_gradient_boost_regression.py
# on custom People dataset

import numpy as np
from sklearn.ensemble import GradientBoostingRegressor

# -----------------------------------------------------------

def accuracy(model, data_X, data_y, pct_close):
  # assumes model has a predict(X)
  n = len(data_X)
  n_correct = 0; n_wrong = 0
  for i in range(n):
    x = data_X[i].reshape(1,-1)  # make it a matrix
    y = data_y[i]
    y_pred = model.predict(x)[0]  # predict() expects 2D

    if np.abs(y - y_pred) "lt" np.abs(y) * pct_close:
      n_correct += 1
    else: 
      n_wrong += 1
  # print("Correct = " + str(n_correct))
  # print("Wrong   = " + str(n_wrong))
  return n_correct / (n_correct + n_wrong)

# -----------------------------------------------------------

def MSE(model, data_X, data_y):
  n = len(data_X)
  sum = 0.0
  for i in range(n):
    x = data_X[i].reshape(1,-1)
    y = data_y[i]
    y_pred = model.predict(x)[0]
    sum += (y - y_pred) * (y - y_pred)

  return sum / n

# -----------------------------------------------------------

def main():
  print("\nBegin scikit gradient boost regression demo ")
  print("Predict Income from Sex, Age, City, Height, Psych ")

  np.set_printoptions(precision=4, suppress=True,
    floatmode='fixed')
  np.random.seed(0)  # not used this version

  # 1. load data
  print("\nLoading People train (200), test (400) data ")
 
  train_file = ".\\Data\\People\\people_train_200.txt"
  test_file = ".\\Data\\People\\people_test_40.txt"
  # sex, age, city, height, psych(3)
  x_cols=[0,1,2,3, 4,5,6]
  target_col = 7  # 1.0000 .6413

  train_X = np.loadtxt(train_file, comments="#",
    usecols=x_cols, delimiter=",",  dtype=np.float64)
  train_y = np.loadtxt(train_file, comments="#",
    usecols=target_col, delimiter=",",  dtype=np.float64)

  test_X = np.loadtxt(test_file, comments="#",
    usecols=x_cols, delimiter=",",  dtype=np.float64)
  test_y = np.loadtxt(test_file, comments="#",
    usecols=target_col, delimiter=",",  dtype=np.float64)
  print("Done ")

  print("\nFirst three raw train X: ")
  print("M,29,Austin,68.5,timid,54600.00 ")
  print("F,51,Juneau,60.0,brave,25700.00 ")
  print("M,32,Boston,73.5,alert,74900.00 ")

  print("\nFirst three X predictors: ")
  print(train_X[0:3,:])
  print("\nFirst three y targets: ")
  for i in range(3):
    print("%0.4f" % train_y[i])

  n_estimators = 200
  max_depth = 4
  min_samples_split = 2
  min_samples_leaf = 1
  lrn_rate = 0.05
  seed = 0

  print("\nSetting n_estimators = " + \
    str(n_estimators))
  print("Setting max_depth = " + \
    str(max_depth))
  print("Setting min_samples_split = " + \
    str(min_samples_split))
  print("Setting min_samples_leaf = " + \
    str(min_samples_leaf))
  print("Setting lrn_rate = %0.4f " % lrn_rate)

  # GradientBoostingRegressor(*, loss='squared_error',
  # learning_rate=0.1, n_estimators=100, subsample=1.0,
  # criterion='friedman_mse', min_samples_split=2,
  # min_samples_leaf=1, min_weight_fraction_leaf=0.0,
  # max_depth=3, min_impurity_decrease=0.0, init=None,
  # random_state=None, max_features=None, alpha=0.9,
  # verbose=0, max_leaf_nodes=None, warm_start=False,
  # validation_fraction=0.1, n_iter_no_change=None,
  # tol=0.0001, ccp_alpha=0.0)

  model = \
    GradientBoostingRegressor(max_depth=max_depth, 
    min_samples_split=min_samples_split,
    n_estimators=n_estimators,
    min_samples_leaf=min_samples_leaf,
    learning_rate=lrn_rate, random_state=seed)

  print("\nTraining gradient boost regression model ")
  model.fit(train_X, train_y)
  print("Done " )

  print("\nEvaluating model ")
  acc_train = accuracy(model, train_X, train_y, 0.10)
  print("\nAcc train (within 0.10): %0.4f " % acc_train)
  acc_test = accuracy(model, test_X, test_y, 0.10)
  print("Acc test (within 0.10): %0.4f " % acc_test)

  mse_train = MSE(model, train_X, train_y)
  print("\nMSE train: %0.4f " % mse_train)
  mse_test = MSE(model, test_X, test_y)
  print("MSE test: %0.4f " % mse_test)

  r2_train = model.score(train_X, train_y)
  print("\nR2 train: %0.4f " % r2_train)
  r2_test = model.score(test_X, test_y)
  print("R2 test: %0.4f " % r2_test)

  print("\nPredicting for x = ")
  print(train_X[0])
  pred_y = model.predict(train_X[0].reshape(1,-1))[0]
  print("Predicted y = %0.4f " % pred_y)

  print("\nEnd demo ")

if __name__ == "__main__":
  main()

Training data:

# people_train_200.txt
#
# M = 0,  F = 1
# city_target_encoding = {
#     "Boston": 30600.00,  0
#     "Casper": 33100.00,  1
#     "Austin": 33961.90,  2
#     "Denver": 35535.00,  3
#     "Eugene": 44090.00,  4
#     "Irvine": 51840.00,  5
#     "Fresno": 52950.00,  6
#     "Gallup": 52960.00,  7
#     "Juneau": 58842.86,  8
#     "Helena": 64194.74,  9
# }
# timid = 1,0,0  alert = 0,1,0  brave = 0,0,1
#
# Sex Age City Height Psych Income
#
0, 0.29, 2, 0.685, 1,0,0, 0.54600
1, 0.51, 8, 0.600, 0,0,1, 0.25700
0, 0.32, 0, 0.735, 0,1,0, 0.74900
0, 0.19, 5, 0.610, 0,0,1, 0.30300
1, 0.45, 4, 0.645, 0,1,0, 0.35600
0, 0.58, 6, 0.720, 1,0,0, 0.52100
1, 0.22, 2, 0.630, 0,1,0, 0.27800
0, 0.64, 8, 0.775, 0,0,1, 0.88200
1, 0.38, 3, 0.660, 1,0,0, 0.32400
0, 0.41, 7, 0.700, 0,1,0, 0.51200
1, 0.27, 1, 0.625, 0,0,1, 0.28900
0, 0.53, 9, 0.740, 1,0,0, 0.62300
1, 0.34, 0, 0.650, 0,1,0, 0.33100
0, 0.60, 5, 0.765, 0,0,1, 0.79800
1, 0.21, 4, 0.605, 1,0,0, 0.22400
0, 0.48, 6, 0.690, 0,1,0, 0.48700
1, 0.62, 2, 0.670, 0,0,1, 0.53600
0, 0.25, 8, 0.715, 1,0,0, 0.45200
1, 0.39, 3, 0.640, 0,1,0, 0.31800
0, 0.55, 7, 0.750, 0,0,1, 0.72100
1, 0.31, 1, 0.615, 1,0,0, 0.24500
0, 0.43, 9, 0.705, 0,1,0, 0.52400
1, 0.50, 0, 0.635, 0,0,1, 0.38900
0, 0.28, 5, 0.730, 0,1,0, 0.58100
1, 0.56, 4, 0.665, 1,0,0, 0.47800
0, 0.20, 6, 0.620, 0,0,1, 0.29100
1, 0.42, 2, 0.655, 0,1,0, 0.37200
0, 0.66, 8, 0.780, 1,0,0, 0.85400
1, 0.36, 3, 0.620, 0,0,1, 0.29800
0, 0.49, 7, 0.725, 0,1,0, 0.59300
1, 0.24, 1, 0.600, 0,1,0, 0.21200
0, 0.57, 9, 0.755, 0,0,1, 0.76400
1, 0.46, 0, 0.640, 1,0,0, 0.36100
0, 0.33, 5, 0.680, 0,1,0, 0.44300
1, 0.61, 4, 0.675, 0,0,1, 0.57900
0, 0.23, 6, 0.710, 1,0,0, 0.42800
1, 0.52, 2, 0.630, 0,1,0, 0.39400
0, 0.65, 8, 0.760, 0,0,1, 0.81200
1, 0.30, 3, 0.610, 1,0,0, 0.23100
0, 0.44, 7, 0.695, 0,1,0, 0.49600
1, 0.59, 1, 0.680, 0,0,1, 0.58300
0, 0.35, 9, 0.720, 1,0,0, 0.49800
1, 0.18, 0, 0.605, 0,1,0, 0.21900
0, 0.51, 5, 0.745, 0,0,1, 0.66700
1, 0.40, 4, 0.650, 1,0,0, 0.33500
0, 0.67, 6, 0.785, 0,1,0, 0.94200
1, 0.26, 2, 0.625, 0,0,1, 0.27400
0, 0.37, 8, 0.690, 0,1,0, 0.45100
1, 0.54, 3, 0.660, 1,0,0, 0.43600
0, 0.30, 7, 0.730, 0,0,1, 0.59800
1, 0.47, 1, 0.645, 0,1,0, 0.38200
0, 0.63, 9, 0.770, 1,0,0, 0.83100
1, 0.28, 0, 0.615, 0,0,1, 0.26300
0, 0.42, 5, 0.700, 0,1,0, 0.48900
1, 0.35, 4, 0.635, 1,0,0, 0.28700
0, 0.59, 6, 0.750, 0,0,1, 0.74600
1, 0.20, 2, 0.610, 0,1,0, 0.22800
0, 0.46, 8, 0.720, 1,0,0, 0.53200
1, 0.58, 3, 0.670, 0,1,0, 0.51400
0, 0.31, 7, 0.715, 0,0,1, 0.52900
1, 0.33, 1, 0.620, 1,0,0, 0.25800
0, 0.54, 9, 0.735, 0,1,0, 0.64500
1, 0.41, 0, 0.655, 0,0,1, 0.36700
0, 0.22, 5, 0.630, 1,0,0, 0.26400
1, 0.67, 4, 0.685, 0,1,0, 0.69100
0, 0.38, 6, 0.695, 0,0,1, 0.47900
1, 0.25, 2, 0.600, 1,0,0, 0.19800
0, 0.61, 8, 0.760, 0,1,0, 0.77400
1, 0.43, 3, 0.640, 0,0,1, 0.35200
0, 0.26, 7, 0.700, 1,0,0, 0.43100
1, 0.48, 1, 0.650, 0,1,0, 0.39700
0, 0.50, 9, 0.725, 0,0,1, 0.60200
1, 0.29, 0, 0.620, 1,0,0, 0.24900
0, 0.68, 5, 0.785, 0,1,0, 0.96800
1, 0.37, 4, 0.630, 0,0,1, 0.30600
0, 0.18, 6, 0.615, 1,0,0, 0.25300
1, 0.65, 2, 0.680, 0,1,0, 0.63200
0, 0.34, 8, 0.685, 0,0,1, 0.43800
1, 0.53, 3, 0.655, 1,0,0, 0.41900
0, 0.45, 7, 0.710, 0,1,0, 0.53700
1, 0.23, 1, 0.610, 0,0,1, 0.23600
0, 0.52, 9, 0.740, 1,0,0, 0.61100
1, 0.32, 0, 0.630, 0,1,0, 0.28400
0, 0.62, 5, 0.770, 0,0,1, 0.82500
1, 0.19, 4, 0.600, 1,0,0, 0.20100
0, 0.40, 6, 0.700, 0,1,0, 0.47300
1, 0.60, 2, 0.675, 0,0,1, 0.56800
0, 0.27, 8, 0.695, 1,0,0, 0.41200
1, 0.44, 3, 0.645, 0,1,0, 0.36900
0, 0.56, 7, 0.755, 0,0,1, 0.75300
1, 0.36, 1, 0.635, 1,0,0, 0.29100
0, 0.47, 9, 0.715, 0,1,0, 0.55800
1, 0.55, 0, 0.660, 0,0,1, 0.46200
0, 0.21, 5, 0.625, 0,1,0, 0.28900
1, 0.63, 4, 0.680, 1,0,0, 0.58400
0, 0.39, 6, 0.690, 0,0,1, 0.46100
1, 0.22, 2, 0.615, 1,0,0, 0.21700
0, 0.64, 8, 0.770, 0,1,0, 0.85900
1, 0.38, 3, 0.635, 0,0,1, 0.33800
0, 0.29, 7, 0.705, 1,0,0, 0.45800
1, 0.51, 1, 0.655, 0,1,0, 0.43100
0, 0.58, 9, 0.760, 0,0,1, 0.78200
1, 0.27, 0, 0.610, 0,1,0, 0.24300
0, 0.36, 5, 0.685, 1,0,0, 0.42100
1, 0.66, 4, 0.685, 0,0,1, 0.67500
0, 0.24, 6, 0.640, 0,1,0, 0.31200
1, 0.49, 2, 0.650, 1,0,0, 0.38400
0, 0.60, 8, 0.755, 0,0,1, 0.76900
1, 0.31, 3, 0.625, 0,1,0, 0.26800
0, 0.43, 7, 0.700, 1,0,0, 0.48200
1, 0.57, 1, 0.665, 0,0,1, 0.51900
0, 0.32, 9, 0.690, 0,1,0, 0.43900
1, 0.20, 0, 0.600, 1,0,0, 0.19200
0, 0.53, 5, 0.750, 0,1,0, 0.68400
1, 0.42, 4, 0.640, 0,0,1, 0.36100
0, 0.65, 6, 0.775, 1,0,0, 0.84700
1, 0.28, 2, 0.620, 0,1,0, 0.26200
0, 0.35, 8, 0.680, 1,0,0, 0.39600
1, 0.50, 3, 0.650, 0,0,1, 0.41800
0, 0.25, 7, 0.690, 0,1,0, 0.41300
1, 0.45, 1, 0.640, 1,0,0, 0.34800
0, 0.59, 9, 0.765, 0,1,0, 0.80100
1, 0.33, 0, 0.625, 0,0,1, 0.30200
0, 0.23, 5, 0.635, 0,1,0, 0.29700
1, 0.64, 4, 0.680, 0,0,1, 0.62100
0, 0.41, 6, 0.705, 1,0,0, 0.47100
1, 0.21, 2, 0.605, 0,0,1, 0.21800
0, 0.67, 8, 0.785, 0,1,0, 0.95600
1, 0.37, 3, 0.630, 1,0,0, 0.29200
0, 0.48, 7, 0.720, 0,0,1, 0.59100
1, 0.26, 1, 0.615, 0,1,0, 0.26100
0, 0.55, 9, 0.745, 1,0,0, 0.65800
1, 0.40, 0, 0.645, 0,1,0, 0.34900
0, 0.19, 5, 0.615, 0,0,1, 0.28400
1, 0.68, 4, 0.690, 1,0,0, 0.70200
0, 0.36, 6, 0.680, 0,1,0, 0.42600
1, 0.29, 2, 0.625, 0,0,1, 0.29300
0, 0.62, 8, 0.765, 1,0,0, 0.79300
1, 0.44, 3, 0.640, 0,1,0, 0.35700
0, 0.28, 7, 0.700, 0,0,1, 0.46900
1, 0.52, 1, 0.660, 1,0,0, 0.42100
0, 0.46, 9, 0.710, 0,1,0, 0.53900
1, 0.24, 0, 0.605, 0,0,1, 0.20700
0, 0.57, 5, 0.755, 1,0,0, 0.75100
1, 0.35, 4, 0.630, 0,1,0, 0.29600
0, 0.66, 6, 0.780, 0,0,1, 0.88900
1, 0.30, 2, 0.620, 0,0,1, 0.29700
0, 0.38, 8, 0.690, 1,0,0, 0.42400
1, 0.56, 3, 0.665, 0,1,0, 0.48300
0, 0.27, 7, 0.695, 0,0,1, 0.44800
1, 0.43, 1, 0.635, 1,0,0, 0.32800
0, 0.61, 9, 0.770, 0,1,0, 0.83600
1, 0.31, 0, 0.620, 0,0,1, 0.28100
0, 0.20, 5, 0.620, 1,0,0, 0.24600
1, 0.65, 4, 0.685, 0,1,0, 0.66400
0, 0.42, 6, 0.700, 0,0,1, 0.51200
1, 0.18, 2, 0.600, 0,1,0, 0.18900
0, 0.63, 8, 0.775, 1,0,0, 0.82900
1, 0.39, 3, 0.635, 0,1,0, 0.32600
0, 0.50, 7, 0.725, 0,0,1, 0.63400
1, 0.25, 1, 0.610, 1,0,0, 0.22300
0, 0.54, 9, 0.740, 0,1,0, 0.65200
1, 0.36, 0, 0.635, 1,0,0, 0.29800
0, 0.31, 5, 0.675, 0,0,1, 0.42900
1, 0.58, 4, 0.670, 0,1,0, 0.52800
0, 0.25, 6, 0.645, 1,0,0, 0.31900
1, 0.47, 2, 0.645, 0,0,1, 0.41100
0, 0.59, 8, 0.750, 0,1,0, 0.73800
1, 0.28, 3, 0.615, 1,0,0, 0.23700
0, 0.40, 7, 0.690, 0,1,0, 0.45600
1, 0.54, 1, 0.650, 0,0,1, 0.46300
0, 0.34, 9, 0.695, 1,0,0, 0.43200
1, 0.23, 0, 0.600, 0,1,0, 0.20800
0, 0.52, 5, 0.740, 0,0,1, 0.66100
1, 0.41, 4, 0.640, 1,0,0, 0.32900
0, 0.68, 6, 0.785, 0,1,0, 0.97800
1, 0.26, 2, 0.610, 0,0,1, 0.23800
0, 0.37, 8, 0.685, 0,1,0, 0.40900
1, 0.50, 3, 0.655, 1,0,0, 0.42700
0, 0.29, 7, 0.700, 0,0,1, 0.46300
1, 0.46, 1, 0.640, 0,1,0, 0.35400
0, 0.60, 9, 0.760, 1,0,0, 0.80800
1, 0.28, 0, 0.610, 0,0,1, 0.26800
0, 0.43, 5, 0.705, 0,1,0, 0.47600
1, 0.34, 4, 0.625, 1,0,0, 0.28300
0, 0.56, 6, 0.745, 0,0,1, 0.70400
1, 0.19, 2, 0.600, 0,1,0, 0.17400
0, 0.65, 8, 0.770, 1,0,0, 0.86800
1, 0.42, 3, 0.635, 0,0,1, 0.35900
0, 0.30, 7, 0.695, 0,1,0, 0.42200
1, 0.53, 1, 0.655, 1,0,0, 0.44800
0, 0.47, 9, 0.720, 0,0,1, 0.59400
1, 0.22, 0, 0.600, 0,1,0, 0.19800
0, 0.66, 5, 0.780, 1,0,0, 0.89200
1, 0.38, 4, 0.630, 0,0,1, 0.31800
0, 0.21, 6, 0.620, 0,1,0, 0.27100
1, 0.62, 2, 0.665, 1,0,0, 0.50300
0, 0.33, 8, 0.675, 0,0,1, 0.41200
1, 0.51, 3, 0.650, 0,1,0, 0.41500
0, 0.44, 7, 0.705, 1,0,0, 0.48800

Test data:

# people_test_40.txt
#
1, 0.27, 1, 0.620, 0,0,1, 0.28100
0, 0.58, 9, 0.750, 0,1,0, 0.75900
1, 0.32, 0, 0.625, 1,0,0, 0.25900
0, 0.24, 5, 0.640, 0,0,1, 0.32400
1, 0.60, 4, 0.660, 0,1,0, 0.52100
0, 0.39, 6, 0.685, 1,0,0, 0.42700
0, 0.26, 2, 0.690, 0,1,0, 0.40200
1, 0.49, 8, 0.645, 0,0,1, 0.38800
0, 0.53, 3, 0.730, 1,0,0, 0.61200
1, 0.30, 7, 0.615, 0,1,0, 0.23400
0, 0.61, 1, 0.760, 0,0,1, 0.78900
1, 0.22, 9, 0.610, 1,0,0, 0.21500
0, 0.45, 0, 0.700, 0,1,0, 0.49300
1, 0.57, 5, 0.670, 0,0,1, 0.51100
0, 0.28, 4, 0.695, 1,0,0, 0.40800
1, 0.39, 6, 0.630, 0,1,0, 0.31100
0, 0.51, 2, 0.720, 0,0,1, 0.57400
1, 0.33, 8, 0.625, 1,0,0, 0.26700
0, 0.64, 3, 0.770, 0,1,0, 0.83500
1, 0.18, 7, 0.600, 0,0,1, 0.16200
0, 0.42, 1, 0.695, 1,0,0, 0.46100
1, 0.59, 9, 0.680, 0,1,0, 0.56300
0, 0.36, 0, 0.685, 0,0,1, 0.43400
1, 0.23, 5, 0.610, 1,0,0, 0.22100
0, 0.67, 4, 0.780, 0,1,0, 0.92400
1, 0.45, 6, 0.645, 0,0,1, 0.37800
0, 0.31, 2, 0.700, 1,0,0, 0.43500
1, 0.52, 8, 0.655, 0,1,0, 0.43900
0, 0.20, 3, 0.620, 0,0,1, 0.26200
1, 0.48, 7, 0.650, 1,0,0, 0.37100
0, 0.55, 1, 0.740, 0,1,0, 0.63800
1, 0.29, 9, 0.620, 0,0,1, 0.27300
0, 0.63, 0, 0.765, 1,0,0, 0.80100
1, 0.41, 5, 0.640, 0,1,0, 0.34200
0, 0.22, 4, 0.630, 0,0,1, 0.28600
1, 0.64, 6, 0.675, 1,0,0, 0.61400
0, 0.37, 2, 0.690, 0,1,0, 0.44100
1, 0.50, 8, 0.650, 0,0,1, 0.42200
0, 0.27, 3, 0.695, 1,0,0, 0.39800
1, 0.44, 7, 0.640, 0,1,0, 0.33900

Data smoothing is a classical statistics technique that removes noise from a dataset so that overall patterns and trends are visible. There are many algorithms for smoothing. I put together a demo of data smoothing using the Nadaraya-Watson (NW) kernel smoothing technique, with the C# language.

The output of the demo program is:

Begin Nadaraya-Watson kernel smoothing demo

Loading simple XXy data
Done

First three data X:
  0.1000  0.1000
  0.2000  0.2000
  0.3000  0.3000

First three data y:
  0.0998
  0.3987
  0.0955

Creating NW kernel smoother with RBF sigma = 0.1600
Done

Computing smoothed y values
Done

   X             y        smooth y
-------------------------------------
  0.10  0.10    0.0998    0.2069
  0.20  0.20    0.3987    0.2303
  0.30  0.30    0.0955    0.2474
  0.50  0.50    0.5794    0.5599
  0.60  0.60    0.5646    0.6681
  0.70  0.70    0.9442    0.7357
  0.90  0.90    0.4833    0.5556
  1.00  1.00    0.5415    0.5255

End NW kernel smoothing demo

The source data is:

0.1000, 0.1000, 0.0998
0.2000, 0.2000, 0.3987
0.3000, 0.3000, 0.0955
0.5000, 0.5000, 0.5794
0.6000, 0.6000, 0.5646
0.7000, 0.7000, 0.9442
0.9000, 0.9000, 0.4833
1.0000, 1.0000, 0.5415

There are eight data items. The first two values on each line are the elements of a two-dimensional X vector, such as (0.10, 0.10). The third value is a dependent y value. The data is artificial but you can imagine the X values are inputs to some system and the y values are outputs.

Data smoothing is most often used with small datasets that have a single x value and a single y value. This allows the source data and the smoothed data to be visualized using a graph. However, I used data with two x values, such as (0.10, 0.10), where both values are artificially the same, to demonstrate that NW kernel smoothing can use multidimensional X data with several values.

Understanding the Radial Basis Function

In order to understand NW kernel smoothing, you must understand the radial basis function (RBF) function. The RBF function computes a measure of similarity between two vectors. There are several RBF versions. The version used by the demo is often called the sigma (bandwidth) version:

rbf(v1, v2) = exp(-1 * ||v1 - v2||^2 / (2 * sigma^2))

The v1 and v2 are vectors. The exp() function is Euler’s constant (approximately 2.71828) raised to a power. The ||v1 – v2||^2 is the squared Euclidean distance between vectors v1 and v2. The sigma, also called the bandwidth, is an arbitrary constant.

If the two vectors passed to RBF are the same, then rbf(v1, v2) = 1.0 (maximum similarity). As the two vectors get more and more different, the value of rbf(v1, v2) gets closer and closer to 0. RBF is symmetric so rbf(v1, v2) = rbf(v2, v1).

Suppose v1 = (2.50, -3.25, 1.00) and v2 = (2.0, -3.0, 1.0) and sigma = 1.5, then the squared distance term is:

||v1 - v2||^2 = (2.50 - 2)^2 + (-3.25 - (-3))^2 + (1.00 - 1)^2
              = (0.50)^2 + (-0.25)^2 + (0.00)^2
              = 0.2500 + 0.0625 + 0.0000
              = 0.3125

and the RBF value is:

rbf(v1, v2) = exp( -1 * ||v1 - v2||^2 / (2 * sigma^2) )
            = exp( -0.3125 / (2 * 1.5^2) )
            = exp( -0.3125 / 4.50 )
            = exp(-0.0694)

Computing a NW Kernel Smoothed y Value

Computing a NW smoothed data point is best explained using a concrete example. Suppose you want to compute the smoothed y’ value for x = (0.70, 0.70) in the demo data. The actual y value is 0.9442. The smoothed y’ value is 0.7357. It is computed as follows:

x = (0.70, 0.70) 

x0 = (0.10, 0.10)  w0 = rbf(x, x0) = 0.0000  y0 = 0.0998  w0 * y0 = 0.0000
x1 = (0.20, 0.20)  w1 = rbf(x, x1) = 0.0001  y1 = 0.3987  w1 * y1 = 0.0000
x2 = (0.30, 0.30)  w2 = rbf(x, x2) = 0.0019  y2 = 0.0955  w2 * y2 = 0.0001
x3 = (0.50, 0.50)  w3 = rbf(x, x3) = 0.2096  y3 = 0.5794  w3 * y3 = 0.1214
x4 = (0.60, 0.60)  w4 = rbf(x, x4) = 0.6766  y4 = 0.5646  w4 * y4 = 0.3820
x5 = (0.70, 0.70)  w5 = rbf(x, x5) = 1.0000  y5 = 0.9442  w5 * y5 = 0.9442
x6 = (0.90, 0.90)  w6 = rbf(x, x6) = 0.2096  y6 = 0.4833  w6 * y6 = 0.1013
x7 = (1.00, 1.00)  w7 = rbf(x, x7) = 0.0297  y7 = 0.5415  w7 * y7 = 0.0161
                                    -------                        -------
                       sum weights = 2.1276        sum weighted y = 1.5653

smooth y' = sum weighted y / sum weights
          = 1.5653 / 2.1276
          = 0.7357

There are eight data points. A weight, wi, is computed for each data point, where the weights are the RBF values between the x value and each data points. Each weight is multiplied by associated actual y value, and those products are summed to give the numerator term of the smoothed value. The denominator is the sum of the weight values.

I miss the days when I used to read a daily newspaper to get all my news and information. Here is some old health advice.

Demo program. Replace “lt” (less than), “gt”, “lte”, “gte” with Boolean operator symbols (my blog editor chokes on symbols).

using System;
using System.IO;
using System.Collections.Generic;

namespace KernelNWSmoothing
{
  internal class KernelNWSmoothingProgram
  {
    static void Main(string[] args)
    {
      Console.WriteLine("\nBegin Nadaraya-Watson" +
        " kernel smoothing demo ");

      Console.WriteLine("\nLoading simple XXy data ");
      string dataFile = "..\\..\\..\\Data\\simple_xxy_8.txt";
      double[][] X = MatLoad(dataFile,
        new int[] { 0, 1 }, ',', "#");
      double[] y = MatToVec(MatLoad(dataFile,
        new int[] { 2 }, ',', "#"));
      Console.WriteLine("Done ");

      Console.WriteLine("\nFirst three data X: ");
      for (int i = 0; i "lt" 3; ++i)
        VecShow(X[i], 4, 8, true);

      Console.WriteLine("\nFirst three data y: ");
      for (int i = 0; i "lt" 3; ++i)
        Console.WriteLine(y[i].ToString("F4").PadLeft(8));

      double sigma = 0.16; // larger sigma  = smoother
      Console.WriteLine("\nCreating NW kernel smoother" +
        " with RBF sigma = " + sigma.ToString("F4"));
      KernelSmoother ks = new KernelSmoother("rbf", sigma);
      Console.WriteLine("Done ");

      Console.WriteLine("\nComputing smoothed y values ");
      ks.Fit(X, y); // copy dataX and dataY
      Console.WriteLine("Done ");

      Console.WriteLine("\n   X         y        smooth y");
      Console.WriteLine("--------------------------------");
      for (int i = 0; i "lt" X.Length; ++i)
      {
        double[] x = X[i];
        double actY = y[i];
        double smoothY = ks.SmoothedY(x);
        VecShow(x, 2, 6, false);
        Console.Write("    " + actY.ToString("F4"));
        Console.WriteLine("    " + smoothY.ToString("F4"));
      }

      Console.WriteLine("\nEnd NW kernel smoothing demo ");
      Console.ReadLine();
    } // Main

    // ------------------------------------------------------
    // helpers for Main(): MatLoad(), MatToVec(), VecShow()
    // ------------------------------------------------------

    static double[][] MatLoad(string fn, int[] usecols,
      char sep, string comment)
    {
      List"lt"double[]"gt" result =
        new List"lt"double[]"gt"();
      string line = "";
      FileStream ifs = new FileStream(fn, FileMode.Open);
      StreamReader sr = new StreamReader(ifs);
      while ((line = sr.ReadLine()) != null)
      {
        if (line.StartsWith(comment) == true)
          continue;
        string[] tokens = line.Split(sep);
        List"lt"double"gt" lst = new List"lt"double"gt"();
        for (int j = 0; j "lt" usecols.Length; ++j)
          lst.Add(double.Parse(tokens[usecols[j]]));
        double[] row = lst.ToArray();
        result.Add(row);
      }
      sr.Close(); ifs.Close();
      return result.ToArray();
    }

    static double[] MatToVec(double[][] M)
    {
      int nRows = M.Length;
      int nCols = M[0].Length;
      double[] result = new double[nRows * nCols];
      int k = 0;
      for (int i = 0; i "lt" nRows; ++i)
        for (int j = 0; j "lt" nCols; ++j)
          result[k++] = M[i][j];
      return result;
    }

    static void VecShow(double[] vec, int dec, int wid,
      bool newLine)
    {
      for (int i = 0; i "lt" vec.Length; ++i)
        Console.Write(vec[i].ToString("F" + dec).
          PadLeft(wid));
      if (newLine == true)
        Console.WriteLine("");
    }

  } // class Program

  // ========================================================

  public class KernelSmoother
  {
    public string kernel; // 
    public double sigma;  // RBF sigma, aka bandwidth
    public double[][] dataX;
    public double[] dataY;

    public KernelSmoother(string kernel, double sigma)
    {
      if (kernel != "rbf")
        Console.WriteLine("Only rbf kernel supported ");
      this.kernel = kernel; // dummy not used this version
      this.sigma = sigma;
    }

    public void Fit(double[][] dataX, double[] dataY)
    {
      this.dataX = dataX; // copy by reference
      this.dataY = dataY;
    }

    public double SmoothedY(double[] x)
    {
      int n = this.dataX.Length;
      double numer = 0.0; // sum of weighted y values
      double denom = 0.0; // sum of weights (RBF values)
      for (int i = 0; i "lt" n; ++i)
      {
        double currRbf = RbfSigma(x, this.dataX[i],
          this.sigma);
        numer += currRbf * this.dataY[i];
        denom += currRbf;
      }
      return numer / denom; // assume no div by 0
    }

    private static double RbfSigma(double[] x1, double[] x2,
      double sigma)
    {
      // rbf(x1,x2) = exp(-1 * ||x1-x2||^2 / (2 * sigma^2))
      int dim = x1.Length;
      double sum = 0.0;
      for (int i = 0; i "lt" dim; ++i)
        sum += (x1[i] - x2[i]) * (x1[i] - x2[i]);
      return Math.Exp(-1 * sum / (2 * sigma * sigma));
    }

    private static double RbfGamma(double[] x1, double[] x2,
      double gamma)
    {
      // rbf(x1,x2) = exp(-1 * gamma * ||x1-x2||^2
      int dim = x1.Length;
      double sum = 0.0;
      for (int i = 0; i "lt" dim; ++i)
        sum += (x1[i] - x2[i]) * (x1[i] - x2[i]);
      return Math.Exp(-1 * gamma * sum);
    }

  } // class KernelSmoother

} // ns

Demo data:

# simple_xxy_8.txt
#
0.1000, 0.1000, 0.0998
0.2000, 0.2000, 0.3987
0.3000, 0.3000, 0.0955
0.5000, 0.5000, 0.5794
0.6000, 0.6000, 0.5646
0.7000, 0.7000, 0.9442
0.9000, 0.9000, 0.4833
1.0000, 1.0000, 0.5415