R/PAMfmado.R

In ExtremalDep: Extremal Dependence Models

PAMfmado <- function(x,K,J=0,threshold=.99,max.min=0){
  #
  # Description 
  #   This function performs the PAM algorithm based on the F-madogram distance
  #
  # Aguments 
  #   x a matrix
  #   K number of clusters
  #   J number of resampling for which the stations are randomly moved to break the dependence 
  #                    J=0 means no resampling
  #   threshold corresponding to the quantile level for the resampling.
  #   the resulting quantile is stored via sil.final<<-quantile(sil.vec[sil.vec>0],probs=threshold)
  #                    J=0 means no resampling
  # 
  Nnb = ncol(x) 
  Tnb = nrow(x)
  
  #--- DISTANCE MATRIX
  #--- F-MADOGRAM  
  V = array(NaN, dim = c(Tnb,Nnb))
  for(p in 1:Nnb) {
    x.vec = as.vector(x[,p])
    x.vec[x.vec < max.min]=NA # thresholding 
    Femp = ecdf(x.vec)(x.vec)
    V[,p] = Femp
  }
  # with Madogram
  DD = dist(t(V),method = "manhattan",diag = TRUE, upper = TRUE)/(2*Tnb)
  # !!! NA are excluded for when computing the distance between two rows but the dist is then multiplied by (Tnb/ nb of non-NA) 
  # in order to have comparable distances (ie wrt the number of data points) among all rows 
  
  #---------- CLUSTERING WITH PAM -------#
  output = pam(DD,K,diss = TRUE,medoids = NULL)
  
  #################################################################################
  ####  resampling step to compute the silhouette coeff for independent case
  #################################################################################
  if(J>0){
    sample.R = function(x) sample(x)
    sil.vec=rep(0,J)
    for(j in 1:J){
      x.sample=apply(x,2,sample.R) # resampling within each column 
      for(p in 1:Nnb) {
        x.bis = as.vector(x.sample[,p])
        x.bis[x.bis < max.min]=NA # thresholding 
        Femp = ecdf(x.bis)(x.bis)
        V[,p] = Femp
      }
      # with Madogram
      DD = dist(t(V),method = "manhattan",diag = TRUE, upper = TRUE)/(2*Tnb)
      sil.vec[j]=pam(DD,K,diss = TRUE,medoids = NULL,keep.diss=FALSE,keep.data=FALSE)$silinfo$avg.width
    }
    sil.final = quantile(sil.vec[sil.vec>0],probs=threshold)
    print("The quantile for resampling is ", sil.final, "\n", sep="")
  }
  
  return(output)
}