############################################################################################################# #### Classification Method using Confidence Set #### #### Part 1: read in the training data set cT #### #### cc: number of classes #### #### pp: feature vector dimension #### #### nni: number of training observations from the i-th class, i=1, ..., cc #### #### Part 2: computer the critical constant lambda, which depends on #### #### alpha: 1-alpha is the proportion of correct classifications of future objects #### #### gamma: confidence about the randomness in training data set cT #### #### cc, pp, nn1,...,nncc: specified in Part 1 above #### #### SS: number of outlayer simulations (affect the accuracy of lambda due to simul randomness) #### #### QQ: number of innerlayer simulations (affect the accuracy of lambda due to simul randomness) #### #### lambda=9.175 for seed 10,alpha=0.05,gamma=0.95,cc=3,pp=2,nni=50,SS=10000,QQ=10000 #### #### Part 3: Classify a future observaed object, with inputs #### #### yy : the vector of measurement(s) for new case #### #### cT data: the training dataset matrix which the 1st column is class #### ############################################################################################################# ################ #### Part 1 #### ################ getwd() setwd("\\\\filestore.soton.ac.uk/users/wl/mydocuments/mydocuments/MyResearch/My papers/inProgress/Classification/2018-8/R code") getwd() #To make sure that the dataset is in the right folder/directory iris<-read.csv("iris.csv",head=T) # Data used = Sepal length and width only irism<-iris[,c(1,2,5)] attach(irism) mean1<-apply(irism[Class=="Iris-setosa",1:2],2,mean) mean2<-apply(irism[Class=="Iris-versicolor",1:2],2,mean) mean3<-apply(irism[Class=="Iris-virginica",1:2],2,mean) sigma1<-cov(irism[Class=="Iris-setosa",1:2]) sigma2<-cov(irism[Class=="Iris-versicolor",1:2]) sigma3<-cov(irism[Class=="Iris-virginica",1:2]) ################ #### Part 2 #### ################ cc <- 3 pp <- 2 nn = c(50,50,50) alpha <- 0.05 gamma <- 0.95 SS <- 10000 # One may increase SS to improve the accuracy of lambda due to simulation randomness QQ <- 10000 # One may increase QQ to improve the accuracy of lambda due to simulation randomness library("MASS") # simulate mutivariate normal random vectors Ip <- diag(x = 1, nrow=pp, ncol=pp) # Start the clock! ptm <- proc.time() set.seed(10) #set.seed(5) #set.seed(18) #set.seed(100) #1. S loop lambda_s <- c() for (k in 1:SS ){ #2. l(cc) loop tl <- c() for (j in 1:cc){ u <- mvrnorm(n=1,rep(0,pp),Ip/nn[j]) v <- mvrnorm(n=(nn[j]-1),rep(0,pp),Ip) #vi=v[i,] #3. Q loop # calculate tlq for each q in Q loop w <- c() t <- c() sumvvt <- matrix(0,nrow=pp,ncol=pp) for (i in 1:(nn[j]-1) ){ sumvvt <- sumvvt + v[i,]%*%t(v[i,]) } for (i in 1:QQ){ w <- mvrnorm(n=1,rep(0,pp),Ip) t[i] <- t(w-u)%*%solve(1/(nn[j]-1)*sumvvt)%*%(w-u) } # ordered --> for finding tl[(1-a)Q] t <- sort(t) # find tl[(1-a)Q] for each l tl[j] <- t[round((1-alpha)*QQ)] } # max t -->lampda_s lambda_s[k] <- max(tl) } # ordered lampda_[s] lambda_s <- sort(lambda_s) lambda <- lambda_s[round((gamma)*SS)] lambda # Stop the clock time<-proc.time() - ptm time ################ #### Part 3 #### ################ ### Classify a given yy by checking whether it is in each of the cc acceptance sets ### yy <- c(4.5, 2.0) #Observed feature of a future object sigma1M <- solve(sigma1) sigma2M <- solve(sigma2) sigma3M <- solve(sigma3) A1<- sign( lambda - t(yy-mean1)%*%sigma1M%*%(yy-mean1) ) # +ve sign means "in the acceptance set" A2<- sign( lambda - t(yy-mean2)%*%sigma2M%*%(yy-mean2) ) A3<- sign( lambda - t(yy-mean3)%*%sigma3M%*%(yy-mean3) ) c(A1,A2,A3) ## When pp=2, the acceptance sets can be drawn as contours ## xliml = 3.8 xlimu = 8.5 yliml = 1.5 ylimu = 4.5 plot(yy[1],yy[2],type="p",pch=8,xlim=c(xliml,xlimu),ylim=c(yliml,ylimu),xlab="", ylab="") # Draw this future object yy #plot(yy[1],yy[2],type="p",pch=8,xlim=c(xliml,xlimu),ylim=c(yliml,ylimu)) # Draw this future object yy text(yy[1],yy[2]-0.2,"(4.5, 2.0)") ##Contour mesh x1 = seq(xliml, xlimu, 0.02) x2 = seq(yliml, ylimu, 0.02) I = length(x1) J = length(x2) Zz = matrix(0, nrow=I, ncol=J) ##Contour of A1 #Function fA1# fA1 <- function(yy){ val <- t(yy-mean1)%*%sigma1M%*%(yy-mean1) return(val) } for (i in 1:I){ for (j in 1:J){ Zz[i,j] = fA1(c(x1[i],x2[j])) } } v = c(lambda) fA1Cont <- contourLines(x1,x2,Zz,levels=v) lines(fA1Cont[[1]][[2]],fA1Cont[[1]][[3]],type="l") #contour(x1,x2,Zz,levels=v, xlab="x_1", ylab="x_2") text(mean1[1],mean1[2]+0.2,"Class 1") points(mean1[1],mean1[2],type="p", pch=3) #Centre of the ellips ##Contour of A2 #Function fA2# fA2 <- function(yy){ val <- t(yy-mean2)%*%sigma2M%*%(yy-mean2) return(val) } for (i in 1:I){ for (j in 1:J){ Zz[i,j] = fA2(c(x1[i],x2[j])) } } fA2Cont <- contourLines(x1,x2,Zz,levels=v) lines(fA2Cont[[1]][[2]],fA2Cont[[1]][[3]],type="l") text(mean2[1],mean2[2]-0.2,"Class 2") points(mean2[1],mean2[2],type="p", pch=3) #Centre of the ellips ##Contour of A3 #Function fA3# fA3 <- function(yy){ val <- t(yy-mean3)%*%sigma3M%*%(yy-mean3) return(val) } for (i in 1:I){ for (j in 1:J){ Zz[i,j] = fA3(c(x1[i],x2[j])) } } fA3Cont <- contourLines(x1,x2,Zz,levels=v) lines(fA3Cont[[1]][[2]],fA3Cont[[1]][[3]],type="l") text(mean3[1],mean3[2]+0.2,"Class 3") points(mean3[1],mean3[2],type="p", pch=3) #Centre of the ellips ######################################################################### ######################################################################### #### Now suppose all the 4 meansurements are used for classification #### #### We go through similar computation as above but without pictures #### ######################################################################### ######################################################################### # Data used = Sepal length and width detach(irism) attach(iris) mean1<-apply(iris[Class=="Iris-setosa",1:4],2,mean) mean2<-apply(iris[Class=="Iris-versicolor",1:4],2,mean) mean3<-apply(iris[Class=="Iris-virginica",1:4],2,mean) sigma1<-cov(iris[Class=="Iris-setosa",1:4]) sigma2<-cov(iris[Class=="Iris-versicolor",1:4]) sigma3<-cov(iris[Class=="Iris-virginica",1:4]) ################ #### Part 2 #### ################ cc <- 3 pp <- 4 nn = c(50,50,50) alpha <- 0.05 gamma <- 0.95 SS <- 10000 QQ <- 10000 library("MASS") # simulate mutivariate normal random vectors Ip <- diag(x = 1, nrow=pp, ncol=pp) # Start the clock! ptm <- proc.time() set.seed(100) #1. S loop lambda_s <- c() for (k in 1:SS ){ #2. l(cc) loop tl <- c() for (j in 1:cc){ u <- mvrnorm(n=1,rep(0,pp),Ip/nn[j]) v <- mvrnorm(n=(nn[j]-1),rep(0,pp),Ip) #vi=v[i,] #3. Q loop # calculate tlq for each q in Q loop w <- c() t <- c() sumvvt <- matrix(0,nrow=pp,ncol=pp) for (i in 1:(nn[j]-1) ){ sumvvt <- sumvvt + v[i,]%*%t(v[i,]) } for (i in 1:QQ){ w <- mvrnorm(n=1,rep(0,pp),Ip) t[i] <- t(w-u)%*%solve(1/(nn[j]-1)*sumvvt)%*%(w-u) } # ordered --> for finding tl[(1-a)Q] t <- sort(t) # find tl[(1-a)Q] for each l tl[j] <- t[round((1-alpha)*QQ)] } # max t -->lampda_s lambda_s[k] <- max(tl) } # ordered lampda_[s] lambda_s <- sort(lambda_s) lambda <- lambda_s[round((gamma)*SS)] lambda # Stop the clock time<-proc.time() - ptm time ################ #### Part 3 #### ################ ### Classify a given yy by checking whether it is in each of the cc acceptance sets ### yy <- c(4.5, 3.5, 1.4, 0.27) #Observed feature of a future object sigma1M <- solve(sigma1) sigma2M <- solve(sigma2) sigma3M <- solve(sigma3) A1<- sign( lambda - t(yy-mean1)%*%sigma1M%*%(yy-mean1) ) # +ve sign means "in the acceptance set" A2<- sign( lambda - t(yy-mean2)%*%sigma2M%*%(yy-mean2) ) A3<- sign( lambda - t(yy-mean3)%*%sigma3M%*%(yy-mean3) ) c(A1,A2,A3)