
#相対尺度法

mat=t(matrix(c(1,3,6,3,5,2,6,3,0),ncol=3))

dim=nrow(mat)

f=apply(mat,2,sum)

D=diag(f);tf=sum(f)

D_n=diag(apply(mat,1,sum))

C0=sqrt(solve(D_n))%*%t(mat)%*%solve(D_n)%*%mat%*%sqrt(solve(D_n))

#iteration

C1=C0-sqrt(D)%*%array(1,dim=c(dim,1))%*%array(1,dim=c(1,dim))%*%sqrt(D)/tf

b0=t(t(c(1,0,-1)))

times=10

b=b0

b_vec=array(0,dim=c(times,dim))

for(j in 1:times){
  
b_dash=C1%*%b

b_dash=b_dash/max(abs(b_dash))

b=b_dash

print(b)

b_vec[j,]=b
  
}

w=c(sqrt(tf/(t(b)%*%b)))*c(b)

x=sqrt(solve(D))%*%w

#相関比を最大にするように計算する(SSt/SSb=eta2)

eta2=t(x)%*%t(mat)%*%solve(D_n)%*%mat%*%x/(t(x)%*%D%*%x)

y=solve(D_n)%*%mat%*%x/c(sqrt(eta2))

SSt=t(x)%*%D%*%x;SSb=t(x)%*%t(mat)%*%solve(D_n)%*%mat%*%x

SSt/SSb

#lagrange function

Q=t(x)%*%t(mat)%*%solve(D_n)%*%mat%*%x-eta2*(t(x)%*%D%*%x-tf)  

b_hat=b

delta1=eta2/sum(diag(C1))
  
SSw=t(x)%*%D%*%x-t(x)%*%t(mat)%*%solve(D_n)%*%mat%*%x

g=t(t(diag(D_n)))

#SSbになるはず

 t(y)%*%(mat%*%solve(D)%*%t(mat)-g%*%t(g)/c(tf))%*%y
 
#SStになるはず
 
 t(y)%*%(D_n-g%*%t(g)/c(tf))%*%y
 
#偏差率相関
 
 p=t(y)%*%mat%*%x/(sqrt(t(x)%*%D%*%x)*sqrt(t(y)%*%D_n%*%y))
 
 t(y)%*%mat%*%x/t(x)%*%D%*%x
 
 #lagrange function p,52 (2.54)
 
 Q2=t(y)%*%mat%*%x-c(p)*(t(x)%*%D%*%x-c(tf))-c(p)*t(y)%*%D_n%*%y
 
 #yとxの関係
 
 floor(c(c(p)*y)*1000)/1000==floor(c(solve(D_n)%*%mat%*%x)*1000)/1000
 
 solve(D_n)%*%mat/(c(p))
 
 
 u=sqrt(D_n)%*%y;v=sqrt(D)%*%x
 

 
B=sqrt(solve(D_n))%*%mat%*%sqrt(solve(D))

eta_mat=(eigen(t(B)%*%B)$values)

W=eigen(t(B)%*%B)$vectors

U=eigen((B)%*%t(B))$vectors

X=sqrt(solve(D))%*%W;Y=sqrt(solve(D_n))%*%U

#U

B%*%W%*%solve(diag(sqrt(eta_mat)))

#W

t(B)%*%U%*%solve(diag(sqrt(eta_mat)))

#singular(3.21)

t(U)%*%B%*%W

t(Y)%*%mat%*%X

sqrt(diag(eta_mat))

iterations=50

x_hat=x;y_hat=y

x_hat_mat=array(0,dim=c(iterations,length(x_hat)))

y_hat_mat=array(0,dim=c(iterations,length(y_hat)))


for(j in 1:iterations){
  
y_hat=solve(D_n)%*%mat%*%x_hat  
  
x_hat=solve(D)%*%t(mat)%*%y_hat

x_hat_mat[j,]=x_hat

y_hat_mat[j,]=y_hat

}


f=mat

iterations=50

x0=x_hat;y0=y_hat

for(j in 1:iterations){
  
f=f-f%*%(x_hat%*%t(y_hat))%*%f/c(t(y_hat)%*%f%*%x_hat)

y0=solve(D_n)%*%f%*%x0

x0=solve(D)%*%t(f)%*%y0

tff=t(y0)%*%D_n%*%y0
  
}

eta=eta2

w_vec=t(t(w))

mat_C=C1

times=30

b0=c(-1,0,1)

b=b0

mat_C=mat_C-c(eta)*(w_vec)%*%t(w_vec)/c(t(w_vec)%*%w_vec)

x_mat=array(0,dim=c(times,length(w_vec)))

y_mat=array(0,dim=c(times,length(w_vec)))

etas=c()

for(j in 1:times){
  
b=mat_C%*%b

b=b/max(abs(b))

w_vec=c(sqrt(tf/(t(b)%*%b)))*b

x_vec=solve(D)%*%w_vec

eta_vec=t(x_vec)%*%t(mat)%*%solve(D_n)%*%mat%*%x_vec/(t(x_vec)%*%D%*%x_vec)

etas=c(etas,eta_vec)

y_vec=solve(D_n)%*%mat%*%x_vec/c(eta_vec)

x_mat[j,]=x_vec;y_mat[j,]=y_vec

print(b)
  
}

w_vec=c(sqrt(tf/(t(b)%*%b)))*b

x_vec=solve(D)%*%w_vec

eta_vec=t(x_vec)%*%t(mat)%*%solve(D_n)%*%mat%*%x_vec/(t(x_vec)%*%D%*%x_vec)

y_vec=solve(D_n)%*%mat%*%x_vec/c(eta_vec)

t(w_vec)%*%w_vec

eta_mat=matrix(rep(etas,length(x_vec)),ncol=length(x_vec))

x_mat_star=x_mat*eta_mat

y_mat_star=y_mat*eta_mat

n=dim(mat)[1];m=dim(mat)[2]

kai2=-(tf-1-(n+m-1)/2)*log(1-etas)

chisq=c()

chisq_sum=c()

for(j in 1:times){

chisq=c(chisq,qchisq(1-0.05,df=(n+m-1-2*j)))

}

kai_sum=c()

for(j in 1:times){
  
kai_sum=c(kai_sum,sum(kai2[j:times])) 
  
chisq_sum=c(chisq_sum,qchisq(1-0.05,df=(n-j)*(m-j)))  
  
}

g=apply(mat,1,sum)


#3.3 定比例の定理

k=5;N=10

mat=array(0,dim=c(N,1))

vecs=c()

for(j in 1:k){
  
sub=array(0,dim=c(N,sample(c(2:5),1)))  

vecs=c(vecs,ncol(sub))
  
for(i in 1:N){
  
val=sample(c(1:10),ncol(sub))

p=val/sum(val)

subs=data.frame(no=c(1:ncol(sub)),sig=ifelse(p==max(p),1,0))

sub[i,]=c(subs$sig)
  
}

colnames(sub)=paste0(rep(j,ncol(sub)))

mat=cbind(mat,sub)
  
}

mat=mat[,2:ncol(mat)]

f=as.matrix(mat)

ff=t(f)%*%f

X=t(sample(c(1:10),k))

no=colnames(mat)

vec=c()

for(j in 1:k){
  
vec=c(vec,rep(X[j],length(no[no==paste0(j)])))  
  
}

X=t(vec);



X_data=data.frame(no=no,X=c(X))

D=diag(diag(ff),nrow(ff))

eta2=(X)%*%t(f)%*%f%*%t((X))/(k*(X)%*%D%*%t((X)))

xffx=c();xdx=c();r=c()

for(j in 1:k){
  
xj=t(X_data$X[X_data$no==j])  
  
fj=f[,colnames(f) %in% paste0(j)]  

Dj=ff[rownames(ff) %in% paste0(j),colnames(ff) %in% paste0(j)]

xffx=c(xffx,xj%*%t(fj)%*%fj%*%t(xj))

xdx=c(xdx,xj%*%Dj%*%t(xj))

et=xj%*%t(fj)%*%fj%*%t(xj)/(k*xj%*%Dj%*%t(xj))

print(xj%*%t(fj)%*%fj%*%t(xj)/(k*xj%*%Dj%*%t(xj)))

r=c(r,(xj%*%t(fj)%*%fj%*%t(xj))^2/(xj%*%Dj%*%t(xj)))

}

sum(xffx)/sum(xdx)

r=(k*et*xffx)/sum(xffx)


#悪夢の頻度と睡眠薬に対する賛否

mat=t(matrix(c(15,8,3,2,0,5,17,4,0,2,6,13,4,3,2,0,7,7,5,9,1,2,6,3,16),ncol=5))

rownames(mat)=c("強く反対","反対","中立","賛成","強く賛成")

colnames(mat)=c("皆無","マレ","ときどき","しばしば","常時")

w=c(-2,-1,0,1,2)

y_all=(apply(mat,1,sum))

y_star=mat%*%w/y_all

x_all=apply(mat,2,sum)

x_star=t(mat)%*%w/x_all

D_n=diag(apply(mat,1,sum))

D=diag(apply(mat,2,sum))

p=t(y_star)%*%mat%*%x_star/(sqrt(t(x_star)%*%D%*%x_star)*sqrt(t(y_star)%*%D_n%*%y_star))

plot(w,x_star,type="o",col=2,ylim=c(min(x_star,y_star),max(x_star,y_star)),xlab="w",ylab="values")

par(new=T)

plot(w,y_star,type="o",col=3,ylim=c(min(x_star,y_star),max(x_star,y_star)),xlab="w",ylab="values")


#一対比較データの数量化

library(dplyr)

mat=t(matrix(c(1,2,1,2,2,1,0,2,2,2,1,1,2,2,2,2,2,1,1,2,0,2,1,1,1,2,2,2,0,1,2,1,1,1,2,1,1,2,1,1,1,1,1,2,2,2,2,2),ncol=8))

rownames(mat)=c(1,2,3,4,5,6,7,8)

colnames(mat)=c("(A,B)","(A,C)","(A,D)","(B,C)","(B,D)","(C,D)")

f=mat;f[f==2]=-1

alphabet=c("A","B","C","D")

n=length(alphabet)

colname=colnames(mat)

A=array(0,dim=c(n*(n-1)/2,n))

ns=0

for(j in 1:(n-1)){
  
ns_sub=ns+1
    
ns=ns+n-j

A[ns_sub:ns,j]=rep(1,(n-j))
    
A[ns_sub:ns,(j+1):ncol(A)]=diag(-1,(n-j))  
  
}

E=f%*%A

apply(E,1,sum)

colnames(E)=alphabet

#刺激数と被験者数

n=length(alphabet);N=nrow(mat)

D_n=diag(n*(n-1),N)

D=diag(N*(n-1),n)

#総反応数

f_t=N*n*(n-1)

H_n=t(E)%*%E/(N*n*(n-1)^2)

X=eigen(H_n)$vectors

eta2=eigen(H_n)$values

SSb=c();SSt=c()

for(j in 1:ncol(X)){
  
SSb=c(SSb,X[,j]%*%t(E)%*%solve(D_n)%*%E%*%X[,j])


SSt=c(SSt,X[,j]%*%D%*%X[,j])    
  
}


#順位データの数量化

mat=t(matrix(c(1,3,2,1,2,3,2,3,1,2,3,1),ncol=4))

E=(max(mat)+1-2*mat)

N=nrow(mat);n=ncol(mat)

D_n=diag(n*(n-1),N)

D=diag(N*(n-1),n)

#総反応数

f_t=N*n*(n-1)

H_n=t(E)%*%E/(N*n*(n-1)^2)

X=eigen(H_n)$vectors

eta2=eigen(H_n)$values

SSb=c();SSt=c()

for(j in 1:ncol(X)){
  
SSb=c(SSb,X[,j]%*%t(E)%*%solve(D_n)%*%E%*%X[,j])


SSt=c(SSt,X[,j]%*%D%*%X[,j])    
  
}

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

for(j in 1:n){

y_mat[,j]=solve(D_n)%*%E%*%c(X[,j])/sqrt(eta2[j])

}


#継時カテゴリーデータの数量化

mat=t(matrix(c(1,3,2,2,3,2,3,2,3,1,3,2,2,2,2),ncol=5))

colnames(mat)=c("A","B","C")

rownames(mat)=c(1:5)

f=array(0,dim=c(nrow(mat),1))

values=order(unique(c(mat)))


for(j in 1:ncol(mat)){
  
mat_sub=array(0,dim=c(nrow(mat),ncol(mat)))  
  
for(i in 1:nrow(mat)){  

vec=array(0,dim=c(1,ncol(mat)))


if(mat[i,j]!=ncol(mat) & mat[i,j]!=1){  
  
vec[1,(mat[i,j]+1):ncol(mat)]=1 
  
vec[1,1:(mat[i,j]-1)]=-1
  
}else{

if(mat[i,j]==1){
  
 vec[1,2:ncol(mat)]=1 
  
}
 
if(mat[i,j]==ncol(mat)){
  
vec[1,1:(ncol(mat)-1)]=-1  
  
}  
  
}

mat_sub[i,]=vec

}

f=cbind(f,mat_sub)

}

t_f=t(f)

f=t(matrix(t_f[t_f!=0],ncol=nrow(mat)))

A=diag(1,length(values)-1)

for(j in 1:(length(values)-1)){
  
A=rbind(A,diag(1,length(values)-1))  

}

A_sub=array(0,dim=c(2*length(values),length(values)))

for(j in 1:length(values)){
  
A_sub[(2*(j-1)+1):(2*j),j]=-1  
  
}

A=cbind(A,A_sub)

E=f%*%A

apply(E,1,sum)