rm(list=ls())

# Initialize GINI first by running the function

GINI<-function(y,i){
  
  TabV<-table(y,i)
  nx<-apply(TabV,2,sum)
  pxy<-TabV/matrix(rep(nx,each=2),2,ncol(TabV))
  vxy<-pxy*(1-pxy)
  zx<-apply(vxy,2,sum)
  n<-sum(TabV)
  -sum(nx/n*zx)
  
}




  #B<-data.frame(Y,X1,X2)
  
  
  plot(X[,1],X[,2],col="white")
  points(X[(Y$Ca<0),1],X[(Y$Ca<0),2],col="blue",pch=".",cex=3)
  points(X[(Y$Ca>=0),1],X[(Y$Ca>=0),2],col="red",pch=".",cex=3)
  
  
  n<-dim(X)[2]
  #X<-X2
  u<-(sort(X[,2])[2:n]+sort(X[,2])[1:(n-1)])/2
  knot<-NULL
  for (s in 1:4){
    
    vgini<-rep(NA,length(u))
    
    for (i in 1:length(u)){
      kn<-c(knot,u[i])
      Fun<-function(x){sum(x<=kn)}
      I<-Vectorize(Fun)(X[,2])
      vgini[i]<-GINI(Y$Ca,I)
    }
    
    plot(u,vgini,type="l",col="green")
    
    k<-which.max(vgini)
    cat("knot",k,u[k],"\n")
    knot<-c(knot,u[k])
    u<-u[-k]
    
  }
  







## Compare to the "tree" package:
#library(tree)
#tree(Y~X2,method="gini")






n<-nrow(B)
u1<-(sort(X1)[2:n]+sort(X1)[1:(n-1)])/2
u2<-(sort(X2)[2:n]+sort(X2)[1:(n-1)])/2
gini<-matrix(NA,nrow(B)-1,2)

for(i in 1:length(u1)){
  I<-(X1<u1[i])
  gini[i,1]<-GINI(Y,I)
  I<-(X2<u2[i])
  gini[i,2]<-GINI(Y,I)
}

mg<-max(gini)
i<-1+sum(mg==max(gini[,2]))
par(mfrow = c(1, 2))

plot(u1,gini[,1],ylim=range(gini),col="green",type="b",xlab="X1",ylab="Gini index")
abline(h=mg,lty=2,col="red")

if(i==1){
  points(u1[which.max(gini[,1])],mg,pch=19,col="red")
  segments(u1[which.max(gini[,1])],mg,u1[which.max(gini[,1])],-100000)
}

plot(u2,gini[,2],ylim=range(gini),col="green",type="b",xlab="X2",ylab="Gini index")
abline(h=mg,lty=2,col="red")

if(i==2){
  points(u2[which.max(gini[,2])],mg,pch=19,col="red")
  segments(u2[which.max(gini[,2])],mg,u2[which.max(gini[,2])],-100000)
}

u2[which.max(gini[,2])]


## Compare to the "tree" package:
library(tree)
arbre<-tree(Y~X1+X2,data=B,method="gini")
arbre$frame[1:4,]
plot(arbre)
text(arbre)
partition.tree(arbre)