行動科学における因子分析法（松浦義之先生著）

#正準相関係数　p.75

R=matrix(c(1,0.859,0.473,0.301,0.859,1,0.436,0.327,0.473,0.436,1,0.73,0.301,0.327,0.73,1),ncol=4)

Vxx=R[1:2,1:2];Vxy=R[1:2,3:4];Vyy=R[3:4,3:4]

eig_values=eigen(solve(Vxx)%*%Vxy%*%solve(Vyy)%*%t(Vxy))$values

eig_vectors=eigen(solve(Vxx)%*%Vxy%*%solve(Vyy)%*%t(Vxy))$vectors

for(i in 1:ncol(eig_vectors)){

eig_vectors[,i]=eig_vectors[,i]/sqrt(apply(eig_vectors^2,2,sum)[i])

}

mu=eig_values[1];lam=eig_values[2]

a=eig_vectors[,1]

b=solve(Vyy)%*%t(Vxy)%*%eig_vectors[,1]/sqrt(mu)

lambda=prod(1-eig_values)

p=ncol(Vxx);q=ncol(Vyy)

m=154

kai2=-(m-(p+q+1)/2)*log(lambda)

qchisq(1-0.01,df=p*q)

#多群因子解析法 p.144

library(dplyr)

#相関行列
R=matrix(c(1,0.846,0.805,0.859,0.473,0.398,0.301,0.382,0.246,0.231,0.146,0,1,0.881,0.826,0.376,0.326,0.277,0.415,0.315,0.326,0.213,0,0,1,0.801,0.386,0.319,0.237,0.345,0.3,0.293,0.314,0,0,0,1,0.436,0.329,0.327,0.365,0.421,0.411,0.213,0,0,0,0,1,0.762,0.730,0.629,0.363,0.321,0.137,0,0,0,0,0,1,0.583,0.577,0.414,0.365,0.195,0,0,0,0,0,0,1,0.539,0.321,0.273,0.332,0,0,0,0,0,0,0,1,0.325,0.421,0.336,0,0,0,0,0,0,0,0,1,0.821,0.754,0,0,0,0,0,0,0,0,0,1,0.832,0,0,0,0,0,0,0,0,0,0,1),ncol=11)

R=R+t(R);diag(R)=diag(R)/2

mat=R

#共通性を対角成分に入れて計算する。各行の最大値で置き換えてもいい。
diag(mat)=c(0.848,0.877,0.821,0.888,0.871,0.654,0.602,0.54,0.761,0.892,0.795)

#多群のクラスを設定する
dat=data.frame(class=c(rep(1,4),rep(2,4),rep(3,3)),mat)

#相関係数の和の行列
summary=array(0,dim=c(3,ncol(mat)))

for(j in 1:nrow(summary)){

dat_sub=dat %>% filter(class==j) %>% select(-class)  

summary[j,]=c(apply(dat_sub,2,sum))  

}

S_dat=data.frame(class=dat$class,t(as.matrix(summary)))

Spq=array(0,dim=c(nrow(summary),nrow(summary)))

for(j in 1:nrow(summary)){

S_dat_sub=S_dat%>%filter(class==j)%>%select(-class)

Spq[j,]=c(apply(S_dat_sub,2,sum))  

}

rt=array(0,dim=dim(Spq))

for(i in 1:nrow(rt)){
for(j in 1:ncol(rt)){

rt[i,j]=Spq[i,j]/(sqrt(Spq[i,i])*sqrt(Spq[j,j]))  

}  
}

pthi_inv=solve(rt)


#因子構造行列
S=t(summary/diag(sqrt(Spq)))

#因子パターン行列
P=S%*%pthi_inv

#再生相関行列
R_star=P%*%rt%*%t(P)

R_star2=R_star;R2=R

diag(R_star2)=0;diag(R2)=0

#残差行列
R_bar=R2-R_star2

#残差行列の最大値
max(abs(R_bar))

#3.4.6 最尤度因子解法　p.199

#相関行列
R=matrix(c(1,0.846,0.805,0.859,0.473,0.398,0.301,0.382,0.246,0.231,0.146,0,1,0.881,0.826,0.376,0.326,0.277,0.415,0.315,0.326,0.213,0,0,1,0.801,0.386,0.319,0.237,0.345,0.3,0.293,0.314,0,0,0,1,0.436,0.329,0.327,0.365,0.421,0.411,0.213,0,0,0,0,1,0.762,0.730,0.629,0.363,0.321,0.137,0,0,0,0,0,1,0.583,0.577,0.414,0.365,0.195,0,0,0,0,0,0,1,0.539,0.321,0.273,0.332,0,0,0,0,0,0,0,1,0.325,0.421,0.336,0,0,0,0,0,0,0,0,1,0.821,0.754,0,0,0,0,0,0,0,0,0,1,0.832,0,0,0,0,0,0,0,0,0,0,1),ncol=11)

R=R+t(R);diag(R)=diag(R)/2

n=nrow(R);m=3

#標本の大きさ
N=100

#勾配上昇法

#初期パラメータの設定
l=cumsum(c(3,3,5))

P=array(0,dim=c(n,m))

for(j in 1:m){

if(j==1){  

P[1:l[1],j]=0.5  

}else{

P[(l[j-1]+1):l[j],j]=0.5  

}

}

U=diag(rep(1,n))

#学習率
eta=0.001

#反復計算回数
ite=10^4


for(l in 1:ite){

V=P%*%t(P)+U

dP=t(-n*t(P)%*%solve(V)+n*t(P)%*%solve(V)%*%R%*%solve(V))

dU=-n*diag(solve(V)-solve(V)%*%R%*%solve(V))/2

#共通性因子パターン行列の更新
P=P+eta*dP

#独立性因子行列の更新
U=U+eta*diag(dU)

#対数尤度関数
L=-N*log(det(V))/2-N*sum(R*solve(V))/2

print(L)

}

library(dplyr);library(stringr)

#(サンプル) alpha因子分析　

#相関行列
R=matrix(c(1,0.846,0.805,0.859,0.473,0.398,0.301,0.382,0.246,0.231,0.146,0,1,0.881,0.826,0.376,0.326,0.277,0.415,0.315,0.326,0.213,0,0,1,0.801,0.386,0.319,0.237,0.345,0.3,0.293,0.314,0,0,0,1,0.436,0.329,0.327,0.365,0.421,0.411,0.213,0,0,0,0,1,0.762,0.730,0.629,0.363,0.321,0.137,0,0,0,0,0,1,0.583,0.577,0.414,0.365,0.195,0,0,0,0,0,0,1,0.539,0.321,0.273,0.332,0,0,0,0,0,0,0,1,0.325,0.421,0.336,0,0,0,0,0,0,0,0,1,0.821,0.754,0,0,0,0,0,0,0,0,0,1,0.832,0,0,0,0,0,0,0,0,0,0,1),ncol=11)

R=R+t(R);diag(R)=diag(R)/2

n=nrow(R);m=4;N=50

#勾配上昇法

#初期パラメータの設定
l=cumsum(c(3,3,2,3))

P=array(0,dim=c(n,m))

for(j in 1:m){

if(j==1){  

P[1:l[1],j]=0.5  

}else{

P[(l[j-1]+1):l[j],j]=0.5  

}

}

U=diag(rep(0.2,n))

#反復計算回数
ite=1000


for(l in 1:ite){

DH=diag(diag(R-U))

sita=diag(eigen(solve(DH)%*%(R-U)%*%solve(DH))$values[1:m])

omega=eigen(solve(DH)%*%(R-U)%*%solve(DH))$vectors[,1:m]  

P=DH%*%omega%*%diag(sqrt(diag(sita)))

U=diag(diag(R-P%*%t(P)))

V=P%*%t(P)+U

alpha=(1/sqrt(diag(sita)[1]))*solve(DH)%*%omega[,1]

coef=(m/(m-1))*(1-as.numeric(t(alpha)%*%(DH^2)%*%alpha))

print(coef)

}