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

```{r}

#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))

}

```

```{r}


#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))
}

```

```{r}


#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))
}

```


```{r}


#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))
}

```

```{r}

#拡張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))
}



```

```{r}

#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)


}

```