% This livescript will give an introduction to the topic of decision/regression trees.  On some level, regression trees are a bit similar to KNN classification/regression.  Recall that in KNN, the feature space is essentially broken up into potentially many non-overlapping regions (based on which K points are the nearest neighbors to the points in that region), and all test points falling into a given region are then classified or regressed to give the same value, which is generally the mode (in classification) or mean (in regression) of the y values for the K nearest neighbor test points that essentially define that region.
% Suppose now that you did not want to pre-specify how the regions of feature space ought to be determined (so, not necessarily by using the nearest neighbors), but still want to construct a method that will break up all of feature space into some number, call it J, of non-overlapping regions, then still classify test points by assigning the value that is the mode/mean of the y values of training points in that region.  In order to move toward developing this method, we should define some kind of cost function that we would like our method to minimize.  Of course, this could well just be the training RSS or classification error, as usual.  In the case of regression, then, we could use the following as our cost function:
% C(\mathbf{R})=\sum_{j=1}^J\sum_{i|x_i\in R_j}(y_i-\hat{y}_j)^2
% Within this formula, \mathbf{R} is the set of our J non-overlapping regions that cover feature space, and \hat{y}_j is the value that is applied to any point within specific region R_j.  Then this cost function just computes the RSS for each region and adds them all up.  The cost function is then fully determined by the regions that are defined.
% Note that, if you have N training datapoints and want to use J=N regions, then one can acheive 0 cost by using 1-NN!  In general, though, finding the set of regions that would provide a global minimum of this cost function for some J<N is seemingly very difficult, if only because these regions could in theory be whatever shapes you like (so long as they are non-overlapping).
% This is really the starting point for decision/regression trees.  They take the above problem and make a variety of simplifications to it, so that:
% a) a "solution" can be readily found and
% b) the solution provided is highly interpretable/actionable by practitioners (possibly no computer needed!)
% The first restriction that trees make is that all of the regions must be hyperrectangles whose roundaries are all parallel to the axes of feature space.  That is, each region can be described as existing within some contiguous band of x values in each dimension.  This helps not only to make finding a minimum easier (by restricting the possibilities, at least), but also helps to make the resulting tree more interpretable.  One need only specify the contiguous band (by giving a lower and upper limit) along each dimension for each region to fully specify the regions, and then any datapoint can be pretty easily placed into its correct region without any calculations being needed.
% However, the restriction above is not really sufficient to make the minimization problem "easy".  This is because there are still a lot of possible ways to break up the feature space into such hyperrectangles.  So, trees construct these regions by essentially abandon the idea of truly solving the minimization problem posed, and instead just try to find an approximate solution using what is called a greedy algorithm.  This algorithm attempts to "solve" the problem by making a sequence of finer and finer binary subdivisions of regions , each one of which is chosen to give the best resulting RSS at that time.  The algorithm is called greedy because it makes no attempt to look ahead into the future and ask how decisions now might impact things further down the line; instead each decision is made only to be the best one at that particular moment.
% To help clarify how this method works, let's consider an example that is lifted from logistic regression.  Here, we will work in two dimensional feature space, and do a regression task where there are only three numerical y values observed.  Do note that this could be performed as a classification task, but we will treat the y values as numerically meaningful and perform regression istead.

nr=100; % number of red, blue, and green points
nb=100;
ng=100;
sd=1; % a value for the standard deviation of the normals
xr=mvnrnd([-1,0],sd*eye(2),nr); % generate data, using different means to make the groups somewhat distinct
xb=mvnrnd([1,0.85],sd*eye(2),nb);
xg=mvnrnd([1,-0.85],sd*eye(2),ng);
scatter(xb(:,1),xb(:,2),'bo') % visualize the data
hold on;
scatter(xr(:,1),xr(:,2),'ro')
scatter(xg(:,1),xg(:,2),'go')
hold off
%find a bounding box for the data
data=[xr;xb;xg];
data=[data,[ones(nr,1);2*ones(nb,1);3*ones(ng,1)]]; %for y values red is 1, blue is 2, green is 3
x1min=min(data(:,1));
x1max=max(data(:,1));
x2min=min(data(:,2));
x2max=max(data(:,2));

% The tree generation algorithm will begin with the full bounding box that contains all of the datapoints within it in feature space.  This full region serves as the root of our tree, and contains all data points.  It then asks for this region: how can I make a binary split parallel to one of the two dimensions so as to minimize the total classification error that occurs given this split?  To write this in more mathematical terms, suppose I choose to split perpendicular to dimension d, which is either 1 or 2 for this 2d example.  Further suppose that I choose the split-point to be at location s_d, which then ends up defining two regions: the "right" region where x_d>s_d and the "left" region where x_d\leq s_d.  After making this choice, the two regions are fully described, and I can use these regions to define what \hat{y}_j is for each region by just taking the mean of the y values for the datapoints in that region (since this is a regression task).  Once that is determined, the current error can be determined:
% C(\mathbf{R})=\sum_{i|x_{id}\leq s_d}(y_i-\hat{y}_1)^2+\sum_{i|x_{id}> s_d}(y_i-\hat{y}_2)^2
% The section below lets you play with selecting a dimension and a cut point for the numerical example, and spits out the error given your choices, along with a picture of the two regions.
%

d=1;
xmin=(d==1)*x1min+(d==2)*x2min;
xmax=(d==1)*x1max+(d==2)*x2max;
delt=(xmax-xmin)/100;
s=-0.142844268222959;

%split data as chosen and compute error
inr1=data(:,d)<=s;
inr2=data(:,d)>s;
yhat1=mean(data(inr1,3));
yhat2=mean(data(inr2,3));
err1=(data(inr1,3)-yhat1)'*(data(inr1,3)-yhat1);
err2=(data(inr2,3)-yhat2)'*(data(inr2,3)-yhat2);
error=err1+err2

scatter(xb(:,1),xb(:,2),'bo') % visualize the data
hold on;
scatter(xr(:,1),xr(:,2),'ro')
scatter(xg(:,1),xg(:,2),'go')
if d==1
    plot([s,s],[x2min,x2max]);
else
    plot([x1min,x1max],[s,s]);
end
hold off

% Mess around with the values for d and s above to see if you can find a combination that leads to the lowest error.
%
% Given our setup, it is easy to show that there must be a (d,s_d) that minimizes this binary error.  Note that the only values of s_d that must be considered are the actual x_{id} values for our datapoints, as it is only at these particular values that the group membership of any points will actually change between the regions.  Since there are a finite number of those, and a finite number of dimensions, at least one of these finite combinations must be the global minimizer for this error.  The code below determines this for our example:

m=4; %define minimum number of points for a region
nt=nr+nb+ng;
x1min=min(data(:,1))-.1;
x1max=max(data(:,1))+.1;
x2min=min(data(:,2))-.1;
x2max=max(data(:,2))+.1;
yhat=mean(data(:,3));
err=(data(:,3)-yhat)'*(data(:,3)-yhat);
regions=[x1min,x1max,x2min,x2max,err,0,nt,true]; %track current regions by bounding box, error before split, index of parent region, number of points, and flag on whether split or not yet

box=[x1min,x1max;x2min,x2max];
errs1=9*nt*ones(nt,2);
errs2=errs1;
n1s=errs1;
n2s=errs1;
for d=1:2
    sortsubdata=sortrows(data,d);
    for i=m:nt-m
        s=(sortsubdata(i,d)+sortsubdata(i+1,d))/2;
        inr1=data(:,d)<=s;
        n1s(i,d)=sum(inr1);
        inr2=data(:,d)>s;
        n2s(i,d)=sum(inr2);
        yhat1=mean(data(inr1,3));
        yhat2=mean(data(inr2,3));
        errs1(i,d)=(data(inr1,3)-yhat1)'*(data(inr1,3)-yhat1);
        errs2(i,d)=(data(inr2,3)-yhat2)'*(data(inr2,3)-yhat2);
    end
end
errors=errs1+errs2;
[mins,locs]=min(errors);
[err0,d0]=min(mins)
sortsubdata=sortrows(data,d0);
ind0=locs(d0);
s0=(sortsubdata(ind0,d0)+sortsubdata(ind0+1,d0))/2
if err0<err %only split if error decreased.  Need to check given minimum size stuff
    if d0==1
        bb=[x1min,s0,x2min,x2max,errs1(ind0,d0),1,n1s(ind0,d0),2*m<=n1s(ind0,d0);s0,x1max,x2min,x2max,errs2(ind0,d0),1,n2s(ind0,d0),2*m<=n2s(ind0,d0)];
    else
        bb=[x1min,x1max,x2min,s0,errs1(ind0,d0),1,n1s(ind0,d0),2*m<=n1s(ind0,d0);x1min,x1max,s0,x2max,errs2(ind0,d0),1,n2s(ind0,d0),2*m<=n2s(ind0,d0)];
    end
    regions=[regions;bb];
end
regions(1,end)=false;

scatter(xb(:,1),xb(:,2),'bo') % visualize the data
hold on;
scatter(xr(:,1),xr(:,2),'ro')
scatter(xg(:,1),xg(:,2),'go')
[nreg,~]=size(regions);
for k=1:nreg
    rectangle('Position',[regions(k,1),regions(k,3),regions(k,2)-regions(k,1),regions(k,4)-regions(k,3)])
end
hold off

%
% At this point, we have made a single binary split of the data so that what was once one region of feature space, the trunk of the tree, has been divided into two branches.  We then need to ask ourselves if we should stop there, or proceed in some way.  In general, having only two branches doesn't constitute much of a tree, so we will continue onward.  The way we do so follows directly from what we've already done.  That is, for each of these two branches, we will again (in principle) perform a binary cut to separate each into two regions, leaving a total of four regions in the end.  Run the code below one time to do this:

tosplit=regions(:,end)==1;
ntosplit=sum(tosplit);
if ntosplit>0
    splitters=find(tosplit);
    for k=1:ntosplit
        region=regions(splitters(k),:);
        err=region(5);
        x1min=region(1);
        x1max=region(2);
        x2min=region(3);
        x2max=region(4);
        box=[x1min,x1max;x2min,x2max];
        inregion=(x1min<data(:,1))&(data(:,1)<x1max)&(x2min<data(:,2))&(data(:,2)<x2max);
        subdata=data(inregion,:);
        nt=sum(inregion);

        errs1=9*nt*ones(nt,2);
        errs2=errs1;
        n1s=errs1;
        n2s=errs1;
        for d=1:2
            sortsubdata=sortrows(subdata,d);
            for i=m:nt-m
                s=(sortsubdata(i,d)+sortsubdata(i+1,d))/2;
                inr1=subdata(:,d)<=s;
                n1s(i,d)=sum(inr1);
                inr2=subdata(:,d)>s;
                n2s(i,d)=sum(inr2);
                yhat1=mean(subdata(inr1,3));
                yhat2=mean(subdata(inr2,3));
                errs1(i,d)=(subdata(inr1,3)-yhat1)'*(subdata(inr1,3)-yhat1);
                errs2(i,d)=(subdata(inr2,3)-yhat2)'*(subdata(inr2,3)-yhat2);
            end
        end
        errors=errs1+errs2;
        [mins,locs]=min(errors);
        [err0,d0]=min(mins);
        sortsubdata=sortrows(subdata,d0);
        ind0=locs(d0);
        s0=(sortsubdata(ind0,d0)+sortsubdata(ind0+1,d0))/2;
        if err0<err
            if d0==1
                bb=[x1min,s0,x2min,x2max,errs1(ind0,d0),splitters(k),n1s(ind0,d0),2*m<=n1s(ind0,d0);s0,x1max,x2min,x2max,errs2(ind0,d0),splitters(k),n2s(ind0,d0),2*m<=n2s(ind0,d0)];
            else
                bb=[x1min,x1max,x2min,s0,errs1(ind0,d0),splitters(k),n1s(ind0,d0),2*m<=n1s(ind0,d0);x1min,x1max,s0,x2max,errs2(ind0,d0),splitters(k),n2s(ind0,d0),2*m<=n2s(ind0,d0)];
            end
            regions=[regions;bb];
        end
        regions(splitters(k),end)=false;
    end
else
    disp('no more splits are possible')
end


scatter(xb(:,1),xb(:,2),'bo') % visualize the data
hold on;
scatter(xr(:,1),xr(:,2),'ro')
scatter(xg(:,1),xg(:,2),'go')
[nreg,~]=size(regions);
for k=1:nreg
    rectangle('Position',[regions(k,1),regions(k,3),regions(k,2)-regions(k,1),regions(k,4)-regions(k,3)])
end
hold off
%below we separate branch regions from leaves, and compute for each branch
%how many leaves eventually arise from it and total error of those leaves
branches=unique(regions(:,6));
branches=branches(2:end); %remove the extra "0" region, which doesn't acutally exist but is listed as the parent of the full dataset region
[nbranch,~]=size(branches);
nleaves=nreg-nbranch
leaves=setdiff((1:nreg)',branches);
nleavesbybranch=zeros(nreg,1);
qleavesbybranch=nleavesbybranch;
for k=1:nleaves
    leaf=leaves(k);
    pbranch=regions(leaf,6);
    nleavesbybranch(pbranch)=nleavesbybranch(pbranch)+1;
    qleavesbybranch(pbranch)=qleavesbybranch(pbranch)+regions(leaf,5);
end
for k=nbranch:-1:2
    branch=branches(k);
    pbranch=regions(branch,6);
    nleavesbybranch(pbranch)=nleavesbybranch(pbranch)+nleavesbybranch(branch);
    qleavesbybranch(pbranch)=qleavesbybranch(pbranch)+qleavesbybranch(branch);
end

% Run the code above as many times as you like to make more and more splits, until it tells you no more splits are possible.  The resulting tree has numerous branches (all the various binary splittings we made along the way) and a number of leaves, which are just whatever regions remain at the end.
%
% It is nice to see this look like an actual tree.  I don't know how to code this, so let's use Matlab's built in (regression) tree constructor:

tree=fitrtree(data(:,1:2),data(:,3));
view(tree,'Mode','graph');
leaves=sum(~tree.IsBranchNode)
histogram(tree.NodeSize(~tree.IsBranchNode)) %visualize the number of points per leaf by a histogram

%
% After performing the method above until we are left with regions (leaves) with only a small number of points each, we will likely have a pretty big tree that has plausibly performed some overfitting.  After all, in the example we are working through here, there should probably not be too many regions, since we really only have three possible y values.  But note that the tree created by Matlab has a lot of leaves.  This is where we can go back and prune our tree.  Pruning the tree involves essentially collapsing certain branches back down, such that all the subregions coming off of some branch point are placed back into a single region that is then just a leaf on the tree.  Note that the Matlab tree viewer shows a "pruning level" tab, where you can choose the level you want to prune the tree by; mess around with this to see it collapse branches down to leaves as you increase the pruning level.
% But, how is this accomplished?  In general, the goal of pruning is the balance accuracy vs complexity; this is essentially the bias-variance tradeoff we've discussed many times before.  One specific method for doing this is called cost-complexity pruning.  To discuss this method, let's introduce a little notation.  Denote our original, unpruned tree as T_0.  Define a subtree of T_0 by T\subset T_0, which is just any tree that might result from the pruning of our original tree (with pruning defined above).  Let the number of leaves on any tree T be denoted by |T|.  Finally, let the RSS within leaf \ell of tree T be Q_\ell(T):
% Q_\ell(T)=\sum_{i|x_i\in\ell}(y_i-\hat{y}_\ell)^2
% Then cost-complexity pruning defines the following object function:
% C_\alpha(T)=\sum_{\ell=1}^{|T|}Q_\ell(T)+\alpha|T|,
% where \alpha\geq0 is a user chosen parameter.  For a given choice of \alpha, the optimal pruned tree according to cost-complexity pruning is that tree that is the smallest one that minimizes the function C_\alpha(T).  So, cost complexity pruning is trying to balance the overall RSS of a given pruned tree with the size of that tree as measured by the number of leaves it possesses, and the parameter \alpha allows the user to move between only a small preference for smaller trees (low \alpha) to a very high preference for small trees (high \alpha).
% While the basic idea and setup of cost-complexity pruning might (hopefully) be clear, how to actually find this optimal subtree for a given \alpha may seem less so.  At least one can note that such an optimal tree must exist: any given original tree has a finite number of subtrees, so there must exist some global minimum of the cost function that at least one subtree  yields.  Of course, there could perhaps be ties between trees, in which case the requirement for the "smallest" tree that minimizes the objective may help decide among those.  It actually turns out that this smallest objective-minimizing subtree is unique, and there is an algorithm for finding it called "weakest link pruning" (see class notes for how this works).
%