本次實驗使用的資料也是主成分分析實體的學生成績資料. pca
data=read.csv("score.csv",header=T)
cor(data)
source("msaR.R") # 呼叫自定義函式 （放在最后了）
fac0=msa.fa(data,2,rotation="none") # 主成分法，且不做因子旋轉
fac0
# 簡單驗證一下結果里都是什么
a = eigen(cor(scale(data))) # 相關系數矩陣的特征值、特征向量
a$values[1:2] # 前兩個特征值
sum((sqrt(3.71)*a$vectors[,1])**2) # 載荷矩陣列平方和
(sqrt(3.71)*0.4121)**2+(sqrt(1.262)*(-0.376))**2 # 行平方和=共同度
sqrt(3.71)*a$vectors[,1] # 載荷矩陣第一列 sqrt(特征值)*特征向量
factanal(data,factors=2,rotation="none") # 極大似然法 不做因子旋轉
fac1=msa.fa(data,2,rotation="varimax") # 用主成分法采用方差最大化作因子正交旋轉
fac1
plot(fac1$loadings,xlab="Factor1",ylab="Factor2") # 因子載荷圖
rnames=c("數學", "物理", "化學", "語文", "歷史", "英語")
text(fac1$loadings[,1],fac1$loadings[,2],labels = rnames, adj=c(0.5, -0.5))
biplot(fac2$scores,fac2$loadings) #畫出各個學生的因子得分圖和原坐標在因子的方向，全面反映了因子與原始資料的關系
樣本相關系數矩陣：
根據樣本的相關系數矩陣可以粗略的認為保留兩個因子.
主成分法估計載荷矩陣：
極大似然法估計載荷矩陣：
從上述極大似然法和主成分法得出的結果可以看出，極大似然法前兩個因子的累積貢獻率為74.5%，而主成分分析法的累計貢獻率則達到了82.87%，說明主成分法效果比極大似然法好. 原因在于，實際資料大多數很難服從多元正態分布的要求，而極大似然法做因子分析要求資料樣本服從多元正態分布.
為了使結果的結構更簡單，更容易解釋，所以對因子進行正交旋轉，結果如下：
語文(
x
4
x_4
x4?)、歷史(
x
5
x_5
x5?)、英語(
x
6
x_6
x6?)在第一個因子的載荷分別為0.8763，0.9174，0.9253，這三個指標都反應了學生的文科水平，所以我們將其命名為“文科因子”；數學(
x
1
x_1
x1?)、物理(
x
2
x_2
x2?)、化學(
x
3
x_3
x3?)、在第二個因子的載荷分別為0.8390，0.7837，0.8967，這三個指標則反應了學生的理科水平，所以我們將其命名為“理科因子”. 文科因子得分越高說明該學生文科成績越好，理科因子得分越高說明該學生理科成績越好.
從因子載荷圖也可以看出，語文、英語、歷史離第一個因子所代表的坐標軸近，數學、物理、化學則離第二個因子所代表的坐標軸近，因子旋轉的作用就體現在這里，讓載荷矩陣的列內方差最大，使得原始變數在因子空間可以離坐標軸更近. 這里做的是正交旋轉，也就是坐標軸在旋轉的程序中始終保持垂直，所以共同度是不會改變的，
各個學生的因子得分圖中可以看出33號同學的文科、理科因子得分都挺高，說明33號的文科理科成績都好（均衡表現好，可以回傳PCA看一下第二主成分的結果. pca）；3號位于右下角，文科、理科因子得分都偏低，說明文、理成績都不好；其他位于左上角和右下角的同學存在偏科的情況.
函式來自《多元統計分析及R語言建模》（王斌會編著，暨南大學出版社，2011年）
# msaR.R
options(digits=4)
#par(mar=c(4,4,2,1),cex=0.8) #鐠佸墽鐤嗛崶鎯ц埌鏉堝綊妾崪灞界摟嫻ｆ挸銇囩亸?
msa.X<-function(df){
X=df[,-1];
rownames(X)=df[,1];
X
}
msa.andrews<-function(x){
# x is a matrix or data frame of data
if (is.data.frame(x)==TRUE)
x<-as.matrix(x)
t<-seq(-pi, pi, pi/30)
m<-nrow(x); n<-ncol(x)
f<-array(0, c(m,length(t)))
for(i in 1:m){
f[i,]<-x[i,1]/sqrt(2)
for( j in 2:n){
if (j%%2==0)
f[i,]<-f[i,]+x[i,j]*sin(j/2*t)
else
f[i,]<-f[i,]+x[i,j]*cos(j%/%2*t)
}
}
#plot(c(-pi,pi), c(min(f),max(f)), type="n", xlab="t", ylab="f(t)")
plot(c(-pi,pi), c(min(f),max(f)), type="n", xlab="", ylab="")
for(i in 1:m) lines(t, f[i,] , col=i)
legend(2,max(f),rownames(x),col=1:nrow(x),lty=1:nrow(x),bty='n',cex=0.8)
}
msa.coef.sd<-function(fm){
b=fm$coef;b
si=apply(fm$model,2,sd);si
bs=b[-1]*si[-1]/si[1]
bs
}
msa.cor.test<-function(X,diag=TRUE){
p=ncol(X);
if(diag){
tp=matrix(1,p,p);
for(i in 1:p){
for(j in 1:i) tp[i,j]=cor.test(X[,i],X[,j])$stat;
for(j in i:p) tp[i,j]=cor.test(X[,i],X[,j])$p.value;
}
cat("corr test: \n");
tp=round(matrix(tp,p,dimnames=list(names(X),names(X))),4)
print(tp)
#return(tp)
cat("lower is t value, upper is p value \n")
} else {
cat("\n corr test: t value, p value \n");
if(is.matrix(X)) var=1:p
else var=names(X);
for(i in 1:(p-1)){
for(j in (i+1):p)
cat(' ',var[i],'-',var[j],cor.test(X[,i],X[,j])$stat,cor.test(X[,i],X[,j])$p.value,"\n")
}
}
}
msa.pca<-function(X,cor=FALSE,m=2,scores=TRUE,ranks=TRUE,sign=TRUE,plot=TRUE){
if(m<1) return
PC=princomp(X,cor=cor)
Vi=PC$sdev^2
Vari=data.frame('Variance'=Vi[1:m],'Proportion'=(Vi/sum(Vi))[1:m],
'Cumulative'=(cumsum(Vi)/sum(Vi))[1:m])
cat("\n")
Loadi=as.matrix(PC$loadings[,1:m])
Compi=as.matrix(PC$scores[,1:m])
if(sign)
for (i in 1:m)
if(sum(Loadi[,i])<0){
Loadi[,i] = -Loadi[,i]
Compi[,i] = -Compi[,i]
}
pca<-NULL
pca$vars=Vari
if(m<=1) pca$loadings = data.frame(Comp1=Loadi)
else pca$loadings = Loadi;
if(scores & !ranks) pca$scores=round(Compi,4)
if(scores & plot){
plot(Compi);abline(h=0,v=0,lty=3)
text(Compi,row.names(X))
# par(mar=c(4,4,2,3))
# biplot(Compi,Loadi); abline(h=0,v=0,lty=3)
# par(mar=c(4,4,1,1))
}
if(scores & ranks){
pca$scores=round(Compi,4)
Wi=Vi[1:m];Wi
Comp=Compi%*%Wi/sum(Wi)
Rank=rank(-Comp)
pca$ranks=data.frame(Comp=round(Comp,4),Rank=Rank)
}
pca
}
msa.fa<-function(X,m=2,scores=TRUE,rotation="varimax",common=TRUE,ranks=TRUE){
if(m<1) return
cat("\n")
S=cor(X);
p<-nrow(S); diag_S<-diag(S); sum_rank<-sum(diag_S)
rowname = names(X)
colname<-paste("Factor", 1:p, sep="")
A<-matrix(0, nrow=p, ncol=p, dimnames=list(rowname, colname))
eig<-eigen(S)
for (i in 1:p)
A[,i]<-sqrt(eig$values[i])*eig$vectors[,i]
for (i in 1:p) { if(sum(A[,i])<0) A[,i] = -A[,i] }
h<-diag(A%*%t(A))
rowname<-c("Variance","Proportion","Cumulative")
B<-matrix(0, nrow=3, ncol=p, dimnames=list(rowname, colname))
for (i in 1:p){
B[1,i]<-sum(A[,i]^2)
B[2,i]<-B[1,i]/sum_rank
B[3,i]<-sum(B[1,1:i])/sum_rank
}
W=B[2,1:m]*100;
Vars=data.frame('Variance'=B[1,],'Proportion'=B[2,]*100,
'Cumulative'=B[3,]*100)
A=A[,1:m]
if(rotation == "varimax" & m>1){
#cat("\n Factor Analysis for Princomp in Varimax: \n\n");
VA=varimax(A); A=VA$loadings;
s2=apply(A^2,2,sum);
k=rank(-s2); s2=s2[k];
W=s2/sum(B[1,])*100;
Vars=data.frame('Variance'=s2,'Proportion'=W,'Cumulative'=cumsum(W))
rownames(Vars) <- paste("Factor", 1:m, sep="")
A=A[,k]
for (i in 1:m) { if(sum(A[,i])<0) A[,i] = -A[,i] }
A=A[,1:m];
colnames(A) <- paste("Factor", 1:m, sep="")
}
fit<-NULL
fit$vars<-round(Vars[1:m,],3)
if(m<=1) fit$loadings <- data.frame("Factor1"=round(A,4))
else fit$loadings <- round(A,4)
if(common){
fit$common <- round(apply(A^2,1,sum),4)
}
Z=as.matrix(scale(X));
PCs=Z%*%solve(S)%*%A
fit$scores <- round(PCs,4)
if(ranks){
W=apply(fit$loadings^2,2,sum)
Wi=W/sum(W);
F=PCs%*%Wi
fit$ranks=data.frame(Factor=round(F,4),Rank=rank(-F))
}
fit
}
msa.KMO<-function(r){
cl <- match.call()
if (nrow(r) > ncol(r))
r <- cor(r, use = "pairwise")
Q <- try(solve(r))
if (class(Q) == as.character("try-error")) {
message("matrix is not invertible, image not found")
Q <- r
}
S2 <- diag(1/diag(Q))
IC <- S2 %*% Q %*% S2
Q <- Image <- cov2cor(Q)
diag(Q) <- 0
diag(r) <- 0
sumQ2 <- sum(Q^2)
sumr2 <- sum(r^2)
MSAi <- colSums(r^2)/(colSums(r^2) + colSums(Q^2))
kmo <- sumr2/(sumr2 + sumQ2)
ans <- list(MSAi = MSAi, KMO = kmo,result = test)
return(ans)
}
msa.bartlett<-function(R, n = NULL, diag = TRUE){
if (dim(R)[1] != dim(R)[2]) {
n <- dim(R)[1]
message("R was not square, finding R from data")
R <- cor(R, use = "pairwise")
}
p <- dim(R)[2]
if (!is.matrix(R))
R <- as.matrix(R)
if (is.null(n)) {
n <- 100
warning("n not specified, 100 used")
}
if (diag)
diag(R) <- 1
detR <- det(R)
statistic <- -log(detR) * (n - 1 - (2 * p + 5)/6)
df <- p * (p - 1)/2
pval <- pchisq(statistic, df, lower.tail = FALSE)
bartlett <- list(chisq = statistic, df = df, p.value = pval)
return(bartlett)
}
msa.cor<-function (X, diag = TRUE){
options(digits = 4)
p = ncol(X)
if (diag) {
tp = matrix(1, p, p)
for (i in 1:p) {
for (j in 1:i) tp[i, j] = cor.test(X[, i], X[, j])$stat
for (j in i:p) tp[i, j] = cor.test(X[, i], X[, j])$p.value
}
cat("corr test: \n")
tp = round(matrix(tp, p, dimnames = list(names(X), names(X))),
4)
print(tp)
cat("lower is t value, upper is p value \n")
}
else {
cat("\n corr test: t value, p value \n")
if (is.matrix(X))
var = 1:p
else var = names(X)
for (i in 1:(p - 1)) {
for (j in (i + 1):p) cat(" ", var[i], "-", var[j],
cor.test(X[,i],X[,j])$stat,cor.test(X[,i],X[,j])$p.value, "\n")
}
}
}
msa.cancor<-function (x, y, pq=min(ncol(x),ncol(y)), plot = FALSE){
x = scale(x)
y = scale(y)
n = nrow(x)
p = ncol(x)
q = ncol(y)
ca = cancor(x, y)
#cat("\n"); print(ca)
r = ca$cor
m <- length(r)
Q <- rep(0, m)
P = rep(0, m)
lambda <- 1
for (k in m:1) {
lambda <- lambda * (1 - r[k]^2)
Q[k] <- -log(lambda)
}
s <- 0
i <- m
for (k in 1:m) {
Q[k] <- (n - k + 1 - 1/2 * (p + q + 3) + s) * Q[k]
P[k] <- 1 - pchisq(Q[k], (p - k + 1) * (q - k + 1))
}
#cat("\n cancor test: \n")
#print(round(data.frame(r, Q, P),4))
cr=round(data.frame(CR=r, Q, P),4)
cat("\n")
u=as.data.frame(ca$xcoef[,1:pq]); colnames(u)=paste('u',1:pq,sep='')
#print(round(u,4))
v=as.data.frame(ca$ycoef[,1:pq]); colnames(v)=paste('v',1:pq,sep='')
#print(round(v,4))
if (plot) {
u1 = as.matrix(x) %*% u[,1]
v1 = as.matrix(y) %*% v[,1]
plot(u1, v1, xlab = "u1", ylab = "v1")
abline(lm(u1 ~ v1))
}
list(cor=cr,xcoef=t(round(u,4)),ycoef=t(round(v,4)))
}
msa.AHP<-function(B){
A=matrix(B,nrow=sqrt(length(B)),ncol=sqrt(length(B)),byrow=TRUE)
print(A)
m=ncol(A)
ai=apply(A,1,prod)^(1/m)
W=ai/sum(ai);
if(m>2){
AW=A%*%W
L_max=sum(AW/W)/m;
CI=(L_max-m)/(m-1);
RI=c(0,0,0.58,0.90,1.12,1.24,1.32,1.41,1.45,1.49,1.51)
CR=CI/RI[m]
cat('\n L_max=',L_max,'\n')
cat(' CI=',CI,'\n')
cat(' CR=',CR,'\n')
if(CR<=0.1) cat(' Consistency test is OK!\n\n')
else cat(' Please adjust the judgment matrix!\n')
}
return(W)
}
轉載請註明出處，本文鏈接：https://www.uj5u.com/qita/397546.html
標籤：AI
上一篇：淺談元宇宙
下一篇：機器學習 簡答題 速記