AHPとコンジョイント分析　木下栄三・大野栄治先生　現代数学社

```{r}
#コンジョイント分析

#数量化1類

data=data.frame(販売率=c(3.0,1.8,1.5,3.3,2.2,2.0,3.5,2.0,1.7,2.3),大都市=c(1,0,0,1,1,1,1,0,1,1),地方都市=c(0,1,1,0,0,0,0,1,0,0),温暖=c(0,1,0,0,0,0,1,1,0,0),普通=c(1,0,1,1,0,0,0,0,0,0),寒冷=c(0,0,0,0,1,1,0,0,1,1))

X=data[,colnames(data)!="温暖"];X=X[,2:ncol(X)]

x=array(0,dim=c(length(colnames(X)),length(colnames(X))))

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


  x[i,j]=sum(X[,i]*X[,j])


}
}

xy=array(0,dim=c(ncol(X),1))

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

xy[i]=sum(X[,i]*data$販売率) 

}

A=solve(x)%*%xy

y=c()

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

y=c(y,sum(X[j,]*A))  

}

U=y



#MONANOVA

for(k in 2:length(y)){
  
if(U[k]<U[k-1]){
  
MU_data=data.frame(num=1:(k-1),MU=0,sig=0)

for(l in 1:(k-1)){  
  
MU=sum(U[k]+U[(k-1):(k-l)])/(l+1)

if(k-l-1==0){

MU_data$MU[l]=MU;MU_data$sig[l]=ifelse(MU>0,1,0)

}else{

MU_data$MU[l]=MU;MU_data$sig[l]=ifelse(MU>U[(k-l-1)],1,0)

}  

}

no=min(MU_data$num[MU_data$sig==1])

MU=MU_data$MU[no]

U[k:(k-no)]=MU

}

U_ave=c()  
    
for(l in 1:length(U)){
  
U_ave=c(U_ave,mean(U[1:l]))  
  
}  
  

eta=sqrt(sum((U-U_ave)^2)/sum((U-mean(U))^2))

print(eta)
  
}  

```

```{r}

#数量化1類　準ニュートン法

#fiacco,mccormick(1968),murtagh,sargent(1969)

data=data.frame(販売率=c(3.0,1.8,1.5,3.3,2.2,2.0,3.5,2.0,1.7,2.3),大都市=c(1,0,0,1,1,1,1,0,1,1),地方都市=c(0,1,1,0,0,0,0,1,0,0),温暖=c(0,1,0,0,0,0,1,1,0,0),普通=c(1,0,1,1,0,0,0,0,0,0),寒冷=c(0,0,0,0,1,1,0,0,1,1))

XX=data[,2:ncol(data)]

Y=data$販売率

XX=matrix(as.matrix(XX),ncol=ncol(XX),nrow=nrow(XX))

f=function(x){

f=-2*t(XX)%*%(Y-XX%*%x)

return(c(f))

}

X=rep(1,ncol(XX))

H=diag(1,length(X))

ite=10000

eta=10^(-2)

h=0.01

for(l in 1:ite){

f_val=f(X)

sigma=t(t(-H%*%f_val))

X=X+eta*sigma

y=f(X)-f_val

H=H+(eta*sigma-H%*%y)%*%t(eta*sigma-H%*%y)/as.numeric(t(y)%*%(eta*sigma-H%*%y))

print(f(X))
 
}





```

```{r}

#効用推定法


X=as.matrix(data[,2:ncol(data)]);Y=c(data$販売率)

a=rep(1,ncol(X));b=rep(1,nrow(X))

ite=10000;delta=0.01

for(j in 1:ite){
  
a_pre=a;b_pre=b  
  
A=sum((U-U_ave)^2);B=sum((U-mean(U))^2)

eta=sqrt(sum((U-U_ave)^2)/sum((U-mean(U))^2))

  
for(i in 1:length(a)){

da=2*sum(X[,i]*(U-U_ave))

db=2*sum(X[,i]*(U-mean(U)))

de=(da*eta/A-db*eta/B)/2

a[i]=a[i]-delta*de
  
}  
  
for(i in 1:length(b)){

da=2*sum((U-U_ave)*Y)

db=2*sum(U-mean(U)*Y)

de=(da*eta/A-db*eta/B)/2

b[i]=b[i]-delta*de

} 

U=c(X%*%a)+c(b*Y)

U_ave=c()

for(l in 1:length(U)){
  
U_ave=c(U_ave,mean(U[1:l]))  
  
}

print(sum(abs(c(a_pre,b_pre)-c(a,b))))
  
}

```



```{r}

#AHPとコンジョイント分析　木下栄蔵・大野栄治先生　　現代数学社

#誤差モデル p.54

X=matrix(c(1,1/5,1/9,5,1,1/3,9,3,1),ncol=3)

X=log(X)

n=nrow(X)

w=c()

k=10

for(i in 1:n){

if(i==1){

w=c(w,sum(X[1,-1])/n+k/n)

}

if(i!=1 & i!=n){

w=c(w,(-sum(X[1:(i-1),i])+sum(X[i,(i+1):n]))/n+k/n)

}

if(i==n){

w=c(w,-sum(X[,n])/n+k/n)

}

}

mat=array(0,dim=c(n,n))

for(i in 1:n){
for(j in 1:n){

mat[i,j]=w[i]-w[j]

}}

#べき乗法 p.58

X=matrix(c(1,1/5,1/9,5,1,1/3,9,3,1),ncol=3)

n=nrow(X)

b=n

v=rep(1/n,n)

ite=100

for(l in 1:ite){

v_pre=v

u=X%*%v

b=sum(u)

v=u/b

cost=sum(abs(v-v_pre))

print(cost)

}







#コンジョイント分析 p.124

data=data.frame(ID=c("A","B","C","D"),price=c("安い","高い","安い","高い"),fuel=c("安い","安い","高い","高い"),order=c(1:4),info=c(3,2,4,1),effect=c(4:1))

delta=cbind(dummy(data$price),dummy(data$fuel))

X=cbind(data$order,data$info);U=data$effect

D=cbind(delta,X)

Z=U

b=solve(t(D)%*%D+diag(rep(1,ncol(D))))%*%t(D)%*%Z

ite=3*10^5

eta=10^(-4)

for(s in 1:ite){

Z_hat=D%*%b

Z_hat_ave=rep(mean(Z_hat),nrow(D))

V=as.numeric(t(c(Z-D%*%b))%*%t(t(c(Z-D%*%b))))

W=as.numeric(t(c(D%*%b-Z_hat_ave))%*%t(t(c(D%*%b-Z_hat_ave))))

S=sqrt(V/W)

g=-(S/W)*t(D)%*%((Z-Z_hat)/(S^2)+Z-Z_hat_ave)

b=b-eta*g

print(S)

}


```