## set seed so the training/test sets can be reproduced
set.seed(1)



## read in data 
MC<-readRDS("speech-MFCCs.rds") 



## make some plots 
for(w in c("up","down")){ 

  png(paste0("speech-",w,".png"),height=4,width=8.5,units="in",res=200)  
  par(family="Times",mfrow=c(1,2),mar=c(3,3,1,.25),mgp=c(1.75,.75,0))

  ### mean plot
  Mw<-apply(MC[[w]],c(2,3),mean)  
  image(1:ncol(Mw),1:nrow(Mw),t(Mw[nrow(Mw):1,]),
        xlab="",ylab="",xaxt="n",yaxt="n") 
  axis(2,labels="time",at=50,line=1,tick=FALSE) 
  axis(1,labels="mfcc",at=6,line=1,tick=FALSE) 
  axis(2,labels=1:nrow(Mw),at=nrow(Mw):1,tick=FALSE,las=1,gap.axis=1) 
  axis(1,labels=1:ncol(Mw),at=1:ncol(Mw),tick=FALSE,las=1,gap.axis=0) 

  ### correlation plot
  Cw<-cor(t(apply(MC[[w]],1,"c")))
  image(t(Cw[nrow(Cw):1,]),xlab="mfcc",ylab="",xaxt="n",yaxt="n")    

  sep<-nrow(Mw)*(0:ncol(Mw))/nrow(Cw)
  abline(h=sep,col=gray(.25)) 
  abline(v=sep,col=gray(.25)) 
  axis(1,labels=1:ncol(Mw),at=sep[-1]-sep[2]/2,tick=FALSE,gap.axis=0)
  axis(2,labels=ncol(Mw):1,at=sep[-1]-sep[2]/2,tick=FALSE,las=2,gap.axis=0) 

  dev.off()
}



## construct training and test set, and vectorize
ntest<-100  
Ytest<-Ytrain<-list()

for(k in seq_along(MC)){ 
  itest<-sample(1:dim(MC[[k]])[1],ntest)  
  Ytest[[k]]<-t(apply(MC[[k]][itest,,],1,c))
  Ytrain[[k]]<-t(apply(MC[[k]][-itest,,],1,c))
}
names(Ytest)<-names(Ytrain)<-names(MC) 



##  mean and variance estimates 
M<-sapply(Ytrain,colMeans) 

V<-list() 
p1<-dim(MC[[1]])[2]
p2<-dim(MC[[1]])[3]
p<-p1*p2

for(k in seq_along(Ytrain)){  
  SMLE<-cov(Ytrain[[k]]) 
  SCSE<-covKCD::covCSE(SMLE,n=nrow(Ytrain[[k]]),p1=p1,p2=p2)
  SSEP<-covKCD::covKCD(SMLE,p1,p2)$K  
  V[[k]]<-list(SMLE=SMLE,SCSE=SCSE,SSEP=SSEP) 
  cat(k,"\n") 
}

colnames(M)<-names(V)<-names(MC)


## pooled variance estimate using Greene and Rayens (1989) estimator
source("speech-lpSAnu.r") 
SA<-array(sapply(V,function(v){ v$SMLE }),dim=c(p,p,length(V)))
nw<-sapply(Ytrain,nrow) 

fit<-optimize(lpSAnu,interval=c(p+3,6*p),SA=SA,n=nw,maximum=TRUE) 
nuHPP<-fit$max 
wHPP<-nw/(nw+nuHPP-p-1) 

### GR shrinkage point
S0<-apply(sweep(SA,3,wHPP/sum(wHPP),"*"),c(1,2),sum ) 

### GR shrinkage estimates 
for(k in seq_along(V)){ V[[k]]$SHPP<-(1-wHPP[k])*S0 + wHPP[k]*V[[k]]$SMLE } 



## classification of test set with QDA 
CA<-array(dim=c(length(V),length(V),length(V[[1]])))
for(j in 1:4){  
  for(k1 in seq_along(Ytest)){   
    Scores<-NULL
    for(k2 in seq_along(Ytest)){ 
      E<-sweep(Ytest[[k1]],2,M[,k2],"-") 
      SSE<-apply( t(E)*solve(V[[k2]][[j]],t(E)),2,sum) 
      dS<-determinant(V[[k2]][[j]])$mod 
      Scores<-cbind(Scores,c(SSE+dS)) 
    } 
  CA[k1,,j]<-(table( c(apply(Scores,1,which.min),1:length(V)))-1)/nrow(Scores)  
  cat(j,k1,"\n") 
  }
}
dimnames(CA)<-list(names(V),names(V),names(V[[1]])) 

## confusion array and summary
CA<-CA[,,c(3,2,1,4)] 
saveRDS(CA,file="speech-CA.rds") 
PCP<-apply(CA,3,diag)
apply(PCP,2,mean) 



## confusion array summaries 

### rate of correct classification
for(i in 1:nrow(PCP)){ 
cat( rownames(PCP)[i] ," & ",paste(PCP[i,],collapse=" & ")," \\\\ \n") 
}

##  plots of confusion matrices
mnames<-c("KMLE","CSE","MLE","PPE")
pdf("speech-confusion.pdf",height=8.5,width=9,family="Times")
par(mfrow=c(2,2),mar=c(3,1,1.35,1),mgp=c(1.75,.75,0),oma=c(.1,2,.1,.1) )
for(k in 1:4){ 
  image(t(CA[10:1,,k]),xaxt="n",yaxt="n")   
  if(is.element(k,c(3,4))) { 
  axis(side=1,at=seq(0,1,length=10),
       labels=dimnames(CA)[[1]],tick=FALSE,las=3 ,cex.axis=1.25)  }
  if(is.element(k,c(1,3))){ 
       axis(side=2,at=seq(1,0,length=10),
       labels=dimnames(CA)[[1]],tick=FALSE,las=1 ,cex.axis=1.25)  } 
  mtext(side=3,mnames[k],line=.2 )
}
dev.off()