Go to Qiita Advent Calendar 2024 Top
LoginSignup
2
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 1 year has passed since last update.

@kozakai-ryouta(亮太 小酒井)

非対称多次元尺度構成法ー行動科学における多変量データ解析　現代数学社　千野直仁先生　

Last updated at Posted at 2021-08-23

#p.28 INDSCAL

I=5;J=K=10

bijk=data.frame(i=0,j=0,k=0,b=0)

for(i in 1:I){
for(j in 1:J){
for(k in 1:K){  
 
  bijk_sub=data.frame(i=i,j=j,k=k,b=i+j+k)
  
  bijk=rbind(bijk,bijk_sub)
  
}  
}
}

bijk=tail(bijk,nrow(bijk)-1)

X=matrix(sample(1:9,J*3,replace=T),ncol=3)

X=matrix(as.matrix(X),ncol=ncol(X),nrow=nrow(X))

w=t(array(0.1,dim=c(I,ncol(X))))

r=ncol(X)

lam=0.01

ite=1000

eta=0.01

for(l in 1:ite){

X2=array(0,dim=c(1,ncol(X)))

for(n in 1:nrow(bijk)){

X2_sub=X[bijk$j[n],]*X[bijk$k[n],]

X2=rbind(X2,X2_sub)
  
}

X2=tail(X2,nrow(X2)-1)

Xw=array(0,dim=c(1,ncol(X)))

for(s in 1:nrow(bijk)){
  
Xw=rbind(Xw,w[,bijk$i[s]]*X[bijk$j[s],]) 
  
}

Xw=tail(Xw,nrow(Xw)-1)

dist=c()

for(s in 1:K){
  
X[s,]=(1-eta)*X[s,]+eta*solve(t(Xw[bijk$k==s,])%*%Xw[bijk$k==s,]+diag(rep(lam,ncol(Xw))))%*%t(Xw[bijk$k==s,])%*%bijk$b[bijk$k==s]  
  
dist=c(dist,sum((bijk$b[bijk$k==s]-Xw[bijk$k==s,]%*%X[s,])^2))

}

for(i in 1:I){

w[,i]=(1-eta)*w[,i]+eta*solve(t(X2[bijk$i==i,])%*%X2[bijk$i==i,]+diag(rep(lam,ncol(X2))))%*%t(X2[bijk$i==i,])%*%bijk$b[bijk$i==i]

dist=c(dist,sum((X2[bijk$i==i,]%*%w[,i]-bijk$b[bijk$i==i])^2))

}

print(sum(dist))

}




#Weeks-Bentler p.36
#p.164 Gauss法

m=5;n=m;r=2
S=array(0,dim=c(m,m))
for(i in 1:m){
for(j in 1:n){
S[i,j]=rpois(1,i)+rpois(1,j)
}}
diag(S)=0


f=function(x){
X=matrix(x[c(1:(m*r))],ncol=r)
vec=x[-c(1:(m*r))]
a=vec[1];b=vec[2];c=vec[-c(1:2)]
d=array(0,dim=c(m,m))
for(j in 1:m){
d[j,]=sqrt(apply((t(X)-c(X[j,]))^2,2,sum))*b+c[j]-c+a
}
return(c(d))
}

Y=c(S)
h=0.001
lam=1
ite=10^5
sita=sample(1:9,m*r+m+2,replace=T)
eta=10^(-3)

for(l in 1:ite){
sita_pre=sita
Z=array(0,dim=c(m*n,length(sita)))
vec=sita
for(j in 1:ncol(Z)){
vec_sub=vec;vec_sub[j]=vec_sub[j]+h
Z_sub=(f(vec_sub)-f(vec))/h
Z[,j]=Z_sub
}
beta=solve(t(Z)%*%Z+diag(lam,ncol(Z)))%*%t(Z)%*%(Y-f(sita))
sita=c(sita+eta*beta)
print(sum((Y-f(sita))^2))
}




#Saito(1991) p.40
#p.164 Gauss法

m=5;n=m;r=2
S=array(0,dim=c(m,m))
for(i in 1:m){
for(j in 1:n){
S[i,j]=rpois(1,i)+rpois(1,j)
}}
diag(S)=0


f=function(x){
X=matrix(x[c(1:(m*r))],ncol=r)
vec=x[-c(1:(m*r))]
a=vec[1];vec=vec[-1]
b=vec[1:m];c=vec[-c(1:m)]
d=array(0,dim=c(m,m))
for(j in 1:m){
d[j,]=sqrt(apply((t(X)-c(X[j,]))^2,2,sum))+b[j]+c+a
}
return(c(d))
}

Y=c(S)
h=0.001
lam=1
ite=10^5
sita=sample(1:9,m*r+2*m+1,replace=T)
eta=10^(-3)

for(l in 1:ite){
sita_pre=sita
Z=array(0,dim=c(m*n,length(sita)))
vec=sita
for(j in 1:ncol(Z)){
vec_sub=vec;vec_sub[j]=vec_sub[j]+h
Z_sub=(f(vec_sub)-f(vec))/h
Z[,j]=Z_sub
}
beta=solve(t(Z)%*%Z+diag(lam,ncol(Z)))%*%t(Z)%*%(Y-f(sita))
sita=c(sita+eta*beta)
print(sum((Y-f(sita))^2))
}




#DEDICOM p.54
#p.164 Gauss法

m=5;n=m;r=2
S=array(0,dim=c(m,m))
for(i in 1:m){
for(j in 1:n){
S[i,j]=rpois(1,i)+rpois(1,j)
}}

f=function(x){
A=matrix(x[c(1:(r*r))],ncol=r)
y=matrix(x[-c(1:(r*r))],ncol=r)
d=y%*%A%*%t(y)
return(c(d))
}

Y=c(S)
h=0.001
lam=1
ite=10^5
sita=sample(1:9,(n+r)*r,replace=T)
eta=10^(-3)

for(l in 1:ite){
sita_pre=sita
Z=array(0,dim=c(m*n,length(sita)))
vec=sita
for(j in 1:ncol(Z)){
vec_sub=vec;vec_sub[j]=vec_sub[j]+h
Z_sub=(f(vec_sub)-f(vec))/h
Z[,j]=Z_sub
}
beta=solve(t(Z)%*%Z+diag(lam,ncol(Z)))%*%t(Z)%*%(Y-f(sita))
sita=c(sita+eta*beta)
print(sum((Y-f(sita))^2))
}



#拡張DEDICOM(Asymscal) p.54
#p.164 Gauss法

m=5;n=m;r=2
S=array(0,dim=c(m,m))
for(i in 1:m){
for(j in 1:n){
S[i,j]=rpois(1,i)+rpois(1,j)
}}

f=function(x){
a=x[1];b=x[2];x=x[-c(1:2)]
A=matrix(x[c(1:(r*r))],ncol=r)
y=matrix(x[-c(1:(r*r))],ncol=r)
d=a*(y%*%A%*%t(y))+b
return(c(d))
}

Y=c(S)
h=0.001
lam=1
ite=10^5
sita=sample(1:9,(n+r)*r+2,replace=T)
eta=10^(-3)

for(l in 1:ite){
sita_pre=sita
Z=array(0,dim=c(m*n,length(sita)))
vec=sita
for(j in 1:ncol(Z)){
vec_sub=vec;vec_sub[j]=vec_sub[j]+h
Z_sub=(f(vec_sub)-f(vec))/h
Z[,j]=Z_sub
}
beta=solve(t(Z)%*%Z+diag(lam,ncol(Z)))%*%t(Z)%*%(Y-f(sita))
sita=c(sita+eta*beta)
print(sum((Y-f(sita))^2))
}





#GIPSCAL p.110

mat=matrix(c(202,176,228,262,440,462,369,504,561),ncol=3)/10

#数量化4類で初期値設定

H=(t(mat)+mat)/2

alpha=-max(H)

H=H+alpha;diag(H)=0

for(j in 1:ncol(H)){

diag(H)[j]=-sum(H[j,])  

}

N=nrow(mat)

r=2

X=eigen(H)$vectors[,1:r]

a=1;b=1;c=1

B=diag(rep(1,N));D=diag(rep(1,N))

ite=100000

eta=0.001

for(l in 1:ite){

Y=t(X)%*%B%*%X;Z=t(X)%*%D%*%X

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

L=array(0,dim=c(r,r))

b_vec=c();l_vec=c()

count=1

for(j in 1:(r-1)){
  
b_vec=c(b_vec,b*count)

l_vec=c(l_vec,count)

count=count*(-1)
  
}

for(j in 2:r){
  
A[(j-1),j:r]=b_vec  

L[(j-1),j:r]=l_vec
  
}

A=A-t(A);diag(A)=a

L=L-t(L)

v=(mat-c)%*%X

u=t(mat-c)%*%X

s_bar=sum(mat)/prod(dim(mat))

p_bar=sum(X%*%t(X))/prod(dim(mat))

o_bar=sum(X%*%L%*%t(X))/prod(dim(mat))

sig_so=sum((mat-s_bar)*(X%*%L%*%t(X)-o_bar))/prod(dim(mat))

sig_po=sum((X%*%t(X)-p_bar)*(X%*%L%*%t(X)-o_bar))/prod(dim(mat))

sig_p=sum((X%*%t(X)-p_bar)^2)/prod(dim(mat))

sig_o=sum((X%*%L%*%t(X)-o_bar)^2)/prod(dim(mat))

sig_sp=sum((mat-s_bar)*(X%*%t(X)-p_bar))/prod(dim(mat))

a=(sig_sp*sig_o-sig_so*sig_po)/(sig_p*sig_o-sig_po^2)

b=(sig_so*sig_p-sig_sp*sig_po)/(sig_p*sig_o-sig_po^2)

c=s_bar-a*p_bar-b*o_bar

Cht=t(A)%*%Y%*%A+A%*%Z%*%t(A)

X=(1-eta)*X+eta*t(solve(Cht)%*%(t(A)%*%t(u)+A%*%t(v)))

Q=sum((mat-X%*%t(X)*a-b*X%*%L%*%t(X)-c)^2)

print(Q)


}
2
0
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
2
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?