

#demo robust-ozone from archetypes
library('archetypes')



### Data set:

data('Ozone', package = 'mlbench')

oz <- Ozone[, -c(1, 2, 3, 9)]
oz <- na.omit(oz)
colnames(oz) <- c('OZONE', '500MH', 'WDSP', 'HMDTY', 'STMP',
                  'INVHT', 'PRGRT', 'INVTMP', 'VZBLTY')

oz0 <- scale(oz)



### The three original archetypes:
nrep=20
set.seed(1234)
a.oz <- stepArchetypes(oz0,k= 3,nrep=nrep)

parameters(bestModel(a.oz))
pb1=barplot(bestModel(a.oz), oz0, percentiles = TRUE)
pb1



### Data set with outliers:

set.seed(1234)
outliers <- t(sapply(runif(5, min = 1.5, max = 2),
                     function(x)
                     x * apply(oz, 2, max) + apply(oz, 2, IQR)))

oz1 <- scale(rbind(oz, outliers))


pairs(oz1)



### Original archetypal algorithm:

set.seed(1234)
a.oz1 <- stepArchetypes(oz1, k=3,nrep=nrep)

parameters(bestModel(a.oz1))
pb2=barplot(bestModel(a.oz1), oz1, percentiles = TRUE)

pb2
pb22=pb2[c(3,1,2),]


### Robust archetypal algorithm:

set.seed(1236)
#ra.oz1 <- robustArchetypes(oz1, 3)

ra.oz1=stepArchetypes(oz1, k=3,nrep=nrep,method=robustArchetypes)


parameters(bestModel(ra.oz1))
pb3=barplot(bestModel(ra.oz1), oz1, percentiles = TRUE)

pb3

###########################

#functions from Archetypalfunctions Directory

nrep=20
set.seed(1234)
ra.oz2=RstepArchetypesMod(oz1, k=3,nrep=nrep,PO=0.5,CT=6)
pb4=barplot(bestModel(ra.oz2), oz1, percentiles = TRUE)

pb4

#changing c in bisquare function
set.seed(1234)
ra.oz22=RstepArchetypesMod(oz1, k=3,nrep=nrep,PO=0.75,CT=1)
pb42=barplot(bestModel(ra.oz22), oz1, percentiles = TRUE)
pb42a=pb42
pb42[2,]=pb42a[1,]
pb42[1,]=pb42a[2,]



#changing c in bisquare function
set.seed(1234)
ra.oz23=RstepArchetypesMod(oz1, k=3,nrep=nrep,PO=0.5,CT=1)
pb43=barplot(bestModel(ra.oz23), oz1, percentiles = TRUE)
pb43a=pb43




#changing c in bisquare function
set.seed(1234)
ra.oz24=RstepArchetypesMod(oz1, k=3,nrep=nrep,PO=0.25,CT=1)
pb44=barplot(bestModel(ra.oz24), oz1, percentiles = TRUE)
pb44a=pb44



#changing c in bisquare function
set.seed(1234)
ra.oz25=RstepArchetypesMod(oz1, k=3,nrep=nrep,PO=0.25,CT=6)
pb45=barplot(bestModel(ra.oz25), oz1, percentiles = TRUE)
pb45a=pb45


#changing c in bisquare function
set.seed(1234)
ra.oz26=RstepArchetypesMod(oz1, k=3,nrep=nrep,PO=0.75,CT=6)
pb46=barplot(bestModel(ra.oz26), oz1, percentiles = TRUE)
pb46a=pb46
pb46[3,]=pb46a[1,]
pb46[1,]=pb46a[3,]


#########################
#Frobenius norm of differences

sqrt(sum(abs(pb1-pb22)^2))
sqrt(sum(abs(pb1-pb3)^2))
sqrt(sum(abs(pb1-pb4)^2))
sqrt(sum(abs(pb1-pb42)^2))
sqrt(sum(abs(pb1-pb43)^2))
sqrt(sum(abs(pb1-pb44)^2))
sqrt(sum(abs(pb1-pb45)^2))
sqrt(sum(abs(pb1-pb46)^2))




##############
#Composition

 bestModel(a.oz1)$beta[,331:335]
 bestModel(ra.oz1)$beta[,331:335]
 bestModel(ra.oz2)$beta[,331:335]
 bestModel(ra.oz22)$beta[,331:335]
