R function with matrix argument to find eigenvalue and eigenvectors

Question

I am just writing a R function to do PCA on certain covariance matrix. Here when I tried to use eigen function to get eigenvalues and eigenvectors the compiler says the object could not be found. Is there any way I can fix the problem？ My code is listed below:

lab3<-function(cov,scale){ 
  if (scale==F)
  cov<-cov2cor(cov)

  dimension<-nrow（cov）
  eig<-eigen(cov)$values
  total<-sum(eig)
  stdev<-sqrt(eig)
  rotation<-eigen(cov)$vectors
  indp<-c(1:dimension)
  cump<-c(1:dimension)
  for （i in 1:dimension）
  {indp[i]=eig[i]/total
   cump[i]=sum(eigenvalue[1:i])/total
  }

  output=list(stdev,rotation,indp,cump)

}

Thanks.

The input can be just d<-matrix(c(1,-2,0,-2,5,0,0,0,2),3) then lab3(d,T) to run the code .. I am really sorry for the confusion caused. I didn't intend to do that and just I didn't know that. Thanks for your time.

Answer 1

The immediate problem is that in

cump[i]=sum(eigenvalue[1:i])/total

you refer to eigenvalue which does not exist. I presume you meant to use eig here instead:

cump[i]=sum(eig[1:i])/total

From the comment the error seems to be:

as.vector(x, mode) : cannot coerce type 'closure' to vector of type 'any'

This results I suspect because you call the function without specifying scale. R will then find the scale function (a closure) and that can't be coerced to a type needed for the if() statement. An easy way to solve this is to do one of:

lab3 <- function(cov, scale = FALSE) {
....

or

lab3 <- function(cov) {
  if(missing(scale))
    scale <- FALSE
....

with the first form preferred.

There are other issues;

1. surely you want

   if(scale)
     cov <- cov2cor(cov)
   

   ? I.e. only if you want all variables scaled to zero mean unit variance is correlation matrix required,

2. the for loop can be done more efficiently with these two lines:

   indp <- eig / total
   cump <- cumsum(indp)
   

   You there don't need the for loop at all, and you don't need to set up indp and cump first either.

3. you call eigen() twice. Better call it once and save the entire returned object. The subset that for the bits you want.

If I solve all these issues then we have the following function:

lab3 <- function(cov, scale=FALSE){ 
  if (scale)
    cov <- cov2cor(cov)

  ed <- eigen(cov)
  eig <- ed$values
  total <- sum(eig)
  stdev <- sqrt(eig)
  rotation <-ed$vectors
  indp <- eig / total
  cump <- cumsum(eig)

  list(stdev, rotation, indp, cump)
}

Which works:

> lab3(cov(iris[, 1:4]))
[[1]]
[1] 2.0562689 0.4926162 0.2796596 0.1543862

[[2]]
            [,1]        [,2]        [,3]       [,4]
[1,]  0.36138659 -0.65658877 -0.58202985  0.3154872
[2,] -0.08452251 -0.73016143  0.59791083 -0.3197231
[3,]  0.85667061  0.17337266  0.07623608 -0.4798390
[4,]  0.35828920  0.07548102  0.54583143  0.7536574

[[3]]
[1] 0.924618723 0.053066483 0.017102610 0.005212184

[[4]]
[1] 4.228242 4.470912 4.549122 4.572957

Finally, I'll note that doing PCA via a singular value decomposition is considered better than via an eigen decomposition for reasons of numerical stability. And you can do all of this via princomp or preferably prcomp, both in base R.