LoginSignup
0
0

More than 3 years have passed since last update.

@kozakai-ryouta(亮太 小酒井)

質的データの数量化―双対尺度法とその応用 (統計ライブラリー) 　朝倉書店

Last updated at Posted at 2019-03-06

#相対尺度法

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)
0
0
4

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
0
0