質的データの判別分析　数量化2類　菅民郎・藤越康祝先生　現代数学社

```{r}

#p.75 

#第1章　数量化2類に用いられる基本解析手法

#1.2 クラメール連関係数

#店舗評価

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

rownames(mat)=c("駅周辺","ビル街","商店街")

colnames(mat)=c("優良店","不良店")

#期待度数(立地特性と店舗評価の関連性がないとしてみなしたもの)

expect_freq=t(t(t(c(apply(mat,2,sum))))%*%t(c(apply(mat,1,sum)))/sum(mat))

#適合度検定

kai2=(mat-expect_freq)^2/expect_freq

#クラメールの連関係数

r=sqrt(sum(kai2)/(sum(mat)*(min(dim(mat))-1)))



```

```{r}

#1.4 数量化2類サンプルスコアの全体変動

library(dplyr)

data=data.frame(group=c(rep(1,5),rep(2,7)),score=c(0.299,1.562,1.957,0.178,0.859,-0.777,-0.799,-0.777,-0.799,-0.404,0.178,-1.48))

ave=data%>%group_by(group)%>%summarise(ave=mean(score))

sy_hat=sum((mean(data$score)-data$score)^2)

#群間変動、郡内変動

dat=data.frame(product=c(rep(1,4),rep(2,5),rep(3,6)),age=c(29,32,35,36,38,40,41,43,48,20,22,24,29,35,38))

ave=mean(dat$age)

aves=dat%>%group_by(product)%>%summarise(ave_age=mean(age))

#群内変動

val=c()

for(j in 1:nrow(aves)){
  
val=c(val,sum((dat$age[dat$product==j]-aves$ave_age[aves$product==j])^2))
  
}

#全体変動

st=sum((dat$age-mean(dat$age))^2)

#群間変動

sb=st-sum(val)

```

```{r}

#1.7 全体変動行列、群間変動行列、群内変動行列

#p.89

data=data.frame(group=c(rep(1,4),rep(2,5),rep(3,6)),x1=c(29,32,35,36,38,40,41,43,48,20,22,24,29,35,38),x2=c(32,29,36,35,38,35,37,45,40,22,20,38,35,29,24),x3=c(0,1,0,1,1,0,1,1,1,0,1,0,0,1,1))

aves=data%>%group_by(group)%>%summarise(x1=mean(x1),x2=mean(x2),x3=mean(x3))

#全体変動行列

mat=as.matrix(data[,colnames(data) %in% c("x1","x2","x3")])

t_mat=t(t(mat)-apply(mat,2,mean))

#郡内変動行列

mat=as.matrix(data[,colnames(data) %in% c("x1","x2","x3")])

for(j in 1:nrow(aves)){
  
ave=aves%>%filter(group==j)%>%select(-group)

ave=c(unlist(ave))

mat[c(1:nrow(data))[data$group==j],]=t(t(mat[c(1:nrow(data))[data$group==j],])-ave)
  
}

W_mat=t(mat)%*%mat

#群間変動行列

mat=as.matrix(data[,colnames(data) %in% c("x1","x2","x3")])

ave_vec=apply(mat,2,mean)

for(j in 1:nrow(aves)){
  
ave=aves%>%filter(group==j)%>%select(-group)

ave=c(unlist(ave))

mat[c(1:nrow(data))[data$group==j],]=t(matrix(rep(ave,length(c(1:nrow(data))[data$group==j])),ncol=length(c(1:nrow(data))[data$group==j])))
  
}

mat=t(t(mat)-c(ave_vec))

B_mat=t(mat)%*%mat



```



```{r}

#質的データの判別分析　数量化2類

#p.137~172 数量化2類

XX=matrix(c(rep(0,5),rep(1,3),0,1,1,rep(0,6),1,0,0,1,1,1,0,0,0,0,0,0,0,1,0,1,0,0,1,1,0,0,0,1,0,1,0,0,0,1,1,0,0,0,0,1,0),ncol=6)

#XXの3列目、6列目を抜く
X=matrix(c(rep(0,5),rep(1,3),0,1,1,rep(0,6),1,0,0,0,1,0,1,0,0,1,1,0,0,0,1,0,1,0,0),ncol=4)

X_childer=X;n=nrow(X)

for(j in 1:ncol(X)){
  
X_childer[,j]=X_childer[,j]-apply(X,2,mean)[j]

}
  
t=t(X_childer)%*%X_childer

X1=X[1:3,];num1=c(1:3)

X2=X[4:5,];num2=c(4:5)

X3=X[6:9,];num3=c(6:9)

x_ave=apply(X,2,mean)

x_ave1=apply(X1,2,mean)

x_ave2=apply(X2,2,mean)

x_ave3=apply(X3,2,mean)

B_vec=cbind(matrix(rep(x_ave1-x_ave,length(num1)),ncol=length(num1)),matrix(rep(x_ave2-x_ave,length(num2)),ncol=length(num2)),matrix(rep(x_ave3-x_ave,length(num3)),ncol=length(num3)))

B=B_vec%*%t(B_vec)

eigen=eigen(solve(t)%*%B)

lam=eigen$values

v=eigen$vectors

mat=solve(t)%*%B

#本の結果での確認

mat_dush=mat;diag(mat_dush)=diag(mat_dush)-lam[1]

vec=c(-9.7434,-4.5811,-1.5811,1)

#3群なので、固有ベクトルは3-1＝2で第2固有ベクトルまで使用する

v_vec=v[,1:2]

c_vec=c()

for(j in 1:ncol(v_vec)){
  
vector=c(v[,j])

val=t(vector)%*%t%*%t(t(vector))

c_val=sqrt(n/val)

c_vec=c(c_vec,c_val)

v_vec[,j]=as.numeric(c_val)*v_vec[,j]
  
}

a=v_vec

a=rbind(a[1:2,],rep(0,ncol(a)),a[3:4,],rep(0,ncol(a)))

#XX%*%a

ns=apply(XX,2,sum)

additiv_w=t(cbind(matrix(rep(apply((a*c(ns))[1:3,],2,sum)/n,3),nrow=ncol(a)),matrix(rep(apply((a*c(ns))[4:6,],2,sum)/n,3),nrow=ncol(a))))

b=a-additiv_w

plot(XX%*%b)

```