因子分析法(Factor Analysis as a Statistical Method, DN Lawley, A E Maxwell 訳本）

```{r}


library(dplyr);library(stringr)

#因子分析法　丘本正監修　p.31~

#最尤法による因子分析　

#相関行列
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=50

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

U=diag(rep(0.2,n))

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


for(l in 1:ite){

sqrt_U=diag(1/sqrt(diag(U)))

omega=eigen(sqrt_U%*%R%*%sqrt_U)$vectors[,1:m]  

sita=eigen(sqrt_U%*%R%*%sqrt_U)$values[1:m]

P=t(diag(sqrt(diag(sita-diag(rep(1,m)))))%*%t(omega)%*%diag(sqrt(diag(U))))

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

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

L=-N*log(det(V))/2-N*sum(R*solve(V))/2

J=t(P)%*%solve(U)%*%P

H=t(P)%*%solve(U)%*%(R-U)%*%solve(U)%*%P

value=sum(solve(diag(rep(1,nrow(J)))+J)*(J^2-H))

print(value)

}

N_dash=N-(2*n+5)/6-2*m/3

kai_mat=(R-V)^2/(t(t(diag(U)))%*%t(diag(U)))

kai2=N_dash*sum(kai_mat[lower.tri(kai_mat)])

qchisq(1-0.05,df=((n-m)^2-(n+m))/2)



```

```{r}

#Factor Analysis as a Statistical Method, DN Lawley, A E Maxwell 

#6.4 The Promax Method p.77

dim=3

Q=abs(P^(dim-1))*P

U=solve(t(P)%*%P)%*%t(P)%*%Q

D=diag(sqrt(diag(solve(t(U)%*%U))))

M=U%*%D

pthi=solve(t(M)%*%M)

cor(c(Q),c(P%*%U))

```

```{r}

#Quantimax p.82

L=P

ite=10^5

for(l in 1:ite){

C=array(0,dim=c(m,m))

D=array(0,dim=c(m,m))

Q=array(0,dim=c(m,m))

for(r in 1:m){
for(s in 1:m){
  
C[r,s]=4*sum(L[,r]*L[,s]*(L[,r]^2-L[,s]^2))  

D[r,s]=sum((L[,r]^2-L[,s]^2)^2-(2*L[,r]*L[,s])^2)

Q[r,s]=sum(L[,r]^4+L[,s]^4)
  
}}

print(sum(log(diag(Q))))

pthi=abs(atan(C/D)/4)*sign(C)



t=diag(rep(1,m))

for(r in 2:(m)){
  
t_sub=diag(rep(1,m))  
  
kaiten=array(0,dim=c(2,2))

kaiten[1,1]=cos(pthi[1,r])
kaiten[1,2]=-sin(pthi[1,r])
kaiten[2,1]=sin(pthi[1,r])
kaiten[2,2]=cos(pthi[1,r])

t_sub[(r-1):(r),(r-1):(r)]=kaiten

t=t%*%t_sub
  
}

L=L%*%t

}

```

```{r}

#Varimax p.85

L=P

L_pre=L

h=sqrt(diag(P%*%t(P)))

ite=10^5

for(l in 1:ite){

A=array(0,dim=c(m,m))

B=array(0,dim=c(m,m))
  
C=array(0,dim=c(m,m))

D=array(0,dim=c(m,m))

V=array(0,dim=c(m,m))

for(r in 1:m){
for(s in 1:m){
  
u=L_pre[,r]/h[r];v=L_pre[,s]/h[s];w=L[,r]/h[r];z=L[,s]/h[s]  
 
A[r,s]=sum(u^2-v^2)
  
B[r,s]=2*sum(u*v)
 
C[r,s]=4*sum(w*z*(w^2-z^2))  

D[r,s]=sum((w^2-z^2)^2-(2*w*z)^2)

V[r,s]=sum(w^4+z^4)-(sum(w^2)^2+sum(z^2)^2)/n
  
}}

print(sum((V)))

pthi=(atan((C-2*A*B/n)/(D-(A^2-B^2)/n))/4)

L_pre=L

t=diag(rep(1,m))

for(r in 2:(m)){
  
t_sub=diag(rep(1,m))  
  
kaiten=array(0,dim=c(2,2))

kaiten[1,1]=cos(pthi[1,r])
kaiten[1,2]=-sin(pthi[1,r])
kaiten[2,1]=sin(pthi[1,r])
kaiten[2,2]=cos(pthi[1,r])

t_sub[(r-1):(r),(r-1):(r)]=kaiten

t=t%*%t_sub
  
}

L=L%*%t

}



```

```{r}

#Factor Analysis as a Statistical Method, DN Lawley, A E Maxwell 

#The Fletcher and Powell Method p.142

g=function(x){
  
x1=x[1];x2=x[2]

f1=x1^2-3*x1

f2=x2^3-4*x2

return(c(f1,f2))

}

sita=rep(10,2)

E=diag(rep(1,2))

ite=3*10^3

sita_pre=rep(0,2)

g_pre=rep(0,2)

eta=0.01

for(l in 1:ite){
  
y=sita-sita_pre

h=g(sita)-g(sita_pre)
  
sita_pre=sita

g_pre=g(sita_pre)
  
z=E%*%h

beta=sum(y*h)

r=sum(z*h)

E=E+t(t(c(y)))%*%t(c(y))/beta-t(t(c(z)))%*%t(c(z))/r
  
sita=sita-eta*E%*%g(sita)

print(g(sita))

}



```

```{r}

#最尤法による因子分析

#Factor Analysis as a Statistical Method, DN Lawley, A E Maxwell 

#4.3 The Iterative Procedure  p.30

R=matrix(c(19.94,10.5,6.59,8.63,10.5,23.56,19.71,7.97,6.59,19.71,20.95,3.93,8.63,7.97,3.93,7.55),ncol=4,nrow=4)

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

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

U=diag(rep(0.2,n))

#The Fletcher and Powell Method p.142

g=function(x,yy){
  
V=diag(c(x))

dx=solve(V)%*%(yy%*%t(yy)+V-R)%*%solve(V)

return(diag(dx))

}

sita=rep(10,ncol(U))

E=diag(rep(1,ncol(U)))

ite=3*10^5

sita_pre=rep(0.01,ncol(U))

g_pre=rep(0,ncol(U))

eta=10^(-3)

for(l in 1:ite){
  
S_star=solve(diag(sqrt(sita)))%*%R%*%solve(diag(sqrt(sita)))

eig_value=eigen(S_star)$values

eig_value=eig_value[eig_value>1]

eig_vector=eigen(S_star)$vectors[,1:length(eig_value)]

P=diag(sqrt(sita))%*%(eig_vector)%*%diag(sqrt(eig_value-1))

y=sita-sita_pre

h=g(sita,P)-g(sita_pre,P)
  
sita_pre=sita

g_pre=g(sita_pre,P)
  
z=E%*%h

beta=sum(y*h)

r=sum(z*h)

E=E+t(t(c(y)))%*%t(c(y))/beta-t(t(c(z)))%*%t(c(z))/r
  
sita=sita-eta*E%*%g(sita,P)

sita=c(sita)

print(sum(abs(g(sita,P))))

}



```