###################################################
# Discriminating spectrometric curves: 
# R/S-plus commandlines
###################################################
#
# Entering spectrometric data 
#############################
SPECDAT <- as.matrix(read.table("npfda-spectrometric.dat"))
attributes(SPECDAT)$dimnames[[1]] <- character(0)
SPECURVES <- SPECDAT[,1:100]           # sample of curves
Dichotomous <- SPECDAT[,101]           # sample of scalar responses
Dichotomous[Dichotomous < 20] <- 1     # Class 1: fat content < 20%
Dichotomous[Dichotomous >= 20] <- 2    # Class 2: fat content >= 20%
##
# Computing predicted class number over 50 samples
##################################################
Misclas.lcv.of.spec.over.50.samples <- 0
subject1 <- (1:215)[Dichotomous==1]
subject2 <- (1:215)[-subject1]	 
for(i in 1:20){
        set.seed(sample(0:1000,1))
        learn1 <- sample(subject1,77)
        learn2 <- sample(subject2,43)
        learning <- c(learn1,learn2)
        testing <- (1:215)[-learning]
        SPECURVESL <- SPECURVES[learning,]
        SPECURVEST <- SPECURVES[testing,]
        Specdichl <- Dichotomous[learning]
        Specdicht <- Dichotomous[testing]
        res <- funopadi.knn.lcv(Specdichl,SPECURVESL,SPECURVEST,2,20,c(0,1),kind.of.kernel = "quadratic",semimetric="deriv")
        Misclas.lcv.of.spec.over.50.samples[i] <- sum(res$Predicted.classnumber != Specdicht)/length(Specdicht)
}
##
# Plotting the 50 misclassification rates
#########################################
boxplot(Misclas.lcv.of.spec.over.50.samples, ylab="Misclassification rates", cex=1)