Mathematica Implementations of the Random Forest algorithm

I'm going to be bold and attempt to edit the Ross code so that it is (a) a little easier to understand and (b) takes the same form of argument as LinearModelFit and other Mathematica prediction creators. I've also added some annotations to the critical code. My variable names are now far longer than the Ross names but perhaps for informative. So far in my testing this code works the same as Ross, but it is possible I have messed something up in my rewrite.

Part 1. Same as Ross but I add a function informationGain that determines the information gain one obtains about y from knowing the ith feature of x.

condEnt = Statistics`Library`NConditionalEntropy;
ent = Statistics`Library`NEntropy;
maxI[v_] := Ordering[v, -1][[1]];
informationGain[x_, y_, i_] := 
  Statistics`Library`NEntropy[x] - 
  Statistics`Library`NConditionalEntropy[x[[i]], y]

Part 2. Same idea as Ross but x now contains not just the features (independent variables in statistics lingo) but the entire data set. It thus has the same structure as the first argument to LinearModelFit and NonlinearModelFit.

  cart[x_?MatrixQ /; ent[Last /@ x] == 0, m_] := First[Last /@ x]

Part 3. The idea here is again to let the first argument just be the entire data set. What I have then done is to use a lot of local variables with hopefully evocative names and some annotations to better explain what is going on in the elegant but terse Ross code. The text below explains the annotations.

cart[x_?MatrixQ, m_Integer] := 
 Block[{dims = Dimensions[x], numberOfInstances, numberOfAttributes, 
allAttributes, y, byFeature,
h, keep, bestVar, ub, allValuesOfBestVar, mask, best, 
xBiggerThanBestSplitQ, bestSplittingValue, 
iXBestGTESplit,
iXBestLTSplit},
numberOfInstances = First[dims];
numberOfAttributes = Last[dims] - 1;(*1*)
allAttributes = Range[numberOfAttributes];(*2*)
y = Part[x, All, -1];(*3*)
byFeature = Transpose[(Most /@ x)];(*4*)
h = ent[y];(*5*)
keep = RandomSample[allAttributes, m];(*6*)
bestVar = 
maxI[Table[
 If[MemberQ[keep, i], informationGain[byFeature, y, i], 0], {i, 
  allAttributes}]];(*7*)
allValuesOfBestVar = byFeature[[bestVar]];(*8*)
ub = Union[allValuesOfBestVar];(*9*)
mask = UnitStep[allValuesOfBestVar - #] & /@ 
ub;(*10*)
best = maxI[(h - condEnt[#, y] & /@ mask)];(*11*)
xBiggerThanBestSplitQ = mask[[best]];(*12*)
bestSplittingValue = ub[[best]];(*13*)
iXBestGTESplit= 
  Pick[x, xBiggerThanBestSplitQ, 1];(*14*)
iXBestLTSplit = 
  Pick[x, xBiggerThanBestSplitQ, 0];(*15*)
If[Min[xBiggerThanBestSplitQ] === Max[xBiggerThanBestSplitQ], 
   RandomChoice[y] , 
{bestVar, bestSplittingValue, 
cart[iXBestGTESplit 
 m], 
cart[
iXBestLTSplit, m]}](*16*)
]

Part 3 Annotations.

1. The number of attributes is just the number of columns minus one, because the last column is the class value.

2. I like the idea of iterating over a set of values, so I create it here.

3. The class values y are the last column of the x matrix.

4. To get the features in the same form as Ross, Transpose x and then take Most of each row.

5. Calculating the entropy of the y values in this level of the decision tree.

6. Use an m-length subset of the attributes. Note that in my code, m is the number to keep, not the number to drop.

7. Find the information gain for each selected attribute. Find the selected attribute that produces the greatest information gain. This produces an index into the attributes. We call it bestVar.

8. Get the features in the bestVar th column of x. Call it allValuesOfBestVar.

9. Find the different values of the features in allValuesOfBestVar. Call it ub. (Union of Best values). These now represent possible splitting points in our best feature.

10. This is the clever part. For each possible splitting point, determine whether each value is above or equal to that splitting point (1) or below it (9). This produces a matrix that is Length[ub] x Length[y]. Call this matrix mask.

11. Rather than compute the information gain from knowing the feature values, compute the information gain from knowing whether the feature for a particular instance is above a splitting point for all the splitting points. Pick the splitting point that produces the most information gain. The index of the splitting point into the list of splitting points is named best. Note that this is the second use of the idea of conditional entropy.

12. Now select the best part of the mask. This gives us the information on >= or < using the best splitting point for the data for the best feature of each instance.

13. We'll need that best splitting value again, so we capture it here as bestSplittingValue.

14. Pick the instances for which the best feature has a value above or equal to the best splitting value. Call those instances iXBestGTESplit

15. Pick the instances for which the best feature has a value less than the best splitting value. Call those instances iXBestLTSplit

16. The recursion step. Run cart again first on the iXBestGTESplit data and then on the iXBestLTSplit data. But before doing so, create a list that shows what feature you split on and what the best splitting value was. Stop when your data has no more information to yield.

Part 4. This code to classify an instance works exactly the same as in Ross. Notice that here, and now consistent with the way it is done above, x is a vector of attributes. It will also work as a vector that contains a list of which most of the elements are attributes and the last element is the class value. Thus, one can now just map over the instances in some testing set to get the predictions.

classify[x_?MatrixQ, Except[_List, d_]] := d
classify[x_?MatrixQ, {best_, value_, d1_, d2_}] := 
  classify[x, If[x[[best]] < value, d2, d1]]

Part 5. The random forest creation function is similar to Ross but, again, x is now a unified matrix containing the instances. I have added an optional argument subsetFactor that can speed up your code (at the risk of losing accuracy) by using just a subset of the instances to create the trees each time.

rf[x_?MatrixQ, m_Integer, ntree_Integer, subsetFactor_: 1] := 
Table[With[{boot = 
 RandomChoice[Range[Length[x]], Round[subsetFactor*Length[x]]]},
 cart[x[[boot]], m]], {ntree}]

I very much enjoy Dan's approach in part because it is so simple both in concept and implementation. I'm taking the liberty here of suggesting a few arguable improvements to his terrific code. For makeForest (a) the data is in the same format as is used in functions such as LinearModelFit (a simple array instead of a list of rules of features onto class); (b) standard machine learning vocabulary in naming variables; (c) direct use of RandomSample on the data.

makeForest[data_, ntrees_, subdim_, sizes_] := 
  Module[{dims = Dimensions[data], numberOfInstances, 
  numberOfAttributes, leaves, slots, subleaves, nf},
  numberOfInstances = First[dims];
  numberOfAttributes = Last[dims] - 1;
  Table[(leaves = RandomSample[data, sizes];
         slots = Sort[RandomSample[Range[dim], subdim]];
         subleaves = Map[#[[slots]] -> #[[-1]] &, leaves];
         nf = Nearest[subleaves];
         {slots, nf}), {ntrees}]
   ]

For classify,(a) clarify with named variables how one extracts the slots and the NearestFunction from the forest; (b) use the built in Commonest instead of the Reverse SortBy Tally composition. This method does lose the number of votes each prediction received; (c) create an optional argument k that effectively implements kNN on the data instead of Dan's required 1NN.

classify[instance_, forest_, n_, k_: 1] := Module[{predictions},
  predictions = Flatten[Table[Module[{tree, nf, slots},
  tree = forest[[j]];
  slots = tree[[1]];
  nf = tree[[2]];
  nf[instance[[slots]], k]], {j, Length[forest]}]];
  Commonest[predictions, n]]

Also, a wacky idea that follows up on Andy Ross's comment: could one not pretty easily create an ensemble of forests and then let the forests that perform well on the training data have breeding rights? Breeding might consist of some sort of reshuffling of the slots. To do this, we might take more liberties with Dan's code by creating a polymorphic makeForest that permits the slots to be input directly.

makeForest[data_, slotList_, sizes_] := 
  Module[{dims = Dimensions[data], numberOfInstances, 
  numberOfAttributes, leaves, slots, subleaves, nf},
  numberOfInstances = First[dims];
  numberOfAttributes = Last[dims] - 1;
 Table[(leaves = RandomSample[data, sizes];
        slots = Sort[slotList[[tree]]];
        subleaves = 
          Map[leaf \[Function] leaf[[slots]] -> leaf[[-1]], leaves];
        nf = Nearest[subleaves];
        {slots, nf}), {tree, Length[slotList]}]
  ]

 makeForest[data_, ntrees_, subdim_, sizes_] := 
   Module[{dims = Dimensions[data], numberOfAttributes},
     numberOfAttributes = Last[dims] - 1;
     makeForest[
        data, 
        Table[Sort[
          RandomSample[Range[numberOfAttributes], subdim]], {ntrees}
         ], 
        sizes
      ]
  ]

Disclaimer: This is not an implementation of the Random Forest Algorithm. Also, while I have on occasion used random florists, until today I had not heard of the Random Forest Algorithm.

I poked around a bit on the Net and learned that these take subsamples of data, subsampling the variables as well, and form decision trees for the subsetted subsamples. These then can be used to classify data points. One tests against each subsamples (tree) and returns the modal value.

Okay there is more that can be done with these, and I probably don't have even the above part correct. Whatever.

For some types of problem one can regard the "space" as having a vector of either continuous, discrete, or perhaps mixed values, with each datum then evaluating to some new value. That is to say, we pretend we have an unknown function f with f(x_1,...,x_n) = y The classifier goal is to find ` given the vector(x_1,...,x_n). The idea I have is to useNearestFunctionas a surrogate for a tree. One creates each tree using a random subset of data points, and a random subset of dimensions. Classification of data can then be done by a voting system:. For a particular input, see what values the forest (set ofNearestFunction`s) gives, and use the most common as the end result. or can use the several most common if a viable set of candidate values is desired.

I will assume the input data is of the form {vec1->val1,...,vecn->valn}. For purposes of terminology to describe the functions below:

rules is the master set ntrees is how many trees to create (size of forest) subdim is the dimension used for each tree sizes is the number of rules to use in each subtree

makeForest[rules_, ntrees_, subdim_, sizes_] := Module[
  {len = Length[rules], dim = Length[rules[[1, 1]]],
   slots, leaves, subleaves, nf},
  Table[
   leaves = rules[[RandomSample[Range[len], sizes]]];
   slots = Sort[RandomSample[Range[dim], subdim]];
   subleaves = Map[#[[1, slots]] -> #[[2]] &, leaves];
   nf = Nearest[subleaves];
   {slots, nf}
   , {ntrees}]
  ]

classify[vector_, forest_, n_] := Module[
  {vals, tree, subvec},
  vals = Flatten[Table[
     tree = forest[[j]];
     tree[[2]][vector[[tree[[1]]]]]
     , {j, Length[forest]}]];
  Reverse[SortBy[Tally[vals], Last][[-n ;; -1]]]
  ]

That's it.

Here is an example. I'll take a million points in the 8-cube {-1,1}^8. We give them ordinal values based on orthant (use base 2...) but before evaluating we "pollute" them first with Gaussian noise.

dim = 8;
size = 6;
pts = RandomReal[{-1, 1}, {10^size, dim}];

ruls = Map[# -> 
     FromDigits[(Sign[# + 
           RandomReal[NormalDistribution[0, .4], Length[#]]] + 1)/2, 
      2] &, pts];

Now we'll create a "forest" (I really should use a different term since I doubt this is what an actual random forest is, but...). I use 300 trees, each taking 4 (of the 8) vector positions, and each using 1000 elements from the full set, chosen at random.

Timing[rForest = makeForest[ruls, 300, 4, 1000];]

(* Out[80]= {4.740000, Null} *)

Now i create an random 8-vector with coordinates in {-1,1}, and also evaluate it according to our actual function.

SeedRandom[11111];
randvec = RandomReal[{-1, 1}, dim]
FromDigits[(Sign[randvec] + 1)/2, 2]

(* Out[82]= {-0.472529240338, 0.534376840246, 0.286266627232, \
0.690979660248, -0.663034535046, 0.715836737332, -0.389931777484, \
0.66884997164} *)

(* Out[83]= 117 *)

We run it through classify[].

Timing[candidates = classify[randvec, rForest, 5]]

(* Out[84]= {1.370000, {{117, 13}, {245, 10}, {116, 10}, {101, 10}, {53, 
   10}}} *)

While I doubt it will always give the "expected" result at the top of the list, it is encouraging to see that it might do so. Also the next ones do tend to share bits with 117.

Of course it might be unrealistic to use half of the full set of dimensions in the trees, so I really do now know how well this might scale to "real world" problems of high dimension. But I thought it might be worth showing this approach since it is simple in terms of code and others might have ideas on tweaking it for better accuracy/performance.