2016-10-27 7 views
0

私はオンラインでPCAをRでコーディングしようとしていますが、このコードは既存の実装が存在しないため、他の人にも便利です。擬似コードはhere(アルゴリズム1)にあります。次のように私はこれまでやったことです:オンラインPCA in R

PCA<-function(X,k,epsilon){ 
    X_f<-norm(as.matrix(X),"f") 
    d<-nrow(X) 
    n<-ncol(X) 
    l<-floor((8*k)/(epsilon^2)) 
    U<-matrix(0,d,l) 
    C<-matrix(0,d,d) 
    Y<-matrix(0,n,l) 
    for(t in 1:n){ 
     r<-X[,t]-(U%*%t(U)%*%X[,t]) 
     n<-C + r%*%t(r) 
     while(norm(n,"2") >= 2*(X_f^2)/l){ 
      lamb<-eigen(C)$values[1] 
      u<-eigen(C)$vectors[,1] 
      U<-cbind(U,u) 
      #U[,which(!apply(U==0,2,all))] 
      C<-C-(lamb*(u%*%t(u))) 
      r<-X[,t]-(U%*%t(U)%*%X[,t]) 
     } 
     C<-C+(r%*%t(r)) 
     y<-matrix(0,1,l)  
     y<-t(U)%*%x_t 
     Y[t,]<-y 
    } 
    return(Y) 
} 

私は有名な漁師のアイリスデータを使用したコードをテストするには:

log.ir <- log(iris[, 1:4]) 
ir.species <- iris[, 5] 

ir.pca <- PCA(log.ir,50,0.2) 

ないコードのバグ、があるようです私には明らかですが、whileループは決して止まらない、何か助けてもらえますか?

+1

アルゴリズムは一般的ではありませんが、同じ論文の付録でもっと便利な2つのアルゴリズムがあります – Jamil

答えて

1

それはあなたが、彼らは一般的に

function(X,k,epsilon){ 
    X_f<-norm(as.matrix(X),"f") 
    d<-nrow(X) 
    n<-ncol(X) 
    l<-floor((8*k)/(epsilon^2)) 
    U<-matrix(0,d,l) 
    C<-matrix(0,d,d) 
    Y<-matrix(0,n,l) 
    for(t in 1:n){ 
    r<-X[,t]-(U%*%t(U)%*%X[,t]) 
    n<-C + r%*%t(r) 
    while(norm(n,"2") >= 2*(X_f^2)/l){ 
     print(norm(n,"2")) 
     print(2*(X_f^2)/l) 
     lamb<-eigen(C)$values[1] 
     u<-eigen(C)$vectors[,1] 
     U<-cbind(U,u) 
     U[,which(!apply(U==0,2,all))] 
     C<-C-(lamb*(u%*%t(u))) 
     r<-X[,t]-(U%*%t(U)%*%X[,t]) 
    } 
    C<-C+(r%*%t(r)) 
    y<-matrix(0,1,l)  
    y<-t(U)%*%x_t 
    Y[t,]<-y 
    } 
    return(Y) 
} 

debug(PCA) 

を変更することはありませんことがわかりますdebug(PCA)while(norm(n,"2") >= 2*(X_f^2)/l)が終了したことがないので、あなたはこれらの値から印刷する場合2*(X_f^2)/l)は実際には常によりも小さいnorm(n,"2")

で、だ、とデバッグしたい関数の中のprintステートメントを使用すると、問題を診断する良い方法です。

+0

読者の方へ:この回答はアルゴリズムの正しい実装を示していません。 'ステートメント。 – knb