因子分析　市川雅教先生　オーム社

```{r}

#因子分析

data=data.frame(days=c(1:28),T1=c(-2.7,-0.3,-1.1,-1.5,-0.5,-1.4,-3.2,-2.4,-1.9,-3.7,-3.1,-1.1,1,1.9,1.7,2.4,3.5,2.9,1.4,0,1,4.4,2.6,-5.3,-4.1,-1.6,2.9,-5.4),T2=c(-3.6,-2.3,-1.8,-0.1,-0.8,-1.6,-3,-4.9,-4.4,-3.5,-3.8,-3.1,-2,-0.4,0.1,0.5,1,1.9,-0.1,-2,0.1,0.4,0.5,-5.3,-6.1,-3.3,0.3,-6.7),T3=c(-3.7,-1.7,-1.4,1.1,0,-3.7,-6.4,-7.1,-6.7,-5.6,-5.1,-7.6,-2.3,-1.9,-0.2,-0.7,1.1,2.3,0.7,-2,-1.1,0.5,-1.2,-4.5,-5.7,-2.9,1.3,-5.7),T4=c(-6.8,-3.4,-2.1,-1.7,-1.7,-7.9,-9.8,-9.1,-6.5,-9.1,-10.8,-9.9,-6.7,-4.4,-2.8,-3.2,-1.5,0,-0.8,-4.3,-2.3,-1.6,-3.1,-3.1,-7.1,-5.3,-2.5,-5.2),T5=c(-7.6,-7.6,-6.2,-0.8,-4.2,-8.2,-13.3,-12.9,-16.9,-15.1,-15.1,-16.1,-11.1,-6.5,-6.9,-5.5,-4.3,-2.4,-4.8,-9.6,-6.1,-5.2,-10.2,-7.1,-8.1,-6.4,-1.9,-5.1))


ave=apply(data,2,mean)[2:ncol(data)]

sd=apply(data,2,sd)[2:ncol(data)]

R=cor(as.matrix(data[2:ncol(data)]))

diag(R)<-0

max_R=apply(R,1,max)

diag(R)<-max_R

mu=sum(R)

a1=apply(R,1,sum)/sqrt(mu)

r1=R-t(t(a1))%*%a1

r1[1,]=r1[1,]*(-1)

r1[2,]=r1[2,]*(-1)

r1[,1]=r1[,1]*(-1)

r1[,2]=r1[,2]*(-1)

diag(r1)=0

diag(r1)=apply(r1,1,max)

r2=r1

a2=apply(r2,1,sum)/sqrt(sum(r2))

a2[1]=a2[1]*(-1);a2[2]=a2[2]*(-1)

r3=r2-t(t(a2))%*%a2

diag(r3)=0

diag(r3)<-apply(r3,1,max)

a3=apply(r3,1,sum)/sqrt(sum(r3))

fai=sum(abs(r3))/sum(abs(r2))

#principal factor analysis

A=cbind(a1,a2)

mat=matrix(sum(a1*a2),ncol=2,nrow=2)

diag(mat)=apply(A^2,2,sum)

lambda1=eigen(mat)$vectors[,1]/sum(eigen(mat)$vectors[,1]^2)

(mat-diag(eigen(mat)$values[1],2))%*%lambda1

lambda2=eigen(mat)$vectors[,2]/sum(eigen(mat)$vectors[,2]^2)

(mat-diag(eigen(mat)$values[2],2))%*%lambda2

lambda=cbind(lambda1,lambda2)

A%*%lambda


```

```{r}

#主因子法

library(dplyr)

n=100;m=20

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

for(j in 1:ncol(X)){
  
  X[,j]=rpois(n,sample(20:50,1))
  
}

cov_mat=cov(X)

eigen_values=eigen(cov_mat)$values

eigen_vectors=eigen(cov_mat)$vector

Ha=eigen_vectors

for(j in 1:ncol(eigen_vectors)){
  
Ha[,j]=Ha[,j]*sqrt(eigen_values)[j]  
  
}

m=17

mu=apply(X,2,mean)

ite=10000

pthi_hat=diag(0.1,ncol(cov_mat))

cost=10;cost_pre=11

cost_vec=c()

for(j in 1:ite){
  
if(sum(pthi_hat>0.0001)==ncol(cov_mat)){  

eigen_values=eigen(cov_mat-pthi_hat)$values    
  
eigen_vectors=eigen(cov_mat-pthi_hat)$vector

tueta=diag(sqrt(eigen_values[1:m]))

omega=eigen_vectors[,1:m]

lambda=omega%*%tueta

H_hat=omega%*%(tueta^2)%*%t(omega)

eigen(H_hat)$values

pthi=cov_mat-H_hat

pthi_hat=diag(diag(cov_mat-lambda%*%t(lambda)))

#print(paste0("times:",j))

#print(sum(diag(pthi_hat)))

cost=sum(abs(cov_mat-(lambda%*%t(lambda)+pthi_hat)))

cost_vec=c(cost_vec,cost)

#print(cost)

}else{
  
print("stop")  
  
print(paste0("cost:",sum(abs(cov_mat-(lambda%*%t(lambda)+pthi_hat)))))  
  
}
  
}  

eta2=diag(lambda%*%t(lambda))

pthi_vec=diag(pthi_hat)

sigma=eta2+pthi_vec

common_factor=(eta2/(eta2+(pthi_vec)))

unique_factor=(pthi_vec/(eta2+(pthi_vec)))

mu+lambda%*%common_factor[1:m]



```

```{r}

#主因子法

#ためしにirisで因子分析

library(dplyr)



X=data.frame(iris)

species=X[,ncol(X)]

X=X[,-ncol(X)]

for(j in 1:ncol(X)){
  
X[,j]=as.numeric(X[,j])  
  
}

X=as.matrix(X)

n=nrow(X);m=ncol(X)


cov_mat=cov(X)

eigen_values=eigen(cov_mat)$values

eigen_vectors=eigen(cov_mat)$vector

Ha=eigen_vectors

for(j in 1:ncol(eigen_vectors)){
  
Ha[,j]=Ha[,j]*sqrt(eigen_values)[j]  
  
}

mu=apply(X,2,mean)

m=2

mu=apply(X,2,mean)

ite=10000

pthi_hat=diag(0.1,ncol(cov_mat))

cost=10;cost_pre=11

cost_vec=c()

for(j in 1:ite){
  
if(sum(pthi_hat>0.0001)==ncol(cov_mat)){  

eigen_values=eigen(cov_mat-pthi_hat)$values    
  
eigen_vectors=eigen(cov_mat-pthi_hat)$vector

tueta=diag(sqrt(eigen_values[1:m]))

omega=eigen_vectors[,1:m]

lambda=omega%*%tueta

H_hat=omega%*%(tueta^2)%*%t(omega)

eigen(H_hat)$values

pthi=cov_mat-H_hat

pthi_hat=diag(diag(cov_mat-lambda%*%t(lambda)))

#print(paste0("times:",j))

#print(sum(diag(pthi_hat)))

cost=sum(abs(cov_mat-(lambda%*%t(lambda)+pthi_hat)))

cost_vec=c(cost_vec,cost)

#print(cost)

}else{
  
print("stop")  
  
print(paste0("cost:",sum(abs(cov_mat-(lambda%*%t(lambda)+pthi_hat)))))  
  
}
  
}  

eta2=diag(lambda%*%t(lambda))

pthi_vec=diag(pthi_hat)

sigma=eta2+pthi_vec

common_factor=(eta2/(eta2+(pthi_vec)))

unique_factor=(pthi_vec/(eta2+(pthi_vec)))

mu+lambda%*%common_factor[1:m]



```

```{r}

#イメージ理論(Guttman(1953))

n=100;m=20

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

for(j in 1:m){
  
X[,j]=rpois(n,sample(1:50,1)) 
  
}

cov_mat=cov(X)

X=rpois(m,sample(1:50,1))

X=X-mean(X)

sigma11=cov_mat[1,1]

sigma12=cov_mat[1,2:ncol(cov_mat)]

sigma21=cov_mat[2:nrow(cov_mat),1]

sigma22=cov_mat[2:nrow(cov_mat),2:ncol(cov_mat)]

alpha=solve(sigma22)%*%sigma21

g1=t(alpha)%*%t(t(X[2:length(X)]))

var_alpha_x2=t(alpha)%*%sigma21

cor_coef2=t(alpha)%*%sigma22%*%alpha/sigma11

cor_coef=sqrt(cor_coef2)

R2=1-1/(diag(cov_mat)*diag(solve(cov_mat)))

R=sqrt(R2)

delta=diag(1/(diag(solve(cov_mat))))

u=delta%*%solve(cov_mat)%*%X

g_vec=(diag(1,m)-delta%*%solve(cov_mat))%*%X

var_g=cov_mat+delta%*%solve(cov_mat)%*%delta-2*delta

var_u=delta%*%solve(cov_mat)%*%delta

cov_g_u=delta-var_u

diag(delta-var_u)

var_x=var_g-var_u+2*delta

R2*diag(cov_mat)

(1-R2)*diag(cov_mat)

result=data.frame(g=g_vec,u=u,var_g=diag(var_g),var_u=diag(var_u),R=R)

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

for(j in 1:m){
for(i in 1:m){  
 
p_cor[i,j]=solve(cov_mat)[i,j]/sqrt(solve(cov_mat)[i,i]*solve(cov_mat)[j,j])  
  
}
}

cov_u=p_cor

for(j in 1:m){
for(i in 1:m){  
 
cov_u[i,j]=cov_u[i,j]*sqrt(diag(var_u)[i]*diag(var_u)[j])  
  
}
}


```

```{r}

#alpha 因子分析

library(dplyr);

vec=sample(1:100,5*5,replace = T)

S=t(t(vec))%*%vec

pthi=diag(0.1,ncol(S))

ite=10000

D_H_vec=c()

for(j in 1:ite){

D_H=diag(diag(S-pthi),ncol(S))

mat=solve(D_H)%*%(S-pthi)%*%solve(D_H)

omega=eigen(mat)$vector

sita=diag(eigen(mat)$values,ncol(S))

sita[(sita)<10^(-3)]=0

lambda=D_H%*%omega%*%sqrt(sita)

pthi=diag(diag(S-lambda%*%t(lambda)),ncol(S))

print(sum(abs(diag(D_H))))

D_H_vec=c(D_H_vec,sum(abs(diag(D_H))))

}

plot(c(1:length(D_H_vec)),D_H_vec,type="o",col=2)

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

delta=t(alpha)%*%(D_H^2)%*%alpha


```
```{r}

#正準因子分析

rm(list=ls())

n=50;m=20

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

for(j in 1:m){

X[,j]=rpois(n,sample(1:50,1)) 

}

S=cov(X)

L=t(chol(S))

W=solve(L)

pthi=diag(rep(0.1,m))

mat=W%*%pthi%*%t(W)

M=m

m=7;d=sum((m-M)^2-(m+M))/2

ite=200

if(d>0){

vec=c()

for(j in 1:ite){

N=diag(sort(eigen(mat)$values[(ncol(mat)-m+1):ncol(mat)]))

N[N>1]=1

Z=eigen(mat)$vector[,(ncol(mat)-m+1):ncol(mat)]

lambda=L%*%Z%*%sqrt(diag(1,m)-(N))

pthi=diag(diag(S-lambda%*%t(lambda)))

print(paste0("times:",j))

print(mean(abs(S-(lambda%*%t(lambda)+pthi))))

mat=W%*%pthi%*%t(W) 

vec=c(vec,mean(abs(S-(lambda%*%t(lambda)+pthi))))

}

ite_data=data.frame(times=c(1:length(vec)),value=vec)

plot(c(1:length(vec)),vec,type="p")

eta2=diag(S)-diag(pthi)

#p.85 likelihood ratio test

estimate=(lambda%*%t(lambda)+pthi)

Y_m=(n-(2*m+4*M+5)/6)*(-log(det(S%*%solve(estimate))))

d=sum((m-M)^2-(m+M))/2

qchisq(1-0.05,d)

print(paste0("kai2:",qchisq(1-0.05,d),",test-statics-value:",Y_m))

F_hat=-log(det(S%*%solve(estimate)))

#近似値 p.89

F_hat_value=(sum((S-estimate)^2/(t(t(diag(pthi)))%*%t(diag(pthi))))-sum(diag((S-estimate)^2/(t(t(diag(pthi)))%*%t(diag(pthi))))))/2

#AIC誤差による検定

p_star=M*(M+1)/2
  
pthi_p_k=M*(m+1)-m*(m-1)/2

#AIC誤差による統計検定量

AIC_df=abs(Y_m-2*(p_star-pthi_p_k))

print(paste0("kai2:",qchisq(1-0.01,floor(p_star-pthi_p_k)),",test-statics-value:",AIC_df))

qchisq(1-0.01,floor(p_star-pthi_p_k))

GFI=1-sum(diag((S-estimate)%*%solve(estimate)))^2/(sum(diag(S%*%solve(estimate)))^2)

AGFI=1-(p_star/d)*(1-GFI)

X=matrix(sample(1:10,30*M,replace=T),nrow=M)

Y=(estimate%*%solve(S))%*%X

}

```

```{r}
#直交回転法 p.137

p=10;m=6

A=matrix(sample(1:10,p*m,replace=T),ncol=m)

w=1

M=diag(1,p)-(w/p)*t(t(rep(1,p)))%*%t(rep(1,p))

lambda=A

ite=100

ORMAX_vec=c()

CFMIN=c()

for(j in 1:ite){

B=lambda*lambda

W_mat=lambda*(M%*%B)

U=svd(t(A)%*%W_mat)$u

D=svd(t(A)%*%W_mat)$d

V=svd(t(A)%*%W_mat)$v

t=U%*%t(V)

lambda_pre=lambda

lambda=A%*%t

#print(sum(abs(lambda-lambda_pre)))

print(paste0("ORMAX:",sum(lambda^4)-(w/p)*sum(apply(lambda^2,2,sum)^2)))


ORMAX_vec=c(ORMAX_vec,sum(lambda^4)-(w/p)*sum(apply(lambda^2,2,sum)^2))

eta=apply(lambda^2,1,sum)

CFMIN=c(CFMIN,-(sum(lambda^4)-(w/p)*sum(apply(lambda^2,2,sum)^2))+(1-w/p)*sum(eta^4))

}



plot(ORMAX_vec,main="ORMAX")

plot(CFMIN,main="CFMIN")

D1=diag(sqrt(diag(solve(t(t)%*%t))))

pthi1=solve(D1)%*%solve(t(t)%*%t)%*%solve(D1)
```

```{r}

#因子分析-その理論と方法　柳井晴夫等

#PARAFAC model p.183

n=20;p=5;m=6

X1=matrix(sample(1:10,n*p,replace=T),ncol=p)

X2=matrix(sample(1:30,n*p,replace=T),ncol=p)

X3=matrix(sample(1:60,n*p,replace=T),ncol=p)

W1=diag(sample(1:10,m,replace=T),m);W2=diag(sample(20:30,m,replace=T),m);W3=diag(1,m)

f=matrix(sample(1,n*m,replace=T),ncol=m)

lambda=matrix(sample(1:20,m*p,replace=T),nrow=p)

ite=10000

for(j in 1:ite){
  
val1=X1%*%lambda%*%W1+X2%*%lambda%*%W2+X3%*%lambda%*%W3

val2=(W1%*%t(lambda)%*%lambda%*%W1+W2%*%t(lambda)%*%lambda%*%W2+W3%*%t(lambda)%*%lambda%*%W3)

f=val1%*%solve(val2)

lambda=(t(X1)%*%f%*%W1+t(X2)%*%f%*%W2+t(X3)%*%f%*%W3)%*%solve(W1%*%t(f)%*%f%*%W1+W2%*%t(f)%*%f%*%W2+t(W3)%*%t(f)%*%f%*%W3)

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

for(s in 1:m){
for(l in 1:m){  
 
t[s,l]=sum(lambda[,s]*lambda[,l])*sum(f[,s]*f[,l])   
  
}
}  

h1=apply(t(X1)%*%f*lambda,2,sum)

h2=apply(t(X2)%*%f*lambda,2,sum)

h3=apply(t(X3)%*%f*lambda,2,sum)

W1=diag(c(solve(t)%*%h1))

W2=diag(c(solve(t)%*%h2))

W3=diag(c(solve(t)%*%h3))

cost=sum(diag(t(X1-f%*%W1%*%t(lambda))%*%(X1-f%*%W1%*%t(lambda))+t(X2-f%*%W2%*%t(lambda))%*%(X2-f%*%W2%*%t(lambda))+t(X3-f%*%W3%*%t(lambda))%*%(X3-f%*%W3%*%t(lambda))))

print(cost)

}


```
```{r}

#generalized latent variable modeling p.177

#探索的因子分析

n=20;m=10;N=50

Y=array(0,dim=c(n,N))

lambda=array(1,dim=c(n,N))

sita=diag(1,n)

for(j in 1:N){
  
Y[,j]=rpois(n,sample(30:40,1))
  
}


S_yy=array(0,dim=c(n,n))

for(j in 1:N){
  
S_yy=S_yy+Y[,j]%*%t(Y[,j])/N  

}


ite=10000

for(j in 1:ite){
  
sita_pre=sita  
  
f=t(lambda)%*%solve(lambda%*%t(lambda)+sita)  
  
gamma=diag(1,N)-f%*%lambda

lambda=S_yy%*%t(f)%*%solve(f%*%S_yy%*%t(f)+gamma)

sita=diag(diag(S_yy-lambda%*%f%*%S_yy))

print(sum(abs(sita_pre-sita)))
  
}

eta=array(0,dim=c(N,N))

for(j in 1:N){
  
eta[,j]=f%*%Y[,j]  
  
}

predict=lambda%*%eta

```

```{r}

#LISREL

m=10;l=12;n=6

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

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

for(j in 1:nrow(X)){
  
X[j,]=rpois(n,sample(10:20,1))  

Y[j,]=rpois(n,sample(30:40,1))
  
}

S=cbind(rbind(cov((X),(X)),(cov((X),(Y)))),rbind((cov((Y),(X))),cov((Y),(Y))))

det(S)

lambda_y=t(array(1,dim=c(m,n)))

lambda_x=t(array(1,dim=c(m,n)))

B=diag(1,m);A=solve(B)

gamma=diag(1,m)

sita_epsilon=diag(1,n);sita_delta=diag(1,n)

D=solve(B)%*%gamma

pthi=diag(1,m);tueta=diag(1,m);sigma_epsilon=diag(1,m)

sigma_yy=lambda_y%*%(solve(B)%*%(gamma%*%tueta%*%t(gamma)+pthi)%*%solve(t(B)))%*%t(lambda_y)+sita_epsilon

sigma_yx=lambda_y%*%solve(B)%*%gamma%*%tueta%*%t(lambda_x)+sita_delta

sigma_xy=lambda_x%*%tueta%*%t(gamma)%*%solve(t(B))%*%t(lambda_y)

sigma_xx=lambda_x%*%tueta%*%t(lambda_x)+sita_delta

sigma=cbind(rbind(sigma_yy,sigma_xy),rbind(sigma_yx,sigma_xx))

omega=solve(sigma)%*%(sigma-S)%*%solve(sigma)

omega_yy=omega[1:n,1:n];omega_yx=omega[1:n,(n+1):ncol(omega)]

omega_xy=omega[(n+1):nrow(omega),1:n];omega_xx=omega[(n+1):nrow(omega),(n+1):ncol(omega)]

C=D%*%tueta%*%t(D)+A%*%pthi%*%t(A)


#iteration

eta=-0.0001

ite=1000000

f=c()

f_val=10^(7);f_val_pre=10^7+1

for(j in 1:ite){
  
if(f_val_pre-f_val>0){
  
f_val_pre=f_val  
  
lambda_y=lambda_y+eta*(omega_yy%*%lambda_y%*%C+t(omega_xy)%*%lambda_x%*%pthi%*%t(D))  

lambda_x=lambda_x+eta*(omega_xy%*%lambda_y%*%D%*%pthi+omega_xx%*%lambda_x%*%pthi)  
  
B=B-eta*(t(A)%*%t(lambda_y)%*%(omega_yy%*%lambda_y%*%C+t(omega_xy)%*%lambda_x%*%pthi%*%t(D)))

gamma=gamma+eta*(t(A)%*%t(lambda_y)%*%(omega_yy%*%lambda_y%*%D+t(omega_xy)%*%lambda_x)%*%pthi)

pthi=pthi+eta*(t(A)%*%t(lambda_y)%*%omega_yy%*%lambda_y%*%A)
  
tueta=tueta+eta*(t(D)%*%t(lambda_y)%*%omega_yy%*%lambda_y%*%D+t(lambda_x)%*%omega_xy%*%lambda_y%*%D+t(D)%*%t(lambda_y)%*%t(omega_xy)%*%lambda_x+t(lambda_x)%*%omega_xx%*%lambda_x)

sita_epsilon=sita_epsilon+eta*diag(diag(omega_yy))

sita_delta=sita_delta+eta*diag(diag(omega_xy))

sigma_yy=lambda_y%*%(solve(B)%*%(gamma%*%tueta%*%t(gamma)+pthi)%*%solve(t(B)))%*%t(lambda_y)+sita_epsilon

sigma_yx=lambda_y%*%solve(B)%*%gamma%*%tueta%*%t(lambda_x)+sita_delta

sigma_xy=lambda_x%*%tueta%*%t(gamma)%*%solve(t(B))%*%t(lambda_y)

sigma_xx=lambda_x%*%tueta%*%t(lambda_x)+sita_delta

sigma=cbind(rbind(sigma_yy,sigma_xy),rbind(sigma_yx,sigma_xx))

omega=solve(sigma)%*%(sigma-S)%*%solve(sigma)

omega_yy=omega[1:n,1:n];omega_yx=omega[1:n,(n+1):ncol(omega)]

omega_xy=omega[(n+1):nrow(omega),1:n];omega_xx=omega[(n+1):nrow(omega),(n+1):ncol(omega)]

A=solve(B);D=solve(B)%*%gamma;C=D%*%tueta%*%t(D)+A%*%pthi%*%t(A)


f_val=sum((S-sigma)^2)

f_val=as.numeric(f_val)

print(f_val)

f=c(f,f_val)

R_x=1-det(sita_delta)/det(sigma_xx)

R_y=1-det(sita_epsilon)/det(sigma_yy)

}else{
  
print("stop")  
  
}

}

```