采用主成分法实现因子分析中的参数估计

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 }