多次元尺度構成法　朝倉書店　斎藤堯幸先生　

```{r}




#Cooper(1972) p.44

mat=matrix(0,nrow=10,ncol=10)

mat[1,]=c(0,8.89,10.16,15.84,16.84,33.01,32.69,43.47,37.78,39.31)
mat[2,]=c(0,0,9.64,15.37,20.45,30.75,34.09,41.61,39.26,41.47)
mat[3,]=c(0,0,0,6.12,11.2,22.98,24.56,33.61,29.75,31.86)
mat[4,]=c(0,0,0,0,8.48,17.20,18.75,27.64,23.91,26.23)
mat[5,]=c(0,0,0,0,0,20.77,16.25,29.52,21.17,22.49)
mat[6,]=c(0,0,0,0,0,0,13.01,10.86,15.61,19.06)
mat[7,]=c(0,0,0,0,0,0,0,15.99,5.19,7.69)
mat[8,]=c(rep(0,8),14.84,17.68)
mat[9,]=c(rep(0,9),3.45)

S=mat+t(mat)
n=nrow(S)
r=2

#数量化4類
e=-S
H=e+t(e)
alpha=-max(H)
H=H+alpha;diag(H)=0
for(j in 1:ncol(H)){
diag(H)[j]=-sum(H[j,])  
}
ys=eigen(H)$vector[,1:r]
alpha=1

ite=10000
eta=0.01

for(l in 1:ite){
  
dys=array(0,dim=dim(ys)) 
pre=array(0,dim=dim(S))
for(i in 1:n){  
value=-(t(ys)-ys[i,])
dys[i,]=-2*apply(t(value)*(S[i,]-sqrt(apply(value^2,2,sum))-alpha)/ifelse(sqrt(apply(value^2,2,sum))==0,1,sqrt(apply(value^2,2,sum))),2,sum)
pre[i,]=sqrt(apply(value^2,2,sum))
}
 
ys=ys-eta*dys
alpha=sum(S-pre)/(n^2)

cost=sum((S-pre-alpha)^2)
print(cost)
 
}


```

```{r}




#Messick & Abelson p.45

#Saitoの方法

mat=matrix(0,nrow=10,ncol=10)

mat[1,]=c(0,8.89,10.16,15.84,16.84,33.01,32.69,43.47,37.78,39.31)
mat[2,]=c(0,0,9.64,15.37,20.45,30.75,34.09,41.61,39.26,41.47)
mat[3,]=c(0,0,0,6.12,11.2,22.98,24.56,33.61,29.75,31.86)
mat[4,]=c(0,0,0,0,8.48,17.20,18.75,27.64,23.91,26.23)
mat[5,]=c(0,0,0,0,0,20.77,16.25,29.52,21.17,22.49)
mat[6,]=c(0,0,0,0,0,0,13.01,10.86,15.61,19.06)
mat[7,]=c(0,0,0,0,0,0,0,15.99,5.19,7.69)
mat[8,]=c(rep(0,8),14.84,17.68)
mat[9,]=c(rep(0,9),3.45)

mat=mat+t(mat)

n=nrow(mat)

r=3

S=mat

S2=S^2

Im=t(array(1,dim=c(n,1)))

H=t(Im)%*%(Im)-diag(1/n,n)

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

alpha=sample(1:50,1)


ite=10^3

eta=0.01

for(l in 1:ite){

S_a=S+alpha*J

S2_a=S2+2*alpha*S+J*alpha^2

B=-(1/2)*t(H)%*%S2_a%*%H

eigen_values=eigen(B)$values

eigen_vectors=eigen(B)$vectors

R_a=sum(eigen_values[eigen_values>0])

R1_a=sum(eigen_values[1:r])

X=eigen_vectors

term2=c()

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

term2=c(term2,t(X[,j])%*%S%*%X[,j])

}

dlam=-term2+alpha

P=sum(eigen_values[1:r]^2)/sum(eigen_values[eigen_values>0]^2)

lam=eigen_values[eigen_values>0]

dP=2*(sum(lam[1:r]*dlam[1:r])-P*sum(lam*dlam[1:length(lam)]))/sum(lam^2)

alpha=alpha*(1-eta)+eta*dP

print(P)

}

```

```{r}



#Messick & Abelson p.45

mat=matrix(0,nrow=10,ncol=10)

mat[1,]=c(0,8.89,10.16,15.84,16.84,33.01,32.69,43.47,37.78,39.31)
mat[2,]=c(0,0,9.64,15.37,20.45,30.75,34.09,41.61,39.26,41.47)
mat[3,]=c(0,0,0,6.12,11.2,22.98,24.56,33.61,29.75,31.86)
mat[4,]=c(0,0,0,0,8.48,17.20,18.75,27.64,23.91,26.23)
mat[5,]=c(0,0,0,0,0,20.77,16.25,29.52,21.17,22.49)
mat[6,]=c(0,0,0,0,0,0,13.01,10.86,15.61,19.06)
mat[7,]=c(0,0,0,0,0,0,0,15.99,5.19,7.69)
mat[8,]=c(rep(0,8),14.84,17.68)
mat[9,]=c(rep(0,9),3.45)

mat=mat+t(mat)

n=nrow(mat)

r=2

S=mat

S2=S^2

Im=t(array(1,dim=c(n,1)))

H=t(Im)%*%(Im)-diag(1/n,n)

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

alpha=sample(1:50,1)


ite=10

for(l in 1:ite){

S_a=S+alpha*J

S2_a=S2+2*alpha*S+J*alpha^2

B=-(1/2)*t(H)%*%S2_a%*%H

eigen_values=eigen(B)$values

eigen_vectors=eigen(B)$vectors

R_a=sum(eigen_values[eigen_values>0])

R1_a=sum(eigen_values[1:r])

X=eigen_vectors

term1=c()

term2=c()

for(j in 1:r){

term1=c(term1,t(X[,j])%*%B%*%X[,j])  

term2=c(term2,t(X[,j])%*%S%*%X[,j])

}

term1=sum(term1);term2=sum(term2)

alpha=Re(polyroot(c(term1-R1_a,term2,r/2)))

alpha=alpha[term1-alpha*term2+r*alpha^2/2==min(term1-alpha*term2+r*alpha^2/2)]

P=sum(eigen_values[1:r]^2)/sum(eigen_values[eigen_values>0]^2)

print(P)

}

```




```{r}

#多次元尺度構成法　朝倉書店　斎藤堯幸先生　

#KL数量化 p.64

mat=matrix(0,nrow=10,ncol=10)

mat[1,]=c(0,8.89,10.16,15.84,16.84,33.01,32.69,43.47,37.78,39.31)
mat[2,]=c(0,0,9.64,15.37,20.45,30.75,34.09,41.61,39.26,41.47)
mat[3,]=c(0,0,0,6.12,11.2,22.98,24.56,33.61,29.75,31.86)
mat[4,]=c(0,0,0,0,8.48,17.20,18.75,27.64,23.91,26.23)
mat[5,]=c(0,0,0,0,0,20.77,16.25,29.52,21.17,22.49)
mat[6,]=c(0,0,0,0,0,0,13.01,10.86,15.61,19.06)
mat[7,]=c(0,0,0,0,0,0,0,15.99,5.19,7.69)
mat[8,]=c(rep(0,8),14.84,17.68)
mat[9,]=c(rep(0,9),3.45)
s=mat+t(mat)
n=nrow(s)

#数量化4類
e=-s
H=e+t(e)
alpha=-max(H)
H=H+alpha;diag(H)=0
for(j in 1:ncol(H)){
diag(H)[j]=-sum(H[j,])  
}
x=eigen(H)$vector[,1]
L=1


ite=10^4
eta=10^(-4)

for(l in 1:ite){
D=(t(matrix(rep(x,n),ncol=n))-c(x))
dx=apply((s-L)*D-D^3,1,sum)
x=x-eta*dx
x=x-mean(x)
L=sum(s-D^2)/(n^2)
W2=sum((s-L-D^2)^2)
print(W2)
}


```


```{r}

#schonemann method --見直し中--

library(dplyr)

m=100;n=50;r=10

X=array(0,dim=c(m,r));Y=array(0,dim=c(n,r))

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

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


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

for(j in 1:n){
for(i in 1:m){  
  
S[j,i]=sqrt(t(X[i,]-Y[j,])%*%t(t(X[i,]-Y[j,])))  
  
}
}


S_hat=S^2

S_hat=S_hat-apply(S^2,1,mean)

S_hat=t(t(S_hat)-apply(S^2,2,mean))

S_hat=S_hat+mean(S^2)

B_star=-S_hat/2

A=B_star

lambda=diag(svd(A)$d);U=svd(A)$u;V=svd(A)$v

B=t(A)%*%A

eigen_values=eigen(B)$values

eigen_values=eigen_values[eigen_values>10^(-5)]

r=length(eigen_values)

eigen_vectors=eigen(B)$vectors[,1:r]

vec=sqrt(apply(eigen_vectors^2,2,sum))

Q=eigen_vectors

for(j in 1:r){
  
  Q[,j]=Q[,j]/vec[j]
  
}

lambda=diag(eigen_values)

P=A%*%Q%*%diag(1/eigen_values)

A_hat=P%*%diag(eigen_values)%*%t(Q)

Y0=P;X0=diag(eigen_values)%*%t(Q)

X0=t(X0)

if(ncol(S)-1<=r*(r+3)/2){

K1=array(0,dim=c(n-1,r*(r+1)/2))

K2=array(0,dim=c(n-1,r))

f=c()

for(j in 1:(n-1)){
  
Y1=Y0[j,];Y2=Y0[n,]  

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

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

for(l in 1:r){
for(k in 1:r){  
if(l>=k){
 
mat1[l,k]=Y1[l]*Y1[k]

mat2[l,k]=Y2[l]*Y2[k]
  
}  
}
}
  
mat1[upper.tri(mat1)]=NA;mat2[upper.tri(mat2)]=NA

K1_vec=as.numeric(mat1[is.na(mat1)!=T])-as.numeric(mat2[is.na(mat2)!=T])
  
K1[j,]=K1_vec

K2[j,]=as.numeric(Y1-Y2)

f=c(f,apply(S^2,1,mean)[j]-apply(S^2,1,mean)[n])

}

K=cbind(K1,K2)

matrix=t(K)%*%K

params=(svd(matrix)$v)%*%diag(1/svd(matrix)$d)%*%t(svd(matrix)$u)%*%t(K)%*%f

n_vec=params[1:(r*(r+1)/2)];z_vec=params[((r*(r+1)/2)+1):length(params)]

u_vec=z_vec/2

n_tt=c(n_vec[1])

n_vec_sub=n_vec[-1]

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

for(j in (r-1):1){
  
no=r-j  
  
M[no,(no+1):r]=n_vec_sub[c(1:j)]  
  
n_vec_sub=n_vec_sub[-c(1:j)]  
  
n_tt=c(n_tt,n_vec_sub[1])

}

M=M+t(M);diag(M)=n_tt/2

lambda2=eigen(M)$values;V=eigen(M)$vector

lambda2=lambda2[lambda2>0]

V=V[,1:length(lambda2)]

t=V%*%diag(sqrt(lambda2))

error=sum(abs(M-t%*%t(t)))/prod(dim(M))

inverse_t=diag(1/sqrt(lambda2))%*%t(V)

a_b=inverse_t%*%u_vec

X_hat=X0%*%t(inverse_t);X_hat=X_hat+matrix(rep(a_b,nrow(X_hat)),nrow=nrow(X_hat))

Y_hat=Y0%*%t

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

for(j in 1:n){
for(i in 1:m){  
  
d[j,i]=sqrt(t(X_hat[i,]-Y_hat[j,])%*%t(t(X_hat[i,]-Y_hat[j,])))  
  
}
}



}

```
```{r}

#Eckart-Young decomposition

library(dplyr)

m=100;n=20

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

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


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

for(j in 1:m){
for(i in 1:m){  
 
dist[i,j]=sum((X[i,]-X[j,])^2)   
  
}
}

S=dist

S_hat=S

for(j in 1:m){
for(i in 1:m){  
  
S_hat[i,j]=S[i,j]^2+mean(S^2)-apply(S^2,1,mean)[i]-apply(S^2,2,mean)[j]  
  
}
}

lambda=diag(svd(A)$d);U=svd(A)$u;V=svd(A)$v

B=t(A)%*%A

eigen_values=eigen(B)$values

eigen_values=eigen_values[eigen_values>10^(-5)]

r=length(eigen_values)

eigen_vectors=eigen(B)$vectors[,1:r]

vec=sqrt(apply(eigen_vectors^2,2,sum))

Q=eigen_vectors

for(j in 1:r){
  
  Q[,j]=Q[,j]/vec[j]
  
}

lambda=diag(eigen_values)

P=A%*%Q%*%diag(1/eigen_values)

A_hat=P%*%diag(eigen_values)%*%t(Q)


```

```{r}

#schonemann method p.121
#p.164 Gauss法

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

f=function(x){
X=matrix(x[c(1:(m*r))],ncol=r)
Y=matrix(x[-c(1:(m*r))],ncol=r)
d=array(0,dim=c(m,n))
for(i in 1:m){
d[i,]=apply((t(Y)-X[i,])^2,2,sum)
}
return(c(sqrt(d)))
}

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

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}

#p.131 B-T-L model (Bradley-Terry-Luce model)

mat=matrix(0,nrow=10,ncol=10)

mat[1,]=c(0,8.89,10.16,15.84,16.84,33.01,32.69,43.47,37.78,39.31)
mat[2,]=c(0,0,9.64,15.37,20.45,30.75,34.09,41.61,39.26,41.47)
mat[3,]=c(0,0,0,6.12,11.2,22.98,24.56,33.61,29.75,31.86)
mat[4,]=c(0,0,0,0,8.48,17.20,18.75,27.64,23.91,26.23)
mat[5,]=c(0,0,0,0,0,20.77,16.25,29.52,21.17,22.49)
mat[6,]=c(0,0,0,0,0,0,13.01,10.86,15.61,19.06)
mat[7,]=c(0,0,0,0,0,0,0,15.99,5.19,7.69)
mat[8,]=c(rep(0,8),14.84,17.68)
mat[9,]=c(rep(0,9),3.45)

d=mat+t(mat)
S=d[1:6,]
m=nrow(S)
n=ncol(S)
r=3

ite=10000
eta=0.0001

xs=array(0.1,dim=c(n,r))
ys=array(0.001,dim=c(m,r))
tau=rep(1,m)


for(l in 1:ite){
  
ys_pre=ys
xs_pre=xs
tau_pre=tau
  
dys=array(0,dim=dim(ys)) 
dtau=rep(0,length(tau))
pre=array(0,dim=dim(S))
for(i in 1:m){
dys[i,]=4*apply(t(t(xs)-ys[i,])*c(S[i,]-apply((t(xs)-ys[i,])^2,2,sum)-rep(tau[i],n)),2,sum)
dtau[i]=-4*sum(S[i,]-apply((t(xs)-ys[i,])^2,2,sum)-rep(tau[i],n))
pre[i,]=apply((t(xs)-ys[i,])^2,2,sum)+rep(tau[i],n)
}

dxs=array(0,dim=dim(xs))
for(j in 1:n){
dxs[j,]=4*apply(t(t(ys)-xs[j,])*c(S[,j]-apply((t(ys)-xs[j,])^2,2,sum)-tau),2,sum)
}

ys=ys-eta*dys
xs=xs-eta*dxs
tau=tau-eta*dtau
print(sum((S-pre)^2))

}


```



```{r}
#p.160 7.4 確率的アプローチ

n=5;m=6

I=n
J=m
K=m

S=list()
for(i in 1:I){
dmat=array(0,dim=c(J,K))  
for(j in 1:J){
for(k in 1:K){
dmat[j,k]=rpois(1,J+1)+rpois(1,K+2) 
}}
dmat=dmat+rpois(1,i)
diag(dmat)=0
S[[i]]=dmat
}



r=3

X=array(runif(m*r),dim=c(m,r))
W=array(1,dim=c(n,r))

ite=10^4
eta=10^(-3)

Svec=c()

for(l in 1:ite){

X_pre=X
W_pre=W

dij=list()
yijk=list()
Wit=array(0,dim=c(n,r))
Sval=0
for(i in 1:n){
d=array(0,dim=c(m,m))  
for(j in 1:m){  
d[j,]=sqrt(apply(((t(X)-X[j,])^2)*W[i,],2,sum)) 
}
dij[[i]]=d
yijk[[i]]=(log(S[[i]]/d)+log(t(S[[i]])/t(d)))/(d^2)
diag(yijk[[i]])=0
mt=S[[i]]/d;diag(mt)=1
Sval=Sval+sum(log(mt)^2)
}


Wit=array(0,dim=c(n,r))
Xjt=array(0,dim=c(m,r))
for(t in 1:r){
Xd=array(0,dim=c(m,m))
Xd2=array(0,dim=c(m,m))
for(j in 1:m){
Xd[j,]=(X[,t]-X[j,t])^2  
Xd2[j,]=(X[,t]-X[j,t])
}
for(i in 1:n){
mat=-Xd*(log(S[[i]])-log(dij[[i]]))/(dij[[i]]^2)
diag(mat)=0
Wit[i,t]=sum(mat)
mat=Xd2*(log(S[[i]])-log(dij[[i]]))/(dij[[i]]^2)
diag(mat)=0
Xjt=t(t(apply(mat,1,sum)))%*%t(c(W[i,]))
}}



W=W-eta*Wit
for(t in 1:r){
W[,t]=sqrt(n)*W[,t]/sqrt(sum(W[,t]^2))  
}
X=X-eta*Xjt

Svec=c(Svec,Sval)

#print(sum(abs(W-W_pre))+sum(abs(X-X_pre)))
print(Sval)
}

plot(Svec)
```