非計量多変量解析法ー主成分分析から多重対応分析へ 足立浩平・村上隆先生

#非計量モデル

n=5;m=4

X=matrix(sample(1:10,n*m,replace=T),ncol=m)

J=diag(1,n)-t(t(rep(1,n)))%*%t(rep(1,n))/n

diag(sqrt(diag(solve(t(X)%*%J%*%X))))

Z=X

K=svd(Z)$u;lambda=diag(svd(Z)$d);L=svd(Z)$v

t=eigen(t(X)%*%X)$vector

f=sqrt(n)*K%*%diag(1,ncol(K))%*%diag(rep(1,ncol(K)))

A=L%*%lambda%*%diag(1,ncol(K))%*%diag(rep(1,ncol(K)))/sqrt(n)

f%*%t(A)

Z



#ムーアペンローズ行列

X_plus=L%*%solve(lambda)%*%t(K)

#ムーアペンローズ行列の性質の確認

#Xに等しい
X%*%X_plus%*%X

#X_plusに等しい
X_plus%*%X%*%X_plus

#X%*%X_plusに等しい
t(X%*%X_plus)

#X_plus%*%Xに等しい
t(X_plus%*%X)

p=2

K_p=K[,1:p];K_rp=K[,(p+1):ncol(K)]

L_p=L[,1:p];L_rp=L[,(p+1):ncol(L)]

t(K_p)%*%K

t(L_p)%*%L

t=diag(1,p)

V=L_p%*%t

g_V=sum(diag(t(V)%*%t(X)%*%X%*%V))


#2.5近似行列の分解

f=K_p%*%t*sqrt(n)

A=L_p%*%lambda[1:p,1:p]%*%t/sqrt(n)


#Xの近似行列
f%*%t(A)

Y=K_p%*%lambda[1:p,1:p]%*%t(L_p)

f2=t((t(f)-c(apply(f,2,mean)))/c(apply(f,2,sd)))

X_p=X[,1:p];X_rp=X[,(p+1):ncol(X)]

D_X=array(0,dim=c(ncol(X),ncol(X)))

D_X[1:p,1:p]=t(X_p)%*%(X_p)

D_X[(p+1):ncol(X),(p+1):ncol(X)]=t(X_rp)%*%(X_rp)

#性質2.3 p.23

N=svd(X%*%solve(sqrt(D_X)))$u[,1:p]

pthi=diag(svd(X%*%solve(sqrt(D_X)))$d[1:p])

M=svd(X%*%solve(sqrt(D_X)))$v[,1:p]

t=diag(1,3)

f=sqrt(n)*N%*%t

W=sqrt(n)*solve(sqrt(D_X))%*%M%*%pthi%*%t

#低次元での近似

f%*%t(W)%*%(sqrt(D_X))/n

X%*%solve(sqrt(D_X))

#等質性基準に基づく多重対応分析

library(dummies)

data=data.frame(パターン=1:13,個数=c(1,3,2,5,3,2,4,3,1,1,3,1,1),無口である=c(0,0,0,0,0,0,0,1,0,1,0,0,0),生真面目である=c(1,0,0,0,1,0,0,0,0,1,0,0,0),空想を好む=c(1,1,0,0,0,0,1,0,0,0,1,0,1),活発である=c(1,0,0,1,0,0,1,0,0,0,1,1,1),冗談をいう=c(1,1,0,1,0,1,1,0,1,0,0,0,0),きれい好きである=c(0,0,1,0,1,0,0,1,1,0,0,0,0),粘り強い=c(0,0,0,0,1,0,0,1,1,1,0,1,1))

n=nrow(data)

G1=dummy(data$無口である);colnames(G1)=c("mukuti","mukutidenai")

G2=dummy(data$生真面目である);colnames(G2)=c("kimajime","kimajimedenai")

G3=dummy(data$粘り強い);colnames(G3)=c("kuusouwokonomu","kuusouwokonomanai")

G=cbind(G1,G2)

D1=t(G1)%*%G1;D2=t(G2)%*%G2;

D=array(0,dim=c(ncol(G1)+ncol(G2),ncol(G1)+ncol(G2)))

D[1:2,1:2]=D1;D[3:4,3:4]=D2;

J=diag(1,n)-t(t(rep(1,n)))%*%t(rep(1,n))/n

mat=J%*%G%*%diag(1/sqrt(diag(D)))

p=2

N=svd(mat)$u[,1:p]

pthi=diag(svd(mat)$d[1:p])

M=svd(mat)$v[,1:p]



f=sqrt(n)*N%*%diag(1,p)

W=sqrt(n)*diag(1/sqrt(diag(D)))%*%M%*%pthi%*%diag(1,p)

G%*%W

f

V=t(W)%*%D%*%W/n


#element:1

mat1=J%*%G1%*%D1

N1=svd(mat1)$u

pthi1=diag(svd(mat1)$d)

M1=svd(mat1)$v

f1=sqrt(n)*N1%*%diag(1,ncol(N1))

W1=sqrt(n)*solve(sqrt(D1))%*%M1%*%pthi1%*%diag(1,ncol(N1))

V1=t(W1)%*%D1%*%W1/n

#element:2

mat2=J%*%G2%*%D2

N2=svd(mat2)$u

pthi2=diag(svd(mat2)$d)

M2=svd(mat2)$v

f2=sqrt(n)*N2%*%diag(1,ncol(N2))

W2=sqrt(n)*solve(sqrt(D2))%*%M2%*%pthi2%*%diag(1,ncol(N2))

V2=t(W2)%*%D2%*%W2/n

#element:3

mat3=J%*%G3%*%D3

N3=svd(mat3)$u

pthi3=diag(svd(mat3)$d)

M3=svd(mat3)$v

f3=sqrt(n)*N3%*%diag(1,ncol(N3))

W3=sqrt(n)*solve(sqrt(D3))%*%M3%*%pthi3%*%diag(1,ncol(N3))

V3=t(W3)%*%D3%*%W3/n

#p.37

library(dummies)

n=5;m=4

X=matrix(sample(1:10,n*m,replace=T),ncol=m)

J=diag(1,n)-t(t(rep(1,n)))%*%t(rep(1,n))/n

diag(sqrt(diag(solve(t(X)%*%J%*%X))))

#得点行列

Z=J%*%X%*%diag(1/sqrt(diag(t(X)%*%J%*%X/n)))

K=svd(Z)$u;lambda=diag(svd(Z)$d);L=svd(Z)$v

f=sqrt(n)*K%*%diag(1,ncol(K))

A=L%*%lambda%*%diag(1,ncol(K))/sqrt(n)

f%*%t(A)

Z

#4.4 成分負荷基準に基づく非計量主成分分析 p.41 (NCA)

f0=f

colnames(X)=paste0("G_",c(1:ncol(X)))

G=rep(1,nrow(X))

for(j in 1:ncol(X)){
  
vec=dummy(X[,j])  

colnames(vec)=paste0("G",j,"_",colnames(vec))
  
G=cbind(G,vec)  
  
}

G=G[,-1]

#p.41 非計量主成分分析の交互SVDアルゴリズム

p=3

f=f0

ite=10000

eta=0.001

for(l in 1:ite){

Bq=c()

for(j in 1:ncol(X)){
  
Dj=t(G[,grep(paste0("G",j),colnames(G))])%*%G[,grep(paste0("G",j),colnames(G))]  
matj=diag(1/sqrt(diag(Dj)))%*%t(G[,grep(paste0("G",j),colnames(G))])%*%f

u=svd(matj)$u
 
sita=diag(svd(matj)$d)

v=svd(matj)$v

qj=sqrt(n)*diag(1/sqrt(diag(Dj)))%*%u[,1]

aj=sita[1]*v[,1]/sqrt(n)

Bq=c(Bq,qj)

}

Bq=diag(Bq)

K=svd(G%*%Bq)$u

lam=diag(svd(G%*%Bq)$d)

L=svd(G%*%Bq)$v

if(l==1){
  
f=sqrt(n)*K[,1:p]%*%diag(1,p)

A=L%*%lam[,1:p]%*%diag(1,p)/sqrt(n)
  
}else{

f=f*(1-eta)+eta*sqrt(n)*K[,1:p]%*%diag(1,p)

A=A*(1-eta)+eta*L%*%lam[,1:p]%*%diag(1,p)/sqrt(n)

LLN=sum((G%*%Bq-f%*%t(A))^2)

print(LLN)

}

}

#p.49 5.2　分割表基準

Dat=matrix(c(28,10,12,19,23,8),ncol=2)

Dr=diag(apply(Dat,1,sum))

Dc=diag(apply(Dat,2,sum))

H=Dat

h=sum(H)

Hmat=diag(c(1/sqrt(diag(Dr))))%*%(H-Dr%*%array(1,dim=c(nrow(Dr),nrow(Dc)))%*%Dc/h)%*%diag(c(1/sqrt(diag(Dc))))

R_childer=svd(Hmat)$u

pthi=diag(svd(Hmat)$d) #通常はこの値によって列数を決める

C_childer=svd(Hmat)$v

pi=t(t(diag(Dr)))%*%t(diag(Dc))/h

#非対称解１

R=sqrt(h)*diag(c(1/sqrt(diag(Dr))))%*%R_childer%*%pthi

C=sqrt(h)*diag(c(1/sqrt(diag(Dc))))%*%C_childer

H_hat=Dr%*%(array(1,dim=c(nrow(Dr),nrow(Dc)))+R%*%t(C))%*%Dc/h


#非対称解２

R=sqrt(h)*diag(c(1/sqrt(diag(Dr))))%*%R_childer

C=sqrt(h)*diag(c(1/sqrt(diag(Dc))))%*%C_childer%*%pthi

H_hat=Dr%*%(array(1,dim=c(nrow(Dr),nrow(Dc)))+R%*%t(C))%*%Dc/h


#対称解

R=sqrt(h)*diag(c(1/sqrt(diag(Dr))))%*%R_childer%*%diag(sqrt(diag(pthi)))

C=sqrt(h)*diag(c(1/sqrt(diag(Dc))))%*%C_childer%*%diag(sqrt(diag(pthi)))

H_hat=Dr%*%(array(1,dim=c(nrow(Dr),nrow(Dc)))+R%*%t(C))%*%Dc/h

#p.52 5.3 分割表基準に基づく多重対応分析

G=matrix(c(0,0,1,1,0,0,1,1,0,0,0,0,0,0,0,0,1,1,1,0,1,0,1,0,0,1,0,1,0,1),ncol=5)

n=nrow(G);k=ncol(G)

m=rep(1,n)%*%G%*%rep(1,k)/n

J=diag(rep(1,n))-array(1,dim=c(n,n))/n

D=diag(c(rep(1,n)%*%G))

mat=J%*%G%*%diag(1/sqrt(diag(D)))

N=svd(mat)$u

pthi=diag(svd(mat)$d)

M=svd(mat)$v

t=diag(rep(1,k))

f=sqrt(n)*N%*%t

w=sqrt(n)*diag(1/sqrt(diag(D)))%*%M%*%pthi%*%t

G_hat=(array(1,dim=c(n,k))+f%*%t(w))%*%D/n