# Local Principal Curves # This R implementation refers to the papers # # [1] Local Principal Curves (Einbeck, Tutz & Evers, 2005), Statistics and Computing, # 15, 301--313. # and # [2] Exploring Multivariate Data Structures with Local Principal Curves # (Einbeck, Tutz & Evers, 2005b), Classification - The Ubiquitous Challenge. # Springer, Heidelberg. Proceedings of the GfKl 2004. # # The function lpc is implemented for branches up to second order and third depth # (for definitions of order and depth see [2]). # For depth =1 (i.e. for the situation in paper [1]), arbitrary dimensions d of the data cloud are allowed. # For depth >1, only bivariate data clouds are allowed (i.e. d=2). # # vs. 0.32 , 16/10/2007 # This software can be used without restrictions with the aim of research # or publication. However, please notify us if you *modify* it and work with # such modified versions: jochen.einbeck at durham.ac.uk # Authors: Jochen Einbeck, with contributions by Ludger Evers and Jo Dwyer. # In a second file (letters.r), we provide the code for the examples with # the simulated letters as in [2]. # Description of variables: # Xi: data matrix with N lines (observations) and d columns (variables). Note that d>2 restricts depth to 1. # h: bandwidth. May be either specified as a single number, then the same bandwidth is used in # all dimensions, or as a d-dimensional bandwidth vector. The default setting is 10% of the range in each direction. # t0: scalar step length. Default setting is t0=h, if h is a scalar, and t0=mean(h), if h is a vector. # iter: maximum number of steps which within each branch. Algorithm will stop automatically when the end of the branch is reached. # Hence, the final number of fitted LPC points per branch is usually much lower. # x0: optionally, one can set one or more starting points here. If the number of starting points specified here is smaller than mult, # then the missing points will be set at random. For example, if d=2, mult=3, and x0=c(58.5, 17.8, 80,20), then # one gets the starting points (58.5, 17.8), (80,20), and a randomly chosen third one. x0 can also be provided as a matrix, # with one row for each starting point. # way: "one": go only in direction of the first local eigenvector, # "back": go only in opposite direction, # "two": go from starting point in both directions. # pen: 'TRUE' activates angle penalization, 'FALSE' deactivates. # penpot: power used for angle penalization (see [1]). # mult: number of starting points. They are chosen at random if no starting point x0[j,] is specified, otherwise they correspond to #the first mult=n rows of the matrix x0. # depth: maximum depth of branches (phimax in [2]), restricted to the values 1,2 or 3, where 2 and 3 are only allowed for d=2. # thresh: threshold for the pruning mechanism (usually, between 0.0 and 0.3). 0.0 means no pruning. # rho0: threshold for the decision to launch a starting point of higher depth (usually, between 0.3 and 0.4). # plotlpc: if 0, no graphical output is given. # if 1, the LPC is given in form of points; # if 2, in form of lines; # if 3, in form of lines and points (only for d=2). # scaled: 'TRUE' activates the scaling of each column of Xi by its range, 'FALSE' (default) deactivates. # ...: All regular graphical arguments apply, for example specifying the limits of each axis with xlim, ylim, or their labels with #xlab, ylab. # Output: The function yields as output for d=2 a list containing # 'LPC': * the LPC itself, including all branches # For each fitted values, this matrix gives # (1) the unit-speed parameter lambda centered at the first fitted LPC point of each branch # (2) a unique branch number (which the software assigns autimatically to the branch) # (3) the depth of the corresponding branch # (4) the order of the corresponding branch # (5) the initialization number. For example, init=4 means that this point stems from the 4th random starting point initialization. # Note that (1) and (2) together form a unique parameterization of the data, which can be used as a dimension-reduced version of the original data # 'rho': * the values of rho, an indicator for the location of junctions on the initial curve (see [2]) # 'starting points depth 1', # 'starting points depth 2', # 'starting points depth 3': starting points of corresponding depth; # and for d >2 just a matrix containing the LPC and the starting points (of depth 1). # Moreover, a graphical output is given. # In the case d=2, one obtains a plot showing the computed local principal curve (+) and the # starting points of depth 1,2,3 , whereby (using the standard R graphics device) # the colors red, green and blue correspond to depths 1,2 and 3, respectively. lpc <- function(Xi,h, t0=mean(h),iter=100, x0, way = "two", pen=TRUE, penpot=2, mult=1, depth=1, thresh=0.0, rho0=0.4, plotlpc=3, scaled=FALSE, ...){ Xi <- as.matrix(Xi) N <- dim(Xi)[1] d <- dim(Xi)[2] X1 <- matrix(0,mult,d) # collects starting points for branches of depth 1. X2<-X3 <- matrix(0,0,d) # collect starting points for branches of depth 2 and 3. S2<-S3 <- matrix(0,0,d) # collects candidates for starting points of depth 2 and 3. saveall <- matrix(0,0,d) # stores the LPC countb <- 0 # counter for branches 17/07/2007 Lambda <- matrix(0,0,5) # stores the parameters (lambda, countb) as well as depth, order, and initialization number of the corresponding branch. s1 <- apply(Xi, 2, function(dat){ diff(range(dat))}) # range if (missing(h)){ if (!scaled){h <- s1/10} else {h<- 0.1}} # bandwidth by default: 10% of range if (length(h)==1){h <- rep(h,d)} if (scaled){ #scales the data to lie by its range Xi <- sweep(Xi, 2, s1, "/") } if (d==1){ stop("Data set needs to consist of at least two variables!") } else if (d > 2 && depth > 1){ cat("Dimension of data set enforces depth=1! \n") depth <- 1 } else if (d == 2 && plotlpc > 0) { plot(Xi[,1],Xi[,2], pch=21, col="grey55", ...) } if (!missing(x0)){ x0 <- matrix(x0, ncol=d, byrow=TRUE) if (dim(x0)[1] < mult) {n <- runif(mult-dim(x0)[1],1,N+1)%/%1; x0 <- rbind(x0,Xi[n,])} # creates the appropriate number of random starting points } else { n <- runif(mult,1,N+1)%/%1; x0 <- matrix(Xi[n,], length(n),d)# corrected 16/10 } #print(x0) for (j in 1:mult){ xo <- x0[j,] # defines appropraite starting point for the jth curve X1[j,] <- xo # adds starting point to list curve0 <- followx(Xi, xo, h, t0, iter, way, pen, penpot, phi =1, 0,rho0, plotlpc) # computes LPC of depth 1. saveall <- rbind(saveall,curve0[[1]]) # stores LPC l <- dim(curve0[[5]])[1] # number of candidates for junctions Lambda <- rbind(Lambda, cbind(curve0[[6]], countb,1,1,j)) countb <- countb + 1 dimnames(Lambda)[[2]]<- c("lambda", "branch", "depth", "order", "init") if (depth >1 && l>=1 ){ # constructs branches of depth 2. for (s in 1:l){ x <- curve0[[5]][s,] # candidates for junctions on initial curve. center.x <- cov.wt(Xi, wt= kernd(Xi,x,h)) # computes local covariance and mean at x eigen.cov <- eigen(center.x[[1]]) # eigenvalues and eigenvec's of local cov matrix eigen.vecd <- eigen.cov[[2]][,2] # second local eigenvector new.x <- center.x[[2]]+ 2*t0 * eigen.vecd # double stepsize to escape from initial curve if (kde2x(Xi,new.x, h) >= thresh){ # Pruning X2 <- rbind(X2, t(new.x)) # adds new starting point to list x <- new.x curve1 <- followx(Xi, x, h, t0, iter, way ="one", pen, penpot, phi=2, lasteigenvector= eigen.vecd, rho0, plotlpc) saveall<- rbind(saveall, curve1[[1]]) # stores LPC S2 <- rbind(S2,curve1[[5]]) # stores candidates for further junctions Lambda <- rbind(Lambda, cbind(curve1[[6]], countb, 2, 2,j)) countb <- countb + 1 if (d==2){ if (plotlpc > 1){ lines(curve1[[1]][,1], curve1[[1]][,2]) } } } new.x <- center.x[[2]]- 2*t0 * eigen.vecd # go in opposite direction if (kde2x(Xi,new.x, h) >= thresh){ # Pruning X2 <- rbind(X2, t(new.x)) x <- new.x curve1 <- followx(Xi, x, h, t0, iter, way="back", pen, penpot, phi=2, lasteigenvector=-eigen.vecd, rho0, plotlpc) saveall<- rbind(saveall, curve1[[1]]) S2 <- rbind(S2, curve1[[5]]) Lambda <- rbind(Lambda, cbind(curve1[[6]], countb, 2,2,j)) countb <- countb + 1 if (d==2){ if (plotlpc > 1){ lines(curve1[[1]][,1], curve1[[1]][,2]) } } # end if d==2 } # end if kde2x..... } # end for (s) }# end if (depth) k <- dim(S2)[1] # no. of candidates for starting points of depth 3. if (depth >2 && k>=1 ){ # constructs branches of depth 3 for (s in 1:k){ x <- S2[s,] # candidates for junctions on curve of depth 2. center.x <- cov.wt(Xi, wt= kernd(Xi,x,h)) eigen.cov <- eigen(center.x[[1]]) eigen.vecd<- eigen.cov[[2]][,2] new.x <- center.x[[2]]+ 2*t0 * eigen.vecd if (kde2x(Xi,new.x, h) >= thresh){ X3 <- rbind(X3, t(new.x)) x <- new.x curve1 <- followx(Xi, x , h, t0, iter, way ="one", pen, penpot, phi=3, lasteigenvector= eigen.vecd, rho0, plotlpc) saveall <- rbind(saveall, curve1[[1]]) S3 <- rbind(S3,curve1[[5]]) Lambda <- rbind(Lambda, cbind(curve1[[6]], countb, 3,2,j)) countb <- countb + 1 if (d==2){ if (plotlpc > 1){ lines(curve1[[1]][,1], curve1[[1]][,2]) } } } new.x <- center.x[[2]]- 2*t0* eigen.vecd if (kde2x(Xi,new.x, h) >= thresh){ X3 <- rbind(X3, t(new.x)) x <- new.x curve1 <- followx(Xi, x, h, t0, iter, way="back", pen, penpot, phi=3, lasteigenvector=-eigen.vecd, rho0, plotlpc) saveall<- rbind(saveall, curve1[[1]]) S3 <- rbind(S3,curve1[[5]]) Lambda <- rbind(Lambda, cbind(curve1[[6]], countb, 3,2,j)) countb <- countb + 1 if (d==2){ if (plotlpc > 1){ lines(curve1[[1]][,1], curve1[[1]][,2]) } } # end if d==2 } # end if kde2x.... } # end for (s) } # end if (depth) if (d==2){ if (plotlpc > 1){ lines(curve0[[1]][,1], curve0[[1]][,2]) } } # end if d==2 } # end for (j) if (d==2){ fit<- list("LPC"=saveall, "Parameterization"= Lambda, "h"=h, "t0"=t0, "dim"=dim(Xi), "rho"=curve0[[4]], "starting points depth 1:" = X1, "starting points depth 2:"= X2, "starting points depth 3:"= X3) } else { if (plotlpc > 0){ pairs(saveall, panel= if (plotlpc==1) "points" else "lines", labels= dimnames(Xi)[[2]]) } fit<- list("LPC"=saveall, "Parameterization"= Lambda, "h"=h, "t0"=t0, "dim"=dim(Xi), "starting points depth 1:" = X1) } class(fit)<-"lpc" fit } # end function # The following function calculates a branch of a LPC from a starting point x0 until it reaches the end of the cloud. # The function lpc.mod only uses the options way = 'one' and way = 'back', # though followx might run as a stand-alone function with the setting way= 'two'. # phi corresponds to the depth of the current branch. followx<-function(Xi, x0, h, t0, iter, way, pen=TRUE, penpot=2, phi =1, lasteigenvector=0, rho0=0.4, plotlpc=1){ #Xi <- as.matrix(Xi) N <- dim(Xi)[1] d <- dim(Xi)[2] if (way=="one"){ b<-1} else{ b<-2} highrhopoints <- matrix(0,0,d) # here points with high local 2nd eigenvalue will be stored save.xd <- eigen.vecd<-matrix(0, b*iter,d) # here the LPC will be stored cos.alt.neu <- rep(1, b*iter) lambda <- rep(NA, b*iter) # Parameterization rho <- rep(0, b*iter) # construction of empty vectors if (plotlpc %in% c(1,3)){ if (phi==1){start.pch="1"} else if (phi ==2){start.pch="2"} else {start.pch="3"} if (d==2){points(x0[1],x0[2],pch=start.pch,col=phi +1)} # plots the starting points of corresponding depth. } x <- x0 if (way == "one" || way =="two"){ for (i in iter:1){ center.x <- cov.wt(Xi, wt= kernd(Xi,x,h)) # local covariance and mean at x mu.x <- center.x[[2]] # Local center of mass at x (mean shift). if (d==2 && plotlpc %in% c(1,3)){ points(mu.x[1],mu.x[2], pch="+", col=phi+1) # plots LPC by sequence of +'s } save.xd[i,] <- mu.x # stores the first point of that branch. if (i==iter){ lambda[i] <- 0 } else { lambda[i] <- lambda[i+1] + sqrt(sum( (mu.x-save.xd[i+1,])^2)) }# measures distance for parameterization eigen.cov <- eigen(center.x[[1]]) # local covariance at x eigen.vecd[i,] <- eigen.cov[[2]][,1] # local first eigenvector at x rho[i]<-eigen.cov[[1]][2]/eigen.cov[[1]][1] # computes rho if ((i< iter) && (rho[i] > rho0) && (rho[i+1]<= rho0)){highrhopoints<-rbind(highrhopoints, x) } # compares with rho0 if (i == iter && lasteigenvector[1] !=0){ cos.alt.neu[i]<- sum(lasteigenvector* eigen.vecd[i,]) } if (i < iter){ cos.alt.neu[i]<- sum(eigen.vecd[i+1,]* eigen.vecd[i,])} # angle between current and previous eigenvector if (cos.alt.neu[i]<0){ eigen.vecd[i,]<- - eigen.vecd[i,]} # signum flipping if (pen==T && i iter +1 && rho[i] > rho0 && rho[i-1]<= rho0){highrhopoints<-rbind(highrhopoints,x)} if (i == (1 + iter) && lasteigenvector [1] != 0){ cos.alt.neu[i]<- -sum(lasteigenvector* eigen.vecd[i,])} if (i >= (2 + iter)){ cos.alt.neu[i]<- sum(eigen.vecd[i-1,]* eigen.vecd[i,])} if (cos.alt.neu[i]<0){ eigen.vecd[i,]<- - eigen.vecd[i,]} if (pen==T && i>=(2+iter)){ a<-(abs(cos.alt.neu[i]))^penpot eigen.vecd[i,]<- a* eigen.vecd[i,]+(1-a)* eigen.vecd[i-1,] } if (pen==T && i == 2 +iter && lasteigenvector[1] !=0 ){ a<-(abs(cos.alt.neu[i]))^penpot eigen.vecd[i,]<-a * eigen.vecd[i,] + (1-a)*lasteigenvector } new.x <- center.x[[2]] - eigen.vecd[i,]*t0 #lambda[i]<- -sqrt(sum( (x0-new.x)^2)) x <- new.x #if (is.na(lambda[i-1])){lambda[i-1]<-0} if ( i!= iter && !is.na(lambda[i-1]) && abs(min(lambda[i],lambda[i-1]))!=0 && abs(lambda[i]-lambda[i-1])/abs(lambda[i]+lambda[i-1]) <0.00001){break} } # for (i) # if (way=="back"){save.xd<-save.xd[(iter+1):(2*iter),] ; # eigen.vecd<- eigen.vecd[(iter+1):(2*iter),] # } # removed 17/07/2007 }# if (way) #print(lambda) filter <- !is.na(lambda); if (way =="two" ){filter[iter+1]<- FALSE}### list(save.xd[filter,],eigen.vecd[filter,],cos.alt.neu[filter], rho[filter], highrhopoints, round(-lambda[filter]+max(lambda[filter]), digits=4)) }# function #------------------------------------------------------------------------------------------ # Gaussian kernel function kern<- function(xi, x = 0, h = 1){ 1/h * dnorm((x - xi)/h) } # Multivariate kernel function kernd <- function(Xi,x,h){ d<-length(x) k<-1 for (j in 1:d){k<- k* kern(Xi[,j],x[j],h[j])} k } # kernel density estimator kde2x<-function(Xi, x, h){ l<-dim(Xi)[1] 1/l*sum(kern(Xi[,1],x[1],h[1])*kern(Xi[,2],x[2],h[2])) }