Daniel@0: function CPD = learn_params(CPD, fam, data, ns, cnodes, varargin)
Daniel@0: % LEARN_PARAMS Construct classification/regression tree given complete data
Daniel@0: % CPD = learn_params(CPD, fam, data, ns, cnodes)
Daniel@0: %
Daniel@0: % fam(i) is the node id of the i-th node in the family of nodes, self node is the last one
Daniel@0: % data(i,m) is the value of node i in case m (can be cell array).
Daniel@0: % ns(i) is the node size for the i-th node in the whold bnet
Daniel@0: % cnodes(i) is the node id for the i-th continuous node in the whole bnet
Daniel@0: %  
Daniel@0: % The following optional arguments can be specified in the form of name/value pairs:
Daniel@0: % stop_cases: for early stop (pruning). A node is not split if it has less than k cases. default is 0.
Daniel@0: % min_gain: for early stop (pruning). 
Daniel@0: %     For discrete output: A node is not split when the gain of best split is less than min_gain. default is 0.  
Daniel@0: %     For continuous (cts) outpt: A node is not split when the gain of best split is less than min_gain*score(root) 
Daniel@0: %                                 (we denote it cts_min_gain). default is 0.006
Daniel@0: % %%%%%%%%%%%%%%%%%%%Struction definition of dtree_CPD.tree%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Daniel@0: % tree.num_node               the last position in tree.nodes array for adding new nodes,
Daniel@0: %                             it is not always same to number of nodes in a tree, because some position in the 
Daniel@0: %                             tree.nodes array can be set to unused (e.g. in tree pruning)  
Daniel@0: % tree.nodes is the array of nodes in the tree plus some unused nodes.
Daniel@0: % tree.nodes(1) is the root for the tree.
Daniel@0: %
Daniel@0: % Below is the attributes for each node
Daniel@0: % tree.nodes(i).used;     % flag this node is used (0 means node not used, it can be removed from tree to save memory)
Daniel@0: % tree.nodes(i).is_leaf;  % if 1 means this node is a leaf, if 0 not a leaf.
Daniel@0: % tree.nodes(i).children; % children(i) is the node number in tree.nodes array for the i-th child node
Daniel@0: % tree.nodes(i).split_id; % the attribute id used to split this node
Daniel@0: % tree.nodes(i).split_threshhold; % the threshhold for continuous attribute to split this node
Daniel@0: % %%%%%attributes specially for classification tree (discrete output)
Daniel@0: % tree.nodes(i).probs     % probs(i) is the prob for i-th value of class node 
Daniel@0: %                         % For three output class, the probs = [0.9 0.1 0.0] means the probability of 
Daniel@0: %                         % class 1 is 0.9, for class 2 is 0.1, for class 3 is 0.0.
Daniel@0: % %%%%%attributes specially for regression tree (continuous output)                          
Daniel@0: % tree.nodes(i).mean      % mean output value for this node
Daniel@0: % tree.nodes(i).std       % standard deviation for output values in this node
Daniel@0: %
Daniel@0: % Author: yimin.zhang@intel.com
Daniel@0: % Last updated: Jan. 19, 2002
Daniel@0: 
Daniel@0: % Want list:
Daniel@0: % (1) more efficient for cts attributes: get the values of cts attributes at first (the begining of build_tree function), then doing bi_search in finding threshhold
Daniel@0: % (2) pruning classification tree using Pessimistic Error Pruning
Daniel@0: % (3) bi_search for strings (used for transform data to BNT format)
Daniel@0: 
Daniel@0: global tree %tree must be global so that it can be accessed in recursive slitting function
Daniel@0: global cts_min_gain
Daniel@0: tree=[]; % clear the tree
Daniel@0: tree.num_node=0;
Daniel@0: cts_min_gain=0;
Daniel@0: 
Daniel@0: stop_cases=0;
Daniel@0: min_gain=0;
Daniel@0: 
Daniel@0: args = varargin;
Daniel@0: nargs = length(args);
Daniel@0: if (nargs>0)
Daniel@0:   if isstr(args{1})
Daniel@0:     for i=1:2:nargs
Daniel@0:       switch args{i},
Daniel@0:         case 'stop_cases', stop_cases = args{i+1};   
Daniel@0:         case 'min_gain', min_gain = args{i+1};
Daniel@0:       end
Daniel@0:     end
Daniel@0:   else
Daniel@0:     error(['error in input parameters']);
Daniel@0:   end
Daniel@0: end
Daniel@0: 
Daniel@0: if iscell(data)
Daniel@0:   local_data = cell2num(data(fam,:));
Daniel@0: else
Daniel@0:   local_data = data(fam, :);
Daniel@0: end
Daniel@0: %counts = compute_counts(local_data, CPD.sizes);
Daniel@0: %CPD.CPT = mk_stochastic(counts + CPD.prior); % bug fix 11/5/01
Daniel@0: node_types = zeros(1,size(ns,2)); %all nodes are disrete
Daniel@0: node_types(cnodes)=1;
Daniel@0: %make the data be BNT compliant (values for discrete nodes are from 1-n, here n is the node size)
Daniel@0: %trans_data=transform_data(local_data,'tmp.dat',[]); %here no cts nodes
Daniel@0: 
Daniel@0: build_dtree (CPD, local_data, ns(fam), node_types(fam),stop_cases,min_gain);
Daniel@0: %CPD.tree=copy_tree(tree);
Daniel@0: CPD.tree=tree; %copy the tree constructed to CPD
Daniel@0: 
Daniel@0: 
Daniel@0: function new_tree = copy_tree(tree)
Daniel@0: % copy the tree to new_tree
Daniel@0: new_tree.num_node=tree.num_node;
Daniel@0: new_tree.root = tree.root;
Daniel@0: for i=1:tree.num_node
Daniel@0:   new_tree.nodes(i)=tree.nodes(i);
Daniel@0: end
Daniel@0: 
Daniel@0: 
Daniel@0: function build_dtree (CPD, fam_ev, node_sizes, node_types,stop_cases,min_gain)
Daniel@0: global tree
Daniel@0: global cts_min_gain
Daniel@0: 
Daniel@0: tree.num_node=0; %the current number of nodes in the tree
Daniel@0: tree.root=1;
Daniel@0: 
Daniel@0: T = 1:size(fam_ev,2) ; %all cases
Daniel@0: candidate_attrs = 1:(size(node_sizes,2)-1); %all attributes
Daniel@0: node_id=1;  %the root node
Daniel@0: lastnode=size(node_sizes,2); %the last element in all nodes is the dependent variable (category node)
Daniel@0: num_cat=node_sizes(lastnode);
Daniel@0: 
Daniel@0: % get minimum gain for cts output (used in stop splitting)
Daniel@0: if (node_types(size(fam_ev,1))==1) %cts output
Daniel@0:   N = size(fam_ev,2);
Daniel@0:   output_id = size(fam_ev,1);
Daniel@0:   cases_T = fam_ev(output_id,:); %get all the output value for cases T
Daniel@0:   std_T = std(cases_T);
Daniel@0:   avg_y_T = mean(cases_T);
Daniel@0:   sqr_T = cases_T - avg_y_T;
Daniel@0:   cts_min_gain = min_gain*(sum(sqr_T.*sqr_T)/N);  % min_gain * (R(root) = 1/N * SUM(y-avg_y)^2)
Daniel@0: end  
Daniel@0: 
Daniel@0: split_dtree (CPD, fam_ev, node_sizes, node_types, stop_cases,min_gain, T, candidate_attrs, num_cat);
Daniel@0:   
Daniel@0: 
Daniel@0: 
Daniel@0: % pruning method
Daniel@0: % (1) Restrictions on minimum node size: A node is not split if it has smaller than k cases.
Daniel@0: % (2) Threshholds on impurity: a threshhold is imposed on the splitting test score. Threshhold can be 
Daniel@0: % imposed on local goodness measure (the gain_ratio of a node) or global goodness.
Daniel@0: % (3) Mininum Error Pruning (MEP), (no need pruning set)
Daniel@0: %     Prune if static error<=backed-up error
Daniel@0: %      Static error at node v: e(v) = (Nc + 1)/(N+k) (laplace estimate, prior for each class equal) 
Daniel@0: %        here N is # of all examples, Nc is # of majority class examples, k is number of classes 
Daniel@0: %      Backed-up error at node v: (Ti is the i-th subtree root)
Daniel@0: %         E(T) = Sum_1_to_n(pi*e(Ti))
Daniel@0: % (4) Pessimistic Error Pruning (PEP), used in Quilan C4.5 (no need pruning set, efficient because of pruning top-down)
Daniel@0: %       Probability of error (apparent error rate)
Daniel@0: %           q = (N-Nc+0.5)/N
Daniel@0: %         where N=#examples, Nc=#examples in majority class
Daniel@0: %     Error of a node v (if pruned)  q(v)= (Nv- Nc,v + 0.5)/Nv
Daniel@0: %     Error of a subtree   q(T)= Sum_of_l_leaves(Nl - Nc,l + 0.5)/Sum_of_l_leaves(Nl)
Daniel@0: %     Prune if q(v)<=q(T)
Daniel@0: % 
Daniel@0: % Implementation statuts:
Daniel@0: % (1)(2) has been implemented as the input parameters of learn_params.
Daniel@0: % (4) is implemented in this function
Daniel@0: function pruning(fam_ev,node_sizes,node_types)
Daniel@0: % PRUNING prune the constructed tree using PEP
Daniel@0: % pruning(fam_ev,node_sizes,node_types)
Daniel@0: %
Daniel@0: % fam_ev(i,j)  is the value of attribute i in j-th training cases (for whole tree), the last row is for the class label (self_ev)
Daniel@0: % node_sizes(i) is the node size for the i-th node in the family
Daniel@0: % node_types(i) is the node type for the i-th node in the family, 0 for disrete node, 1 for continous node
Daniel@0: % the global parameter 'tree' is for storing the input tree and the pruned tree
Daniel@0: 
Daniel@0: 
Daniel@0: function split_T = split_cases(fam_ev,node_sizes,node_types,T,node_i, threshhold)
Daniel@0: % SPLIT_CASES split the cases T according to values of node_i in the family
Daniel@0: % split_T = split_cases(fam_ev,node_sizes,node_types,T,node_i)
Daniel@0: %
Daniel@0: % fam_ev(i,j)  is the value of attribute i in j-th training cases (for whole tree), the last row is for the class label (self_ev)
Daniel@0: % node_sizes(i) is the node size for the i-th node in the family
Daniel@0: % node_types(i) is the node type for the i-th node in the family, 0 for disrete node, 1 for continous node
Daniel@0: % node_i is the attribute we need to split
Daniel@0: 
Daniel@0: if (node_types(node_i)==0) %discrete attribute
Daniel@0:   %init the subsets of T
Daniel@0:   split_T = cell(1,node_sizes(node_i)); %T will be separated into |node_size of i| subsets according to different values of node i
Daniel@0:   for i=1:node_sizes(node_i)   % here we assume that the value of an attribute is 1:node_size
Daniel@0:     split_T{i}=zeros(1,0);
Daniel@0:   end
Daniel@0: 
Daniel@0:   size_t = size(T,2);
Daniel@0:   for i=1:size_t
Daniel@0:     case_id = T(i);
Daniel@0:     %put this case into one subset of split_T according to its value for node_i
Daniel@0:     value = fam_ev(node_i,case_id); 
Daniel@0:     pos = size(split_T{value},2)+1;
Daniel@0:     split_T{value}(pos)=case_id;  % here assumes the value of an attribute is 1:node_size 
Daniel@0:   end
Daniel@0: else %continuous attribute
Daniel@0:   %init the subsets of T
Daniel@0:   split_T = cell(1,2); %T will be separated into 2 subsets (<=threshhold) (>threshhold)
Daniel@0:   for i=1:2   
Daniel@0:     split_T{i}=zeros(1,0);
Daniel@0:   end
Daniel@0: 
Daniel@0:   size_t = size(T,2);
Daniel@0:   for i=1:size_t
Daniel@0:     case_id = T(i);
Daniel@0:     %put this case into one subset of split_T according to its value for node_i
Daniel@0:     value = fam_ev(node_i,case_id); 
Daniel@0:     subset_num=1;
Daniel@0:     if (value>threshhold)
Daniel@0:       subset_num=2;
Daniel@0:     end  
Daniel@0:     pos = size(split_T{subset_num},2)+1;
Daniel@0:     split_T{subset_num}(pos)=case_id;  
Daniel@0:   end
Daniel@0: end
Daniel@0: 
Daniel@0: 
Daniel@0:   
Daniel@0: function new_node = split_dtree (CPD, fam_ev, node_sizes, node_types, stop_cases, min_gain, T, candidate_attrs, num_cat)
Daniel@0: % SPLIT_TREE Split the tree at node node_id with cases T (actually it is just indexes to family evidences).
Daniel@0: % new_node = split_dtree (fam_ev, node_sizes, node_types, T, node_id, num_cat, method)
Daniel@0: %
Daniel@0: % fam_ev(i,j)  is the value of attribute i in j-th training cases (for whole tree), the last row is for the class label (self_ev)
Daniel@0: % node_sizes{i} is the node size for the i-th node in the family
Daniel@0: % node_types{i} is the node type for the i-th node in the family, 0 for disrete node, 1 for continous node
Daniel@0: % stop_cases is the threshold of number of cases to stop slitting
Daniel@0: % min_gain is the minimum gain need to split a node
Daniel@0: % T(i) is the index of i-th cases in current decision tree node, we need split it further
Daniel@0: % candidate_attrs(i) the node id for the i-th attribute that still need to be considered as split attribute 
Daniel@0: %%%%% node_id is the index of current node considered for a split
Daniel@0: % num_cat is the number of output categories for the decision tree
Daniel@0: % output:
Daniel@0: % new_node is the new node created
Daniel@0: global tree
Daniel@0: global cts_min_gain
Daniel@0: 
Daniel@0: size_fam = size(fam_ev,1);            %number of family size
Daniel@0: output_type = node_types(size_fam);   %the type of output for the tree (0 is discrete, 1 is continuous)
Daniel@0: size_attrs = size(candidate_attrs,2); %number of candidate attributes
Daniel@0: size_t = size(T,2);                   %number of training cases in this tree node
Daniel@0: 
Daniel@0: %(1)computeFrequenceyForEachClass(T)
Daniel@0: if (output_type==0) %discrete output
Daniel@0:   class_freqs = zeros(1,num_cat);
Daniel@0:   for i=1:size_t
Daniel@0:     case_id = T(i);
Daniel@0:     case_class = fam_ev(size_fam,case_id); %get the class label for this case
Daniel@0:     class_freqs(case_class)=class_freqs(case_class)+1;
Daniel@0:   end
Daniel@0: else  %cts output
Daniel@0:   N = size(fam_ev,2);
Daniel@0:   cases_T = fam_ev(size(fam_ev,1),T); %get the output value for cases T
Daniel@0:   std_T = std(cases_T);
Daniel@0: end
Daniel@0: 
Daniel@0: %(2) if OneClass (for discrete output) or same output value (for cts output) or Class With #examples < stop_cases
Daniel@0: %         return a leaf;
Daniel@0: %    create a decision node N;
Daniel@0: 
Daniel@0: % get majority class in this node
Daniel@0: if (output_type == 0)
Daniel@0:   top1_class = 0;       %the class with the largest number of cases
Daniel@0:   top1_class_cases = 0; %the number of cases in top1_class
Daniel@0:   [top1_class_cases,top1_class]=max(class_freqs);
Daniel@0: end
Daniel@0:   
Daniel@0: if (size_t==0)     %impossble
Daniel@0:   new_node=-1;
Daniel@0:   fprintf('Fatal error: please contact the author. \n');
Daniel@0:   return;
Daniel@0: end
Daniel@0: 
Daniel@0: % stop splitting if needed
Daniel@0:   %for discrete output: one class 
Daniel@0:   %for cts output, all output value in cases are same
Daniel@0:   %cases too little
Daniel@0: if ( (output_type==0 & top1_class_cases == size_t) | (output_type==1 & std_T == 0) | (size_t < stop_cases))             
Daniel@0:   %create one new leaf node
Daniel@0:   tree.num_node=tree.num_node+1;
Daniel@0:   tree.nodes(tree.num_node).used=1; %flag this node is used (0 means node not used, it will be removed from tree at last to save memory)
Daniel@0:   tree.nodes(tree.num_node).is_leaf=1;
Daniel@0:   tree.nodes(tree.num_node).children=[];
Daniel@0:   tree.nodes(tree.num_node).split_id=0;  %the attribute(parent) id to split this tree node
Daniel@0:   tree.nodes(tree.num_node).split_threshhold=0;  
Daniel@0:   if (output_type==0)
Daniel@0:     tree.nodes(tree.num_node).probs=class_freqs/size_t; %the prob for each value of class node 
Daniel@0: 
Daniel@0:     %  tree.nodes(tree.num_node).probs=zeros(1,num_cat); %the prob for each value of class node 
Daniel@0:     %  tree.nodes(tree.num_node).probs(top1_class)=1; %use the majority class of parent node, like for binary class, 
Daniel@0:                                                    %and majority is class 2, then the CPT is [0 1]
Daniel@0:                                                    %we may need to use prior to do smoothing, to get [0.001 0.999]
Daniel@0:     tree.nodes(tree.num_node).error.self_error=1-top1_class_cases/size_t; %the classfication error in this tree node when use default class
Daniel@0:     tree.nodes(tree.num_node).error.all_error=1-top1_class_cases/size_t;  %no total classfication error in this tree node and its subtree
Daniel@0:     tree.nodes(tree.num_node).error.all_error_num=size_t - top1_class_cases;
Daniel@0:     fprintf('Create leaf node(onecla) %d. Class %d Cases %d Error %d \n',tree.num_node, top1_class, size_t, size_t - top1_class_cases );
Daniel@0:   else
Daniel@0:     avg_y_T = mean(cases_T);
Daniel@0:     tree.nodes(tree.num_node).mean = avg_y_T; 
Daniel@0:     tree.nodes(tree.num_node).std = std_T;
Daniel@0:     fprintf('Create leaf node(samevalue) %d. Mean %8.4f Std %8.4f Cases %d \n',tree.num_node, avg_y_T, std_T, size_t);
Daniel@0:   end  
Daniel@0:   new_node = tree.num_node;
Daniel@0:   return;
Daniel@0: end
Daniel@0:     
Daniel@0: %create one new node
Daniel@0: tree.num_node=tree.num_node+1;
Daniel@0: tree.nodes(tree.num_node).used=1; %flag this node is used (0 means node not used, it will be removed from tree at last to save memory)
Daniel@0: tree.nodes(tree.num_node).is_leaf=1;
Daniel@0: tree.nodes(tree.num_node).children=[];
Daniel@0: tree.nodes(tree.num_node).split_id=0;
Daniel@0: tree.nodes(tree.num_node).split_threshhold=0;  
Daniel@0: if (output_type==0)
Daniel@0:   tree.nodes(tree.num_node).error.self_error=1-top1_class_cases/size_t; 
Daniel@0:   tree.nodes(tree.num_node).error.all_error=0;
Daniel@0:   tree.nodes(tree.num_node).error.all_error_num=0;
Daniel@0: else
Daniel@0:   avg_y_T = mean(cases_T);
Daniel@0:   tree.nodes(tree.num_node).mean = avg_y_T; 
Daniel@0:   tree.nodes(tree.num_node).std = std_T;
Daniel@0: end
Daniel@0: new_node = tree.num_node;
Daniel@0: 
Daniel@0: %Stop splitting if no attributes left in this node
Daniel@0: if (size_attrs==0) 
Daniel@0:   if (output_type==0)
Daniel@0:     tree.nodes(tree.num_node).probs=class_freqs/size_t; %the prob for each value of class node 
Daniel@0:     tree.nodes(tree.num_node).error.all_error=1-top1_class_cases/size_t;  
Daniel@0:     tree.nodes(tree.num_node).error.all_error_num=size_t - top1_class_cases;
Daniel@0:     fprintf('Create leaf node(noattr) %d. Class %d Cases %d Error %d \n',tree.num_node, top1_class, size_t, size_t - top1_class_cases );
Daniel@0:   else
Daniel@0:     fprintf('Create leaf node(noattr) %d. Mean %8.4f Std %8.4f Cases %d \n',tree.num_node, avg_y_T, std_T, size_t);
Daniel@0:   end
Daniel@0:   return;
Daniel@0: end
Daniel@0:       
Daniel@0:   
Daniel@0: %(3) for each attribute A
Daniel@0: %        ComputeGain(A);
Daniel@0: max_gain=0;  %the max gain score (for discrete information gain or gain ration, for cts node the R(T))
Daniel@0: best_attr=0;  %the attribute with the max_gain
Daniel@0: best_split = []; %the split of T according to the value of best_attr
Daniel@0: cur_best_threshhold = 0; %the threshhold for split continuous attribute
Daniel@0: best_threshhold=0;
Daniel@0: 
Daniel@0: % compute Info(T) (for discrete output)
Daniel@0: if (output_type == 0)
Daniel@0:   class_split_T = split_cases(fam_ev,node_sizes,node_types,T,size(fam_ev,1),0); %split cases according to class
Daniel@0:   info_T = compute_info (fam_ev, T, class_split_T);
Daniel@0: else % compute R(T) (for cts output)
Daniel@0: %  N = size(fam_ev,2);
Daniel@0: %  cases_T = fam_ev(size(fam_ev,1),T); %get the output value for cases T
Daniel@0: %  std_T = std(cases_T);
Daniel@0: %  avg_y_T = mean(cases_T);
Daniel@0:   sqr_T = cases_T - avg_y_T;
Daniel@0:   R_T = sum(sqr_T.*sqr_T)/N;  % get R(T) = 1/N * SUM(y-avg_y)^2
Daniel@0:   info_T = R_T;
Daniel@0: end
Daniel@0: 
Daniel@0: for i=1:(size_fam-1)
Daniel@0:   if (myismember(i,candidate_attrs))  %if this attribute still in the candidate attribute set
Daniel@0:     if (node_types(i)==0) %discrete attibute
Daniel@0:       split_T = split_cases(fam_ev,node_sizes,node_types,T,i,0); %split cases according to value of attribute i
Daniel@0:       % For cts output, we compute the least square gain.
Daniel@0:       % For discrete output, we compute gain ratio
Daniel@0:       cur_gain = compute_gain(fam_ev,node_sizes,node_types,T,info_T,i,split_T,0,output_type); %gain ratio
Daniel@0:     else %cts attribute
Daniel@0:       %get the values of this attribute
Daniel@0:       ev = fam_ev(:,T);
Daniel@0:       values = ev(i,:);
Daniel@0:       sort_v = sort(values); 
Daniel@0:         %remove the duplicate values in sort_v
Daniel@0:       v_set = unique(sort_v);  
Daniel@0:       best_gain = 0;
Daniel@0:       best_threshhold = 0;
Daniel@0:       best_split1 = [];
Daniel@0:       
Daniel@0:       %find the best split for this cts attribute
Daniel@0:       % see "Quilan 96: Improved Use of Continuous Attributes in C4.5"
Daniel@0:       for j=1:(size(v_set,2)-1)
Daniel@0:         mid_v = (v_set(j)+v_set(j+1))/2; 
Daniel@0:         split_T = split_cases(fam_ev,node_sizes,node_types,T,i,mid_v); %split cases according to value of attribute i (<=mid_v)
Daniel@0:         % For cts output, we compute the least square gain.
Daniel@0:         % For discrete output, we use Quilan 96: use information gain instead of gain ratio to select threshhold
Daniel@0:         cur_gain = compute_gain(fam_ev,node_sizes,node_types,T,info_T,i,split_T,1,output_type); 
Daniel@0:         %if (i==6)
Daniel@0:         %  fprintf('gain %8.5f threshhold %6.3f spliting %d\n', cur_gain, mid_v, size(split_T{1},2));
Daniel@0:         %end
Daniel@0: 
Daniel@0:         if (best_gain < cur_gain)
Daniel@0:           best_gain = cur_gain;
Daniel@0:           best_threshhold = mid_v;
Daniel@0:           %best_split1 = split_T;     %here we need to copy array, not good!!! (maybe we can compute after we get best_attr
Daniel@0:         end
Daniel@0:       end
Daniel@0:       %recalculate the gain_ratio of the best_threshhold
Daniel@0:       split_T = split_cases(fam_ev,node_sizes,node_types,T,i,best_threshhold);
Daniel@0:       best_gain = compute_gain(fam_ev,node_sizes,node_types,T,info_T,i,split_T,0,output_type); %gain_ratio
Daniel@0:       if (output_type==0) %for discrete output
Daniel@0:         cur_gain = best_gain-log2(size(v_set,2)-1)/size_t; % Quilan 96: use the gain_ratio-log2(N-1)/|D| as the gain of this attr
Daniel@0:       else                %for cts output
Daniel@0:         cur_gain = best_gain;
Daniel@0:       end
Daniel@0:     end
Daniel@0:     
Daniel@0:     if (max_gain < cur_gain)
Daniel@0:       max_gain = cur_gain;
Daniel@0:       best_attr = i;
Daniel@0:       cur_best_threshhold=best_threshhold;  %save the threshhold
Daniel@0:       %best_split = split_T;        %here we need to copy array, not good!!! So we will recalculate in below line 313
Daniel@0:     end
Daniel@0:   end
Daniel@0: end
Daniel@0: 
Daniel@0: % stop splitting if gain is too small
Daniel@0: if (max_gain==0 | (output_type==0 & max_gain < min_gain) | (output_type==1 & max_gain < cts_min_gain)) 
Daniel@0:   if (output_type==0)
Daniel@0:     tree.nodes(tree.num_node).probs=class_freqs/size_t; %the prob for each value of class node 
Daniel@0:     tree.nodes(tree.num_node).error.all_error=1-top1_class_cases/size_t;  
Daniel@0:     tree.nodes(tree.num_node).error.all_error_num=size_t - top1_class_cases;
Daniel@0:     fprintf('Create leaf node(nogain) %d. Class %d Cases %d Error %d \n',tree.num_node, top1_class, size_t, size_t - top1_class_cases );
Daniel@0:   else
Daniel@0:     fprintf('Create leaf node(nogain) %d. Mean %8.4f Std %8.4f Cases %d \n',tree.num_node, avg_y_T, std_T, size_t);
Daniel@0:   end
Daniel@0:   return;
Daniel@0: end
Daniel@0: 
Daniel@0: %get the split of cases according to the best split attribute
Daniel@0: if (node_types(best_attr)==0) %discrete attibute
Daniel@0:   best_split = split_cases(fam_ev,node_sizes,node_types,T,best_attr,0);  
Daniel@0: else  
Daniel@0:   best_split = split_cases(fam_ev,node_sizes,node_types,T,best_attr,cur_best_threshhold);
Daniel@0: end
Daniel@0:   
Daniel@0: %(4) best_attr = AttributeWithBestGain;
Daniel@0: %(5) if best_attr is continuous             ???? why need this? maybe the value in the decision tree must appeared in data
Daniel@0: %       find threshhold in all cases that <= max_V
Daniel@0: %    change the split of T
Daniel@0: tree.nodes(tree.num_node).split_id=best_attr;
Daniel@0: tree.nodes(tree.num_node).split_threshhold=cur_best_threshhold; %for cts attribute only
Daniel@0: 
Daniel@0: %note: below threshhold rejust is linera search, so it is slow. A better method is described in paper "Efficient C4.5"
Daniel@0: %if (output_type==0)
Daniel@0: if (node_types(best_attr)==1)  %is a continuous attribute
Daniel@0:   %find the value that approximate best_threshhold from below (the largest that <= best_threshhold)
Daniel@0:   best_value=0;
Daniel@0:   for i=1:size(fam_ev,2)  %note: need to search in all cases for all tree, not just in cases for this node
Daniel@0:     val = fam_ev(best_attr,i);
Daniel@0:     if (val <= cur_best_threshhold & val > best_value) %val is more clear to best_threshhold
Daniel@0:       best_value=val;
Daniel@0:     end
Daniel@0:   end
Daniel@0:   tree.nodes(tree.num_node).split_threshhold=best_value; %for cts attribute only
Daniel@0: end
Daniel@0: %end
Daniel@0:   
Daniel@0: if (output_type == 0)
Daniel@0:   fprintf('Create node %d split at %d gain %8.4f Th %d. Class %d Cases %d Error %d \n',tree.num_node, best_attr, max_gain, tree.nodes(tree.num_node).split_threshhold, top1_class, size_t, size_t - top1_class_cases );
Daniel@0: else
Daniel@0:   fprintf('Create node %d split at %d gain %8.4f Th %d. Mean %8.4f Cases %d\n',tree.num_node, best_attr, max_gain, tree.nodes(tree.num_node).split_threshhold, avg_y_T, size_t );
Daniel@0: end
Daniel@0:   
Daniel@0: %(6) Foreach T' in the split_T
Daniel@0: %        if T' is Empty
Daniel@0: %            Child of node_id is a leaf
Daniel@0: %        else
Daniel@0: %            Child of node_id = split_tree (T')
Daniel@0: tree.nodes(new_node).is_leaf=0; %because this node will be split, it is not leaf now
Daniel@0: for i=1:size(best_split,2)
Daniel@0:   if (size(best_split{i},2)==0) %T(i) is empty
Daniel@0:     %create one new leaf node
Daniel@0:     tree.num_node=tree.num_node+1;
Daniel@0:     tree.nodes(tree.num_node).used=1; %flag this node is used (0 means node not used, it will be removed from tree at last to save memory)
Daniel@0:     tree.nodes(tree.num_node).is_leaf=1;
Daniel@0:     tree.nodes(tree.num_node).children=[];
Daniel@0:     tree.nodes(tree.num_node).split_id=0;
Daniel@0:     tree.nodes(tree.num_node).split_threshhold=0;  
Daniel@0:     if (output_type == 0)
Daniel@0:       tree.nodes(tree.num_node).probs=zeros(1,num_cat); %the prob for each value of class node 
Daniel@0:       tree.nodes(tree.num_node).probs(top1_class)=1; %use the majority class of parent node, like for binary class, 
Daniel@0:                                                    %and majority is class 2, then the CPT is [0 1]
Daniel@0:                                                    %we may need to use prior to do smoothing, to get [0.001 0.999]
Daniel@0:       tree.nodes(tree.num_node).error.self_error=0; 
Daniel@0:       tree.nodes(tree.num_node).error.all_error=0;  
Daniel@0:       tree.nodes(tree.num_node).error.all_error_num=0;
Daniel@0:     else
Daniel@0:       tree.nodes(tree.num_node).mean = avg_y_T; %just use parent node's mean value
Daniel@0:       tree.nodes(tree.num_node).std = std_T;
Daniel@0:     end
Daniel@0:     %add the new leaf node to parents
Daniel@0:     num_children=size(tree.nodes(new_node).children,2);
Daniel@0:     tree.nodes(new_node).children(num_children+1)=tree.num_node;
Daniel@0:     if (output_type==0)
Daniel@0:       fprintf('Create leaf node(nullset) %d. %d-th child of Father %d Class %d\n',tree.num_node, i, new_node, top1_class );
Daniel@0:     else
Daniel@0:       fprintf('Create leaf node(nullset) %d. %d-th child of Father %d \n',tree.num_node, i, new_node );
Daniel@0:     end
Daniel@0: 
Daniel@0:   else
Daniel@0:     if (node_types(best_attr)==0)  % if attr is discrete, it should be removed from the candidate set  
Daniel@0:       new_candidate_attrs = mysetdiff(candidate_attrs,[best_attr]);
Daniel@0:     else
Daniel@0:       new_candidate_attrs = candidate_attrs;
Daniel@0:     end
Daniel@0:     new_sub_node = split_dtree (CPD, fam_ev, node_sizes, node_types, stop_cases, min_gain, best_split{i}, new_candidate_attrs, num_cat);  
Daniel@0:     %tree.nodes(parent_id).error.all_error += tree.nodes(new_sub_node).error.all_error;
Daniel@0:     fprintf('Add subtree node %d to %d. #nodes %d\n',new_sub_node,new_node, tree.num_node );
Daniel@0: 
Daniel@0: %   tree.nodes(new_node).error.all_error_num = tree.nodes(new_node).error.all_error_num + tree.nodes(new_sub_node).error.all_error_num;
Daniel@0:     %add the new leaf node to parents
Daniel@0:     num_children=size(tree.nodes(new_node).children,2);
Daniel@0:     tree.nodes(new_node).children(num_children+1)=new_sub_node;
Daniel@0:   end
Daniel@0: end   
Daniel@0:   
Daniel@0: %(7) Compute errors of N; for doing pruning
Daniel@0: %    get the total error for the subtree
Daniel@0: if (output_type==0)
Daniel@0:   tree.nodes(new_node).error.all_error=tree.nodes(new_node).error.all_error_num/size_t;
Daniel@0: end
Daniel@0: %doing pruning, but doing here is not so efficient, because it is bottom up.
Daniel@0: %if tree.nodes()
Daniel@0: %after doing pruning, need to update the all_error to self_error
Daniel@0: 
Daniel@0: %(8) Return N
Daniel@0:   
Daniel@0: 
Daniel@0: 
Daniel@0: 
Daniel@0: %(1) For discrete output, we use GainRatio defined as below
Daniel@0: %  			         Gain(X,T)
Daniel@0: % 	GainRatio(X,T) = ----------
Daniel@0: % 			         SplitInfo(X,T)
Daniel@0: %   where
Daniel@0: %   Gain(X,T) = Info(T) - Info(X,T)
Daniel@0: %    				                       |Ti|
Daniel@0: % 	Info(X,T) = Sum for i from 1 to n of ( ---- * Info(Ti))
Daniel@0: %                                          |T|
Daniel@0:  			 
Daniel@0: %   SplitInfo(D,T) is the information due to the split of T on the basis
Daniel@0: %    of the value of the categorical attribute D. Thus SplitInfo(D,T) is
Daniel@0: %  		 I(|T1|/|T|, |T2|/|T|, .., |Tm|/|T|)
Daniel@0: %    where {T1, T2, .. Tm} is the partition of T induced by the value of D.
Daniel@0: 
Daniel@0: %   Definition of Info(Ti)
Daniel@0: %     If a set T of records is partitioned into disjoint exhaustive classes C1, C2, .., Ck on the basis of the 
Daniel@0: %     value of the categorical attribute, then the information needed to identify the class of an element of T 
Daniel@0: %     is Info(T) = I(P), where P is the probability distribution of the partition (C1, C2, .., Ck): 
Daniel@0: %     	P = (|C1|/|T|, |C2|/|T|, ..., |Ck|/|T|)
Daniel@0: %     Here I(P) is defined as
Daniel@0: %       I(P) = -(p1*log(p1) + p2*log(p2) + .. + pn*log(pn))
Daniel@0: % 
Daniel@0: %(2) For continuous output (regression tree), we use least squares score (adapted from Leo Breiman's book "Classification and regression trees", page 231
Daniel@0: %    The original support only binary split, we further extend it to permit multiple-child split
Daniel@0: %                                        
Daniel@0: %     Delta_R = R(T) - Sum for all childe nodes Ti (R(Ti))
Daniel@0: %     Where R(Ti)= 1/N * Sum for all cases i in node Ti ((yi - avg_y(Ti))^2)
Daniel@0: %     here N is the number of all training cases for construct the regression tree
Daniel@0: %          avg_y(Ti) is the average value for output variable for the cases in node Ti
Daniel@0: 
Daniel@0: function gain_score = compute_gain (fam_ev, node_sizes, node_types, T, info_T, attr_id, split_T, score_type, output_type)
Daniel@0: % COMPUTE_GAIN Compute the score for the split of cases T using attribute attr_id
Daniel@0: % gain_score = compute_gain (fam_ev, T, attr_id, node_size, method)
Daniel@0: %
Daniel@0: % fam_ev(i,j)  is the value of attribute i in j-th training cases, the last row is for the class label (self_ev)
Daniel@0: % T(i) is the index of i-th cases in current decision tree node, we need split it further
Daniel@0: % attr_id is the index of current node considered for a split
Daniel@0: % split_T{i} is the i_th subset in partition of cases T according to the value of attribute attr_id
Daniel@0: % score_type if 0, is gain ratio, 1 is information gain (only apply to discrete output)
Daniel@0: % node_size(i) the node size of i-th node in the family
Daniel@0: % output_type: 0 means discrete output, 1 means continuous output.
Daniel@0: gain_score=0;
Daniel@0: % ***********for DISCRETE output*******************************************************
Daniel@0: if (output_type == 0)
Daniel@0:   % compute Info(T)
Daniel@0:   total_cnt = size(T,2);
Daniel@0:   if (total_cnt==0)
Daniel@0:     return;
Daniel@0:   end;
Daniel@0:   %class_split_T = split_cases(fam_ev,node_sizes,node_types,T,size(fam_ev,1),0); %split cases according to class
Daniel@0:   %info_T = compute_info (fam_ev, T, class_split_T);
Daniel@0: 
Daniel@0:   % compute Info(X,T)
Daniel@0:   num_class = size(split_T,2); 
Daniel@0:   subset_sizes = zeros(1,num_class);
Daniel@0:   info_ti = zeros(1,num_class);
Daniel@0:   for i=1:num_class
Daniel@0:     subset_sizes(i)=size(split_T{i},2);
Daniel@0:     if (subset_sizes(i)~=0)
Daniel@0:       class_split_Ti = split_cases(fam_ev,node_sizes,node_types,split_T{i},size(fam_ev,1),0); %split cases according to class
Daniel@0:       info_ti(i) = compute_info(fam_ev, split_T{i}, class_split_Ti);
Daniel@0:     end
Daniel@0:   end    
Daniel@0:   ti_ratios = subset_sizes/total_cnt;  %get the |Ti|/|T|
Daniel@0:   info_X_T = sum(ti_ratios.*info_ti);
Daniel@0: 
Daniel@0:   %get Gain(X,T)
Daniel@0:   gain_X_T = info_T - info_X_T;
Daniel@0: 
Daniel@0:   if (score_type == 1) %information gain
Daniel@0:     gain_score=gain_X_T;
Daniel@0:     return;
Daniel@0:   end
Daniel@0:   %compute the SplitInfo(X,T)   //is this also for cts attr, only split into two subsets
Daniel@0:   splitinfo_T = compute_info (fam_ev, T, split_T);
Daniel@0:   if (splitinfo_T~=0)
Daniel@0:     gain_score = gain_X_T/splitinfo_T;
Daniel@0:   end
Daniel@0: 
Daniel@0: % ************for continuous output**************************************************
Daniel@0: else 
Daniel@0:   N = size(fam_ev,2);
Daniel@0: 
Daniel@0:   % compute R(Ti)
Daniel@0:   num_class = size(split_T,2); 
Daniel@0:   R_Ti = zeros(1,num_class);
Daniel@0:   for i=1:num_class
Daniel@0:     if (size(split_T{i},2)~=0)
Daniel@0:       cases_T = fam_ev(size(fam_ev,1),split_T{i});
Daniel@0:       avg_y_T = mean(cases_T);
Daniel@0:       sqr_T = cases_T - avg_y_T;
Daniel@0:       R_Ti(i) = sum(sqr_T.*sqr_T)/N;  % get R(Ti) = 1/N * SUM(y-avg_y)^2
Daniel@0:     end
Daniel@0:   end
Daniel@0:   %delta_R = R(T) - SUM(R(Ti))
Daniel@0:   gain_score = info_T - sum(R_Ti);
Daniel@0: 
Daniel@0: end
Daniel@0: 
Daniel@0: 
Daniel@0: %   Definition of Info(Ti)
Daniel@0: %     If a set T of records is partitioned into disjoint exhaustive classes C1, C2, .., Ck on the basis of the 
Daniel@0: %     value of the categorical attribute, then the information needed to identify the class of an element of T 
Daniel@0: %     is Info(T) = I(P), where P is the probability distribution of the partition (C1, C2, .., Ck): 
Daniel@0: %     	P = (|C1|/|T|, |C2|/|T|, ..., |Ck|/|T|)
Daniel@0: %     Here I(P) is defined as
Daniel@0: %       I(P) = -(p1*log(p1) + p2*log(p2) + .. + pn*log(pn))
Daniel@0: function info = compute_info (fam_ev, T, split_T)
Daniel@0: % COMPUTE_INFO compute the information for the split of T into split_T
Daniel@0: % info = compute_info (fam_ev, T, split_T)
Daniel@0: 
Daniel@0: total_cnt = size(T,2);
Daniel@0: num_class = size(split_T,2);
Daniel@0: subset_sizes = zeros(1,num_class);
Daniel@0: probs = zeros(1,num_class);
Daniel@0: log_probs = zeros(1,num_class);
Daniel@0: for i=1:num_class
Daniel@0:   subset_sizes(i)=size(split_T{i},2);
Daniel@0: end    
Daniel@0: 
Daniel@0: probs = subset_sizes/total_cnt;
Daniel@0: %log_probs = log2(probs);  % if probs(i)=0, the log2(probs(i)) will be Inf
Daniel@0: for i=1:size(probs,2)
Daniel@0:   if (probs(i)~=0)
Daniel@0:     log_probs(i)=log2(probs(i));
Daniel@0:   end
Daniel@0: end
Daniel@0: 
Daniel@0: info = sum(-(probs.*log_probs));
Daniel@0: