script_clustering_by_pam.r

   1 ##################################################################
   2 ## File: script_clustering_by_pam.r
   3 ##
   4 ## Description:   Using PAM to clustering conso using pam
   5 ## Last modified: jan 2012 by JC
   6 ##
   7 ##################################################################
   8
   9 ## Function: Converts a matrix of curves to the discret wavelet domain
  10 toDWT <- function(x, filter.number = 6, family = "DaubLeAsymm"){
  11   x2   <- spline(x, n = 2^ceiling( log(length(x), 2) ),
  12                  method = 'natural')$y
  13   Dx2 <- wd(x2, family = family, filter.number = filter.number)$D
  14   return(Dx2)
  15 }
  16
  17 ## Function: Computes the absolute contribution of the wavelet's scale
  18 ##           to the total total of the curve.
  19 contrib <- function(x) {
  20   J   <- log( length(x)+1, 2)
  21   nrj <- numeric(J)
  22   t0  <- 1
  23   t1  <- 0
  24   for( j in 1:J ) {
  25     t1     <- t1 + 2^(J-j)
  26     nrj[j] <- sqrt( sum( x[t0:t1]^2 ) )
  27     t0     <- t1 + 1
  28   }
  29   return(nrj)
  30 }
  31
  32 ## 1. Load libraries & data ####
  33 library(wavethresh)
  34 library(cluster)
  35
  36
  37 # powerload is a matrix that contains on each line one observation
  38 #           of length delta
  39
  40 ## 2. DWT  ##################
  41 delta <- ncol(powerload)
  42 n     <- nrow(powerload)
  43 Xdwt  <- t(apply(powerload, 1, toDWT)) # DWT over the lines of powerload.
  44 Xnrj  <- t(apply(Xdwt, 1, contrib))    # Absolute contribution to the energy.
  45
  46 ## 3. Cluster  ##############
  47 K <- 8 # Number of clusters
  48 Xnrj_dist <- dist(Xnrj)
  49 Xnrj_pam  <- pam(as.dist(Xnrj_dist), K)