数量化2類（多変量解析入門Ⅱ　河田至商先生）

```{r}



#数量化２類

library(dplyr)

data1=data.frame(terms=c("自社製を購入"),価格_考慮=c(1,0,1,0),価格_不考慮=c(0,1,0,1),概観_考慮=c(1,1,0,0),概観_不考慮=c(0,0,1,1),性能_検討した=c(0,1,1,1),性能_一応検討=c(1,0,0,0),性能_検討しない=c(0,0,0,0))

data2=data.frame(terms=c("他社製を購入"),価格_考慮=c(1,0,0,0,0,0),価格_不考慮=c(0,1,1,1,1,1),概観_考慮=c(1,0,0,1,1,1),概観_不考慮=c(0,1,1,0,0,0),性能_検討した=c(0,0,0,0,0,0),性能_一応検討=c(0,0,1,0,1,0),性能_検討しない=c(1,1,0,1,0,1))

data=rbind(data1,data2)

split=c(2,2,3)

X=data[,-1]

class=unique(data$terms)

X_B=array(0,dim=c(1,ncol(X)-length(split)))

num=c(1,cumsum(split)[1:(length(split)-1)]+1)

for(i in 1:length(class)){
  
ave=apply(X[data$terms==class[i],],2,mean)

sub=t(matrix(rep(ave[-num],sum(data$terms==class[i])),ncol=sum(data$terms==class[i])))  
 
X_B=rbind(X_B,sub)
 
}

X_B=tail(X_B,nrow(X_B)-1)


X_bar=t(matrix(rep(c(apply(X[,-num],2,mean)),nrow(X)),nrow=ncol(X)-length(split)))

XX=matrix(as.numeric(unlist(X)),nrow=nrow(X),ncol=ncol(X))

X=X[,-num]

X=matrix(as.numeric(unlist(X)),nrow=nrow(X),ncol=ncol(X))

S=t(X-X_bar)%*%(X-X_bar);S_B=t(X_B-X_bar)%*%(X_B-X_bar)

#Sを三角行列に分解

F_matrix=t(chol(S));F_inverse=chol(S)

eta2=(eigen(solve(F_matrix)%*%S_B%*%solve(t(F_matrix)))$values)

eigen_vectors=eigen(solve(F_matrix)%*%S_B%*%solve(t(F_matrix)))$vectors

A=solve(t(F_matrix))%*%eigen_vectors[,1]

A=c(0,A)

for(i in 2:length(num)){
  
A=append(A, c(0), after=num[i]-1)  
  
}

result=data.frame(term=data$term,predict=XX%*%A)

b=c()

for(j in 1:length(A)){

if(j %in% c(1:2)){  

b=c(b,A[j]-sum(A[1:2]*XX_bar[1,1:2]))  

}else{

if(j %in% c(3:4)){

b=c(b,A[j]-sum(A[3:4]*XX_bar[1,3:4]))  

}else{

b=c(b,A[j]-sum(A[5:7]*XX_bar[1,5:7]))  

}  

} 

}

result=data.frame(terms=data$terms,predict=0)

for(j in 1:nrow(result)){

result$predict[j]=sum(b*XX[j,])  

}

plot(c(1:length(result$predict[result$terms=="自社製を購入"])),result$predict[result$terms=="自社製を購入"],type="p",col=2,xlab="",ylab="自社製を購入",ylim=c(min(result$predict),max(result$predict)))



plot(c(1:length(result$predict[result$terms=="他社製を購入"])),result$predict[result$terms=="他社製を購入"],type="p",col=3,xlab="",ylab="他社製を購入",ylim=c(min(result$predict),max(result$predict)))

```

```{r}
#数量化２類
data1=data.frame(terms=c("自社製を購入"),価格_考慮=c(1,0,1,0),価格_不考慮=c(0,1,0,1),概観_考慮=c(1,1,0,0),概観_不考慮=c(0,0,1,1),性能_検討した=c(0,1,1,1),性能_一応検討=c(1,0,0,0),性能_検討しない=c(0,0,0,0))
data2=data.frame(terms=c("他社製を購入"),価格_考慮=c(1,0,0,0,0,0),価格_不考慮=c(0,1,1,1,1,1),概観_考慮=c(1,0,0,1,1,1),概観_不考慮=c(0,1,1,0,0,0),性能_検討した=c(0,0,0,0,0,0),性能_一応検討=c(0,0,1,0,1,0),性能_検討しない=c(1,1,0,1,0,1))
data=rbind(data1,data2)
split=c(2,2,3)
X=data[,-1]
class=unique(data$terms)
X_B=array(0,dim=c(1,ncol(X)-length(split)))
num=c(1,cumsum(split)[1:(length(split)-1)]+1)
for(i in 1:length(class)){
ave=apply(X[data$terms==class[i],],2,mean)
sub=t(matrix(rep(ave[-num],sum(data$terms==class[i])),ncol=sum(data$terms==class[i])))  
X_B=rbind(X_B,sub)
}
X_B=tail(X_B,nrow(X_B)-1)
X_bar=t(matrix(rep(c(apply(X[,-num],2,mean)),nrow(X)),nrow=ncol(X)-length(split)))
XX_bar=t(matrix(rep(c(apply(X,2,mean)),nrow(X)),nrow=ncol(X)))
XX=matrix(as.numeric(unlist(X)),nrow=nrow(X),ncol=ncol(X))
X=X[,-num]
X=matrix(as.numeric(unlist(X)),nrow=nrow(X),ncol=ncol(X))
S=t(X-X_bar)%*%(X-X_bar);S_B=t(X_B-X_bar)%*%(X_B-X_bar)

#Sを三角行列に分解(cholesky分解)
mat=S
A=mat
n=nrow(mat)
U=array(0,dim=dim(A))
for(k in 1:n){
U[k,k]=sqrt(A[k,k])
if(k!=n){
for(i in (k+1):n){
U[k,i]=A[k,i]/U[k,k]  
}
for(i in (k+1):n){
for(j in (k+1):n){  
A[i,j]=A[i,j]-U[k,i]*U[k,j]    
}}
}}
F_matrix=t(U);F_inverse=U

#QR法　
#固有値(A_matの対角成分)と固有ベクトル(Q)
A=solve(F_matrix)%*%S_B%*%solve(t(F_matrix))
A_mat=A
ite=100
Q=diag(rep(1,nrow(A_mat)))
for(l in 1:ite){
q_mat=array(0,dim=c(nrow(A_mat),ncol(A_mat)))
r_mat=array(0,dim=c(nrow(A_mat),ncol(A_mat)))
q1=A_mat[,1]/sqrt(sum(A_mat[,1]^2))
q_mat[,1]=q1
r_mat[1,1]=sqrt(sum(A_mat[,1]^2))
for(j in 2:ncol(A_mat)){
u=A_mat[,j]-apply(t(t(q_mat[,1:(j-1)])*c(t(q_mat[,1:(j-1)])%*%A_mat[,j])),1,sum)
r_mat[j,j]=sqrt(sum(u^2))
q_mat[,j]=u/sqrt(sum(u^2))
}
r_mat2=t(A_mat)%*%q_mat
r_mat2[upper.tri(r_mat2)]=0
diag(r_mat2)=diag(r_mat)
r_mat=t(r_mat2)
A_mat=r_mat%*%q_mat
Q=Q%*%q_mat
}


eta2=diag(A_mat)
eigen_vectors=Q
v=solve(t(F_matrix))%*%eigen_vectors[,1]
v=c(0,v)
for(i in 2:length(num)){
v=append(v, c(0), after=num[i]-1)  
}
result=data.frame(term=data$term,predict=XX%*%v)
b=c()
for(j in 1:length(num)){
if(j!=length(num)){  
number=c(1:length(v))[num[j]:(num[j+1]-1)]  
}else{
number=c(1:length(v))[num[j]:length(v)]
}
b=c(b,v[number]-sum(v[number]*XX_bar[1,number]))
}
result=data.frame(terms=data$terms,predict=XX%*%b)
plot(c(1:length(result$predict[result$terms=="自社製を購入"])),result$predict[result$terms=="自社製を購入"],type="p",col=2,xlab="",ylab="自社製を購入",ylim=c(min(result$predict),max(result$predict)))

plot(c(1:length(result$predict[result$terms=="他社製を購入"])),result$predict[result$terms=="他社製を購入"],type="p",col=3,xlab="",ylab="他社製を購入",ylim=c(min(result$predict),max(result$predict)))
```


```{r}
#数量化２類
data1=data.frame(terms=c("自社製を購入"),価格_考慮=c(1,0,1,0),価格_不考慮=c(0,1,0,1),概観_考慮=c(1,1,0,0),概観_不考慮=c(0,0,1,1),性能_検討した=c(0,1,1,1),性能_一応検討=c(1,0,0,0),性能_検討しない=c(0,0,0,0))
data2=data.frame(terms=c("他社製を購入"),価格_考慮=c(1,0,0,0,0,0),価格_不考慮=c(0,1,1,1,1,1),概観_考慮=c(1,0,0,1,1,1),概観_不考慮=c(0,1,1,0,0,0),性能_検討した=c(0,0,0,0,0,0),性能_一応検討=c(0,0,1,0,1,0),性能_検討しない=c(1,1,0,1,0,1))
data=rbind(data1,data2)
split=c(2,2,3)
X=data[,-1]
class=unique(data$terms)
X_B=array(0,dim=c(1,ncol(X)-length(split)))
num=c(1,cumsum(split)[1:(length(split)-1)]+1)
for(i in 1:length(class)){
ave=apply(X[data$terms==class[i],],2,mean)
sub=t(matrix(rep(ave[-num],sum(data$terms==class[i])),ncol=sum(data$terms==class[i])))  
X_B=rbind(X_B,sub)
}
X_B=tail(X_B,nrow(X_B)-1)
X_bar=t(matrix(rep(c(apply(X[,-num],2,mean)),nrow(X)),nrow=ncol(X)-length(split)))
XX_bar=t(matrix(rep(c(apply(X,2,mean)),nrow(X)),nrow=ncol(X)))
XX=matrix(as.numeric(unlist(X)),nrow=nrow(X),ncol=ncol(X))
X=X[,-num]
X=matrix(as.numeric(unlist(X)),nrow=nrow(X),ncol=ncol(X))
S=t(X-X_bar)%*%(X-X_bar);S_B=t(X_B-X_bar)%*%(X_B-X_bar)
#Sを三角行列に分解
F_matrix=t(chol(S));F_inverse=chol(S)

#QR法　
#固有値(A_matの対角成分)と固有ベクトル(Q)
A=solve(F_matrix)%*%S_B%*%solve(t(F_matrix))
A_mat=A
ite=100
Q=diag(rep(1,nrow(A_mat)))
for(l in 1:ite){
q_mat=array(0,dim=c(nrow(A_mat),ncol(A_mat)))
r_mat=array(0,dim=c(nrow(A_mat),ncol(A_mat)))
q1=A_mat[,1]/sqrt(sum(A_mat[,1]^2))
q_mat[,1]=q1
r_mat[1,1]=sqrt(sum(A_mat[,1]^2))
for(j in 2:ncol(A_mat)){
u=A_mat[,j]-apply(t(t(q_mat[,1:(j-1)])*c(t(q_mat[,1:(j-1)])%*%A_mat[,j])),1,sum)
r_mat[j,j]=sqrt(sum(u^2))
q_mat[,j]=u/sqrt(sum(u^2))
}
r_mat2=t(A_mat)%*%q_mat
r_mat2[upper.tri(r_mat2)]=0
diag(r_mat2)=diag(r_mat)
r_mat=t(r_mat2)
A_mat=r_mat%*%q_mat
Q=Q%*%q_mat
}


eta2=diag(A_mat)
eigen_vectors=Q
v=solve(t(F_matrix))%*%eigen_vectors[,1]
v=c(0,v)
for(i in 2:length(num)){
v=append(v, c(0), after=num[i]-1)  
}
result=data.frame(term=data$term,predict=XX%*%v)
b=c()
for(j in 1:length(num)){
if(j!=length(num)){  
number=c(1:length(v))[num[j]:(num[j+1]-1)]  
}else{
number=c(1:length(v))[num[j]:length(v)]
}
b=c(b,v[number]-sum(v[number]*XX_bar[1,number]))
}
result=data.frame(terms=data$terms,predict=XX%*%b)
plot(c(1:length(result$predict[result$terms=="自社製を購入"])),result$predict[result$terms=="自社製を購入"],type="p",col=2,xlab="",ylab="自社製を購入",ylim=c(min(result$predict),max(result$predict)))

plot(c(1:length(result$predict[result$terms=="他社製を購入"])),result$predict[result$terms=="他社製を購入"],type="p",col=3,xlab="",ylab="他社製を購入",ylim=c(min(result$predict),max(result$predict)))
```


```{r}
#数量化２類
data1=data.frame(terms=c("自社製を購入"),価格_考慮=c(1,0,1,0),価格_不考慮=c(0,1,0,1),概観_考慮=c(1,1,0,0),概観_不考慮=c(0,0,1,1),性能_検討した=c(0,1,1,1),性能_一応検討=c(1,0,0,0),性能_検討しない=c(0,0,0,0))
data2=data.frame(terms=c("他社製を購入"),価格_考慮=c(1,0,0,0,0,0),価格_不考慮=c(0,1,1,1,1,1),概観_考慮=c(1,0,0,1,1,1),概観_不考慮=c(0,1,1,0,0,0),性能_検討した=c(0,0,0,0,0,0),性能_一応検討=c(0,0,1,0,1,0),性能_検討しない=c(1,1,0,1,0,1))
data=rbind(data1,data2)
split=c(2,2,3)
X=data[,-1]
class=unique(data$terms)
X_B=array(0,dim=c(1,ncol(X)-length(split)))
num=c(1,cumsum(split)[1:(length(split)-1)]+1)
for(i in 1:length(class)){
ave=apply(X[data$terms==class[i],],2,mean)
sub=t(matrix(rep(ave[-num],sum(data$terms==class[i])),ncol=sum(data$terms==class[i])))  
X_B=rbind(X_B,sub)
}
X_B=tail(X_B,nrow(X_B)-1)
X_bar=t(matrix(rep(c(apply(X[,-num],2,mean)),nrow(X)),nrow=ncol(X)-length(split)))
XX_bar=t(matrix(rep(c(apply(X,2,mean)),nrow(X)),nrow=ncol(X)))
XX=matrix(as.numeric(unlist(X)),nrow=nrow(X),ncol=ncol(X))
X=X[,-num]
X=matrix(as.numeric(unlist(X)),nrow=nrow(X),ncol=ncol(X))
S=t(X-X_bar)%*%(X-X_bar);S_B=t(X_B-X_bar)%*%(X_B-X_bar)
#Frame 法
A=solve(S)%*%S_B
row=nrow(A)
a=c(-sum(diag(A)))
A2=A
for(j in 1:(row-1)){
As=A2+diag(a[length(a)],row)
A2=A%*%As
a=c(a,-sum(diag(A2))/(j+1))
}
lam=Re(sort(polyroot(c(a[length(a):1],1)),decreasing=T)) 
Q=function(x){
mat=diag(x^(row-1),row)
a=c(-sum(diag(A)))
A2=A
for(j in 1:(row-1)){
As=A2+diag(a[length(a)],row)
mat=mat+As*x^(row-1-j)
A2=A%*%As
a=c(a,-sum(diag(A2))/(j+1))
}
return(mat)
}
v=array(0,dim=dim(A))
for(j in 1:row){
v[,j]=Q(lam[j])[,1]
}
v=t(t(v)/sqrt(apply(v^2,2,sum)))[,1]
v=c(0,v)
for(i in 2:length(num)){
v=append(v, c(0), after=num[i]-1)  
}
result=data.frame(term=data$term,predict=XX%*%v)
b=c()
for(j in 1:length(num)){
if(j!=length(num)){  
number=c(1:length(v))[num[j]:(num[j+1]-1)]  
}else{
number=c(1:length(v))[num[j]:length(v)]
}
b=c(b,v[number]-sum(v[number]*XX_bar[1,number]))
}
result=data.frame(terms=data$terms,predict=XX%*%b)
plot(c(1:length(result$predict[result$terms=="自社製を購入"])),result$predict[result$terms=="自社製を購入"],type="p",col=2,xlab="",ylab="自社製を購入",ylim=c(min(result$predict),max(result$predict)))

plot(c(1:length(result$predict[result$terms=="他社製を購入"])),result$predict[result$terms=="他社製を購入"],type="p",col=3,xlab="",ylab="他社製を購入",ylim=c(min(result$predict),max(result$predict)))
```


```{r}
#数量化２類(2次元以上)
library(dummies)
data=data.frame(terms=c(rep(1,10),rep(2,10),rep(3,10),rep(4,10)),akarui=dummy(c(2,1,2,1,1,2,1,1,2,1,2,1,1,2,1,1,1,2,1,1,2,2,1,1,1,1,1,1,2,2,2,1,1,1,2,1,2,1,1,1)),nonbiri=dummy(c(1,rep(2,9),1,rep(2,7),1,2,1,2,1,1,2,2,2,1,1,2,1,1,2,2,1,1,2,1,1,1)),kibunya=dummy(c(2,2,1,2,2,1,2,2,2,2,2,1,2,2,2,1,1,2,2,2,1,1,1,2,2,2,1,1,1,1,1,2,2,2,2,1,2,2,1,2)),kikubari=dummy(c(1,2,1,2,1,2,1,1,1,1,1,1,2,2,2,2,2,1,2,2,1,2,2,2,2,1,2,2,2,1,1,2,1,2,2,1,2,1,2,2)),yumoa=dummy(c(1,2,2,2,1,2,2,2,2,2,2,2,1,1,2,2,2,2,2,2,1,2,2,1,2,1,1,2,1,2,1,1,1,2,1,1,2,1,1,2)),cool=dummy(c(2,2,2,2,2,2,1,2,2,1,1,2,1,2,2,2,2,1,2,2,2,1,1,1,2,1,1,1,2,1,2,2,2,1,2,2,1,2,2,2)),shintyou=dummy(c(1,2,1,1,2,1,2,1,1,1,2,2,2,1,2,2,2,2,2,2,1,2,2,2,1,2,2,1,1,2,1,2,2,2,2,2,2,2,2,2)),syakouteki=dummy(c(2,2,2,2,1,2,2,2,2,2,2,2,2,1,2,2,1,2,2,2,1,2,1,2,1,1,1,1,2,1,1,2,2,1,2,2,1,2,1,1)),koseiteki=dummy(c(1,2,2,2,2,1,2,2,2,2,1,1,2,1,1,2,2,1,1,1,2,2,2,1,1,2,2,2,2,2,2,1,2,2,2,2,2,2,2,2)))
X=data[,-1]
split=c(2,2,2,2,2,2,2,2,2)
class=unique(data$terms)
X_B=array(0,dim=c(1,ncol(X)-length(split)))
num=c(1,cumsum(split)[1:(length(split)-1)]+1)
for(i in 1:length(class)){
ave=apply(X[data$terms==class[i],],2,mean)
sub=t(matrix(rep(ave[-num],sum(data$terms==class[i])),ncol=sum(data$terms==class[i])))  
X_B=rbind(X_B,sub)
}
X_B=tail(X_B,nrow(X_B)-1)
X_bar=t(matrix(rep(c(apply(X[,-num],2,mean)),nrow(X)),nrow=ncol(X)-length(split)))
XX_bar=t(matrix(rep(c(apply(X,2,mean)),nrow(X)),nrow=ncol(X)))
XX=matrix(as.numeric(unlist(X)),nrow=nrow(X),ncol=ncol(X))
X=X[,-num]
X=matrix(as.numeric(unlist(X)),nrow=nrow(X),ncol=ncol(X))
S=t(X-X_bar)%*%(X-X_bar);S_B=t(X_B-X_bar)%*%(X_B-X_bar)
#Sを三角行列に分解
F_matrix=t(chol(S));F_inverse=chol(S)

#QR法　
#固有値(A_matの対角成分)と固有ベクトル(Q)
A=solve(F_matrix)%*%S_B%*%solve(t(F_matrix))
A_mat=A
ite=100
Q=diag(rep(1,nrow(A_mat)))
for(l in 1:ite){
q_mat=array(0,dim=c(nrow(A_mat),ncol(A_mat)))
r_mat=array(0,dim=c(nrow(A_mat),ncol(A_mat)))
q1=A_mat[,1]/sqrt(sum(A_mat[,1]^2))
q_mat[,1]=q1
r_mat[1,1]=sqrt(sum(A_mat[,1]^2))
for(j in 2:ncol(A_mat)){
u=A_mat[,j]-apply(t(t(q_mat[,1:(j-1)])*c(t(q_mat[,1:(j-1)])%*%A_mat[,j])),1,sum)
r_mat[j,j]=sqrt(sum(u^2))
q_mat[,j]=u/sqrt(sum(u^2))
}
r_mat2=t(A_mat)%*%q_mat
r_mat2[upper.tri(r_mat2)]=0
diag(r_mat2)=diag(r_mat)
r_mat=t(r_mat2)
A_mat=r_mat%*%q_mat
Q=Q%*%q_mat
}


eta2=diag(A_mat)
eigen_vectors=Q
v=solve(t(F_matrix))%*%eigen_vectors[,1:3]

v=rbind(rep(0,ncol(v)),v)
for(i in 2:length(num)){
v=rbind(rbind(v[1:(num[i]-1),],rep(0,ncol(v))),v[num[i]:nrow(v),])
}
result=data.frame(term=data$term,predict=XX%*%v)

b=array(0,dim=dim(v))
for(j in 1:length(num)){
if(j!=length(num)){  
number=c(1:nrow(v))[num[j]:(num[j+1]-1)]  
}else{
number=c(1:nrow(v))[num[j]:nrow(v)]
}
b[number,]=t(t(v[number,])-c(t(v[number,])%*%XX_bar[1,number]))
}

library(scatterplot3d)

result=data.frame(terms=data$terms,predict=XX%*%b)

colors <- c("#ff0000", "#0000ff", "#008000","#000000")
colors_pred <- colors[as.numeric(data$terms)]
scatterplot3d(result[,2],result[,3],result[,4],color=colors_pred,pch=20)

library(dplyr)

summary=result%>%group_by(terms)%>%summarise(predict.1=mean(predict.1),predict.2=mean(predict.2),predict.3=mean(predict.3))


d=array(0,dim=c(nrow(result),length(class)))


for(j in class){
  
d[,j]=apply((t(result[,-1])-unlist(summary[j,-1]))^2,2,sum) 
}

pred_class=c()

for(j in 1:nrow(d)){
  
pred_class=c(pred_class,sum(ifelse(d[j,]==min(d[j,]),1,0)*c(1:length(class))))  
  
}


```