


#Errors from ICA as Epifanio 2008:
eica=c(0.0560    ,0.0627   , 0.0613   , 0.0480  ,  0.0507,    0.0640 ,   0.0627,    0.0600    ,0.0613  ,  0.0773  ,  0.0533 ,
   0.0573  ,  0.0613   , 0.0547,    0.0933,
    0.0613,    0.0667  ,  0.0640  ,  0.0653 ,   0.0813   , 0.0520 ,   0.0720 ,   0.0573 ,   0.0640 ,   0.0680  ,  0.0507 ,   0.0840 ,
       0.0653  ,  0.0667  ,  0.0533,
    0.0640 ,   0.0627,    0.0653  ,  0.0840 ,   0.0880   , 0.0667 ,   0.0627   , 0.0467,    0.0600 ,   0.0773   , 0.0760 ,   0.0720 ,
       0.0627   , 0.0813  ,  0.0547,    0.0747 ,   0.0787  ,  0.0733 ,   0.0493   , 0.0613)

#ese: FAA-Basis
t.test(ese,eica,paired=T,alternative='less')

#        Paired t-test

#data:  ese and eica
#t = -0.88143, df = 49, p-value = 0.1912


#ese: FAA-Discretizing
 t.test(ese,eica,paired=T,alternative='less')

#        Paired t-test

#data:  ese and eica
#t = -0.92138, df = 49, p-value = 0.1807




#ese: FAA-Basis, esec: Constatini
t.test(ese,esec,paired=T,alternative='less')
#t = -2.4842, df = 49, p-value = 0.008225



##########################
#paired proportions

#basis

#points to FDA, coefficient basis

numbase=32#;%12,32,64,256
mibspline15=create.bspline.basis(c(1,21), numbase,6);
miLfdobj = int2Lfd(2);  
la = 10#;   %tb. 1
mifdPar = fdPar(mibspline15, miLfdobj, la);

norep=3 #20

ese=c()

tb=c()
K=3
start.time <- Sys.time()
for (i in 1:50){

we=read.table(casose[i])


 ex=we[,1:101]
 #nombre=['wt' num2str(i) '.txt'];
#wt=load(nombre);
#tex=wt(:,1:101);




 ti=seq(1,21,by=0.2)

direcsbspline15 = smooth.basis(ti,t(ex),mifdPar)




X=t(direcsbspline15$fd$coefs) #coefficeinets

PM=eval.penalty(mibspline15, Lfdobj=miLfdobj, rng=c(1,21))

set.seed(2015)
lass10 <- FstepArchetypesMod(data=X,k=K,verbose=FALSE,nrep=norep,PM=PM,saveHistory=F) 
ai <- bestModel(lass10[[1]]) #archetypes



a3=ai$archetypes #son coeficientes en base de archetypes; archetypes

##convertir en functions
#yaf=fd(t(a3),mibspline15)
#plot(yaf)


  data=X
huge=200
k=3 #3 archetypes

  n <- ncol(t(data))
  x_gvv <- rbind(t(data), rep(huge, n))
  zs <- rbind(t(a3), rep(huge, k))  
  zs <- as.matrix(zs)
  alphascce <- matrix(0, nrow = k, ncol = n)
  
  for (j in 1 : n){
   alphascce[, j] = coef(nnls(zs, x_gvv[,j]))
  }



#test
wt=read.table(casost[i])

tex=wt[,1:101]
direcsbspline15t = smooth.basis(ti,t(tex),mifdPar)




data=t(direcsbspline15t$fd$coefs) #coefficients

  
huge=200
k2=k

  n <- ncol(t(data))
  x_gvv <- rbind(t(data), rep(huge, n))
  
  alphascct <- matrix(0, nrow = k2, ncol = n)
  
  for (j in 1 : n){
   alphascct[, j] = coef(nnls(zs, x_gvv[,j]))
  }



tae=t(alphascce)
tat=t(alphascct)

#mirar donde esté minimo en tae para cada clase, para saber que clase en tat


#en entrenamiento
cle=max.col(-tae)

permuta=c()
#coinci=c()

for (ii in 1:3){
aux=0
 for ( jj in 1:3){

s=sum((cle==ii)== (we[,102]==jj))
if(aux<s){
aux=s
permuta[ii]=jj
}
}
}



cl=max.col(-tat)


clc=cl
clc[cl==1]=permuta[1]
 clc[cl==2]=permuta[2]
 clc[cl==3]=permuta[3]

tb=c(tb,clc==wt[,102])
ese[i]=sum(clc!=wt[,102])/750

cat("iteration = ", i,"mean ese",  mean(ese), "\n")
}
end.time <- Sys.time()
time.taken <- end.time - start.time


summary(ese)


#################
#discretizing, after smoothing


numbase=32#;%12,32,64,256
mibspline15=create.bspline.basis(c(1,21), numbase,6);
miLfdobj = int2Lfd(2);  
la = 10#;   %tb. 1
mifdPar = fdPar(mibspline15, miLfdobj, la);

norep=3 #20

ese=c()
td=c()

K=3
start.time <- Sys.time()
for (i in 1:50){

we=read.table(casose[i])


 ex=we[,1:101]
 #nombre=['wt' num2str(i) '.txt'];
#wt=load(nombre);
#tex=wt(:,1:101);




 ti=seq(1,21,by=0.2)

direcsbspline15 = smooth.basis(ti,t(ex),mifdPar)


fp=eval.fd(ti,direcsbspline15$fd) #evaluate the smoothed 


set.seed(2015)
lass10 <- stepArchetypesRawData(data = t(fp), numArch=K, numRep = norep,verbose = FALSE)
#lass10 <- FstepArchetypesMod(data=X,k=K,verbose=FALSE,nrep=norep,PM=PM,saveHistory=F) 
ai <- bestModel(lass10[[1]]) #archetypes



a3=ai$archetypes #son coeficientes en base de archetypes; archetypes



  data=t(fp)
huge=200
k=3 #3 archetypes

  n <- ncol(t(data))
  x_gvv <- rbind(t(data), rep(huge, n))
  zs <- rbind(t(a3), rep(huge, k))  
  zs <- as.matrix(zs)
  alphascce <- matrix(0, nrow = k, ncol = n)
  
  for (j in 1 : n){
   alphascce[, j] = coef(nnls(zs, x_gvv[,j]))
  }



#test
wt=read.table(casost[i])

tex=wt[,1:101]

direcsbspline15t = smooth.basis(ti,t(tex),mifdPar)


fpt=eval.fd(ti,direcsbspline15t$fd) #evaluate the smoothed 


data=t(fpt)

  
huge=200
k2=k

  n <- ncol(t(data))
  x_gvv <- rbind(t(data), rep(huge, n))
  
  alphascct <- matrix(0, nrow = k2, ncol = n)
  
  for (j in 1 : n){
   alphascct[, j] = coef(nnls(zs, x_gvv[,j]))
  }



tae=t(alphascce)
tat=t(alphascct)

#mirar donde esté minimo en tae para cada clase, para saber que clase en tat


#en entrenamiento
cle=max.col(-tae)

permuta=c()
#coinci=c()

for (ii in 1:3){
aux=0
 for ( jj in 1:3){

s=sum((cle==ii)== (we[,102]==jj))
if(aux<s){
aux=s
permuta[ii]=jj
}
}
}



cl=max.col(-tat)


clc=cl
clc[cl==1]=permuta[1]
 clc[cl==2]=permuta[2]
 clc[cl==3]=permuta[3]

td=c(td,clc==wt[,102])
ese[i]=sum(clc!=wt[,102])/750

cat("iteration = ", i,"mean ese",  mean(ese), "\n")
}
end.time <- Sys.time()
time.taken <- end.time - start.time

summary(ese)


###discretizing without smoothing



norep=3 #20

ese=c()
tr=c()

K=3
for (i in 1:50){

we=read.table(casose[i])


 ex=we[,1:101]



set.seed(2015)
lass10 <- stepArchetypesRawData(data = ex, numArch=K, numRep = norep,verbose = FALSE)
#lass10 <- FstepArchetypesMod(data=X,k=K,verbose=FALSE,nrep=norep,PM=PM,saveHistory=F) 
ai <- bestModel(lass10[[1]]) #archetypes



a3=ai$archetypes #son coeficientes en base de archetypes; archetypes



  data=ex
huge=200
k=3 #3 archetypes

  n <- ncol(t(data))
  x_gvv <- rbind(t(data), rep(huge, n))
  zs <- rbind(t(a3), rep(huge, k))  
  zs <- as.matrix(zs)
  alphascce <- matrix(0, nrow = k, ncol = n)
  
  for (j in 1 : n){
   alphascce[, j] = coef(nnls(zs, x_gvv[,j]))
  }



#test
wt=read.table(casost[i])

tex=wt[,1:101]

data=tex

  
huge=200
k2=k

  n <- ncol(t(data))
  x_gvv <- rbind(t(data), rep(huge, n))
  
  alphascct <- matrix(0, nrow = k2, ncol = n)
  
  for (j in 1 : n){
   alphascct[, j] = coef(nnls(zs, x_gvv[,j]))
  }



tae=t(alphascce)
tat=t(alphascct)

#mirar donde esté minimo en tae para cada clase, para saber que clase en tat


#en entrenamiento
cle=max.col(-tae)

permuta=c()
#coinci=c()

for (ii in 1:3){
aux=0
 for ( jj in 1:3){

s=sum((cle==ii)== (we[,102]==jj))
if(aux<s){
aux=s
permuta[ii]=jj
}
}
}



cl=max.col(-tat)


clc=cl
clc[cl==1]=permuta[1]
 clc[cl==2]=permuta[2]
 clc[cl==3]=permuta[3]

tr=c(tr,clc==wt[,102])
ese[i]=sum(clc!=wt[,102])/750

cat("iteration = ", i,"mean ese",  mean(ese), "\n")
}


summary(ese)


#################Constanini


la = 10#;   %tb. 1
numbase=32#;%12,32,64,256
mibspline15=create.bspline.basis(c(1,21), numbase,6);
miLfdobj = int2Lfd(2);  

mifdPar = fdPar(mibspline15, miLfdobj, la);

norep=3 #20

ese=c()
tc=c()

K=3
for (i in 1:50){

we=read.table(casose[i])


 ex=we[,1:101]
 #nombre=['wt' num2str(i) '.txt'];
#wt=load(nombre);
#tex=wt(:,1:101);




 ti=seq(1,21,by=0.2)

direcsbspline15 = smooth.basis(ti,t(ex),mifdPar)




X=t(direcsbspline15$fd$coefs) #coefficeinets



set.seed(2015)
lass10 <- stepArchetypes(X, k=K, verbose=FALSE,nrep=norep) 
ai <- bestModel(lass10[[1]]) #archetypes



a3=ai$archetypes #son coeficientes en base de archetypes; archetypes

##convertir en functions
#yaf=fd(t(a3),mibspline15)
#plot(yaf)


  data=X
huge=200
k=3 #3 archetypes

  n <- ncol(t(data))
  x_gvv <- rbind(t(data), rep(huge, n))
  zs <- rbind(t(a3), rep(huge, k))  
  zs <- as.matrix(zs)
  alphascce <- matrix(0, nrow = k, ncol = n)
  
  for (j in 1 : n){
   alphascce[, j] = coef(nnls(zs, x_gvv[,j]))
  }



#test
wt=read.table(casost[i])

tex=wt[,1:101]
direcsbspline15t = smooth.basis(ti,t(tex),mifdPar)




data=t(direcsbspline15t$fd$coefs) #coefficeinets

  
huge=200
k2=k

  n <- ncol(t(data))
  x_gvv <- rbind(t(data), rep(huge, n))
  
  alphascct <- matrix(0, nrow = k2, ncol = n)
  
  for (j in 1 : n){
   alphascct[, j] = coef(nnls(zs, x_gvv[,j]))
  }



tae=t(alphascce)
tat=t(alphascct)

#mirar donde esté minimo en tae para cada clase, para saber que clase en tat


#en entrenamiento
cle=max.col(-tae)

permuta=c()
#coinci=c()

for (ii in 1:3){
aux=0
 for ( jj in 1:3){

s=sum((cle==ii)== (we[,102]==jj))
if(aux<s){
aux=s
permuta[ii]=jj
}
}
}



cl=max.col(-tat)


clc=cl
clc[cl==1]=permuta[1]
 clc[cl==2]=permuta[2]
 clc[cl==3]=permuta[3]

tc=c(tc,clc==wt[,102])
ese[i]=sum(clc!=wt[,102])/750

cat("iteration = ", i,"mean ese",  mean(ese), "\n")
}


summary(ese)

 esec=ese


#Mcnemar tests

x=table(tb,td)

mcnemar.test(x)

#        McNemar's Chi-squared test with continuity correction

#data:  x
#McNemar's chi-squared = 0.0093168, df = 1, p-value = 0.9231


ti=read.table('ti.txt') #ti results ICA, fom matlab


 x2=table(tb,as.matrix(ti))
x3=table(td,as.matrix(ti))
x4=table(tr,as.matrix(ti))
x5=table(tb,tr)

library(exact2x2)

 exact2x2(x2, alternative = "less", paired=TRUE)

 #       Exact McNemar-type test

#data:  x2
#b = 762, c = 802, p-value = 0.162

 
exact2x2(x3, alternative = "less", paired=TRUE)

#        Exact McNemar-type test

#data:  x3
#b = 835, c = 879, p-value = 0.1495


exact2x2(x4, alternative = "greater", paired=TRUE)

#       Exact McNemar-type test

#data:  x4
#b = 972, c = 804, p-value = 3.662e-05




exact2x2(x5, alternative = "less", paired=TRUE)

#       Exact McNemar-type test

#data:  x5
#b = 765, c = 973, p-value = 3.331e-07



y=table(td,tc)
exact2x2(y, alternative = "less", paired=TRUE)

#        Exact McNemar-type test

#data:  y
#b = 521, c = 625, p-value = 0.001166


 y=table(tb,tc)
 exact2x2(y, alternative = "less", paired=TRUE)

#        Exact McNemar-type test

#data:  y
#b = 556, c = 656, p-value = 0.00222


