採用主成分法實現因子分析中的參數估計

factpca<-function(x,score="Bartlett",rotation="varimax") {     if(!is.matrix(x)){     x<-as.matrix(x)       #x爲原始的數據矩陣     }     z<-scale(x,center=TRUE,scale=TRUE)   #將原始數據矩陣標準化     p<-ncol(x)    #求觀測變量的個數     if(p<3){         stop("factor analysis requires at least three variables")     }     cr<-cor(z)    #求相關係數矩陣     eig<-eigen(cr)  #求矩陣的特徵值與特徵向量     s=sum(eig$values)     tmp=0.0     flag=0     for(i in 1:length(eig$values)){        tmp=tmp+eig$values[i]        flag=i        if(tmp/s>=0.8)        break     }     rowname<-paste("X",1:p,sep="")     colname<-paste("Factor",1:flag,sep="")     A<-matrix(0,nrow=p,ncol=flag,dimnames=list(rowname,colname)) #構造因子載荷矩陣，初始化爲0     for(j in 1:flag){       A[,j]=sqrt(eig$values[j])*eig$vectors[,j]  #填充矩陣A的值     }     A2<-A                  #未旋轉的載荷矩陣     var.A<-diag(A%*%t(A))  #公共因子的方差     D1<- 1-var.A          #特殊因子的方差     rotmat<-matrix(0,ncol=flag,nrow=flag)     if(rotation=="varimax"){      #採用最大方差旋轉法，進行因子旋轉         #rot<-varimax(A)         rot<-do.call(rotation,c(list(A)))         A1<-rot$loadings           #A1爲旋轉後的載荷矩陣         rotmat<-rot$rotmat        #rotmat爲旋轉矩陣         A<-A1     }     D<-diag(D1)      #特殊因子方差矩陣     D2<-solve(D)     #特殊因子方差矩陣的逆矩陣     F<-0     if(score=="Bartlett"){         F<-solve(t(A)%*%D2%*%A)%*%t(A)%*%D2     }     else if(score=="regression"){         F<-t(A)%*%solve(cr)     }     scores<-F%*%t(z)     #list(rotloadings=A,loadings=A2,eigenvalues=eig$values,eigenvalue=eig$values(1:flag),comdgree=var.A,F=F)     result<-list(rotloadings=rot,loadings=A2,eigen=eig,comdgree=var.A,F=F,scores=t(scores))     result }