数理決定法入門-キャンパスのOR　朝倉書店　今野浩先生

```r

#p.40 通学ルート決定問題

f=c("最寄り駅までの所要時間","乗車時間","乗車金額片道","乗り換え回数","総合所要時間","電車の混雑度")

R=matrix(c(10,62,760,11210,3,7,85.5,6.5,10,72,730,11930,2,8,97,7,10,72,710,9750,3,6,99,11,10,64,710,12460,3,7,89,8,7,74,870,12720,3,10,97.5,6.5),nrow=8)

rownames(R)=c("最寄り駅までの所要時間","乗車時間","乗車片道切符金額","1カ月定期金額","乗換回数","乗換時間","総合所要時間","乗り換え待ち時間")

A=matrix(c(1,3,6,4,8,4,0,1,0.5,5,4,3,0,0,1,1/3,3,2,0,0,0,1,1/3,1/5,rep(0,4),1,1/6,rep(0,5),1),nrow=6)

invA=t(1/A)
invA[lower.tri(invA)]=0
diag(invA)=0
A=A+invA

#べき乗法
u=rep(1,nrow(A))/nrow(A)

ite=100
for(l in 1:ite){
upre=u
v=A%*%u
u=v/sum(v)
a=sum(v)
print(sum(abs(u-upre)))
}

#重要度係数
w=u
names(w)=f



#ルート(R)の相対評価 p.46

#要因f1に対するルートRi,Rjの相対評価の行列
B1=cbind(rbind(array(1,dim=c(4,4)),rep(4,4)),c(rep(1/4,4),1))

B2=array(0,dim=c(5,5))
B2[lower.tri(B2)]=c(1/6,1/6,1/2,1/7,1,5,1/2,5,1/2,1/6)
invB2=t(1/B2)
invB2[lower.tri(invB2)]=0
diag(invB2)=1
B2=B2+invB2

B3=array(0,dim=c(5,5))
B3[lower.tri(B3)]=c(1/2,4,1/4,1/4,6,1/2,1/3,1/7,1/8,1/2)
invB3=t(1/B3)
invB3[lower.tri(invB3)]=0
diag(invB3)=1
B3=B3+invB3

B4=array(0,dim=c(5,5))
B4[lower.tri(B4)]=c(4,2,1,1/4,1/3,1/4,1/7,0.5,1/5,1/4)
invB4=t(1/B4)
invB4[lower.tri(invB4)]=0
diag(invB4)=1
B4=B4+invB4

B5=array(0,dim=c(5,5))
B5[lower.tri(B5)]=c(1/7,1/8,1/3,1/7,0.5,5,1/2,6,3,1/6)
invB5=t(1/B5)
invB5[lower.tri(invB5)]=0
diag(invB5)=1
B5=B5+invB5

B6=array(0,dim=c(5,5))
B6[lower.tri(B6)]=c(5,3,3,9,1/3,1/3,5,1,6,6)
invB6=t(1/B6)
invB6[lower.tri(invB6)]=0
diag(invB6)=1
B6=B6+invB6


#べき乗法
eigen_calc=function(AA){
u=rep(1,nrow(AA))/nrow(AA)
ite=100
for(l in 1:ite){
upre=u
v=AA%*%u
u=v/sum(v)
a=sum(v)
print(paste0("誤差:",sum(abs(u-upre))))
}
return(c(a,u))
}

#べき乗法による結果
f1=eigen_calc(B1)[-1]
f2=eigen_calc(B2)[-1]
f3=eigen_calc(B3)[-1]
f4=eigen_calc(B4)[-1]
f5=eigen_calc(B5)[-1]
f6=eigen_calc(B6)[-1]

fij=rbind(f1,f2,f3,f4,f5,f6)
res=cbind(w,fij)　#表3.3
colnames(res)=c("w",paste0("R",c(1:5)))

#総合評価点
s=t(fij)%*%w


```

```r

#p.64 ファクターモデル
#かなり計算遅いし、なんかな～

J=4
K=3
vj=sample(1:20,K)
tj=sample(1:20,J)
vrs=matrix(sample(1:20,K^2,replace=T),nrow=K)*10
h=10^(-5)

data("EuStockMarkets")
R=matrix(as.numeric(unlist(EuStockMarkets)),nrow=nrow(EuStockMarkets))
R=R[1:100,]

V=function(a,b){
V=array(0,dim=c(J,J))
for(i in 1:J){
for(j in 1:J){
V[i,j]=c(b[i,])%*%vrs%*%c(b[j,]) 
}}
diag(V)=diag(V)+c(tj)  
return(V)
}

loglik=function(a,b){
Ra=t(R)-c(a)
loglik=0
V=V(a,b)
for(i in 1:nrow(R)){
loglik=loglik-as.numeric(Ra[,i])%*%solve(V)%*%c(Ra[,i])
}
detV=det(V)
if(detV>0){
loglik=loglik-nrow(R)*log((det(V))+1)/4
}
return(as.numeric(loglik))
}

f=function(x){
x=matrix(x,nrow=J)  
a=x[,1]
b=x[,-1]
Ra=t(R)-c(a)
Vmat=V(a,b)
dalpha=-apply((solve(Vmat)+t(solve(Vmat)))%*%Ra,1,sum)/2

Vb=function(p){
value=0  
for(i in 1:nrow(R)){
value=value-as.numeric(Ra[,i])%*%solve(V(a,p))%*%c(Ra[,i])/2
}
return(value)
}

dVbeta=array(0,dim=dim(b))  
for(j in 1:J){
for(k in 1:K){
b_sub=b
b_sub[j,k]=b_sub[j,k]+h
dVbeta[j,k]=(Vb(b_sub)-Vb(b))/h-sum(diag(solve(V(a,b))%*%(V(a,b_sub)-V(a,b))/h))
}} 
  
return(c(cbind(dalpha,dVbeta)))
}

beta=solve(matrix(rep(diag(vrs),J),nrow=K)%*%t(matrix(rep(diag(vrs),J),nrow=K))+diag(1,K))%*%(matrix(rep(diag(vrs),J),nrow=K)%*%(diag(cov(R))-tj^2))

#準ニュートン法
X=c(apply(R,2,mean),rep(beta,J))　#初期パラメータ
eta=10^(-3)
H=diag(f(X))
loglikvec=c()
paramdiffvec=c()

#更新が合計で10^(-6)より大きければ計算する
while(sum(abs(H%*%f(X))>10^(-6)))({
X_pre=X  
X=X-eta*H%*%f(X)  
s=X-X_pre
y=f(X)-f(X_pre)
if(as.numeric(t(s)%*%H%*%y)==0){
H=H+((s-H%*%y)/as.numeric(t(s)%*%H%*%y+0.0001))%*%t(s)%*%H
}else{
H=H+((s-H%*%y)/as.numeric(t(s)%*%H%*%y))%*%t(s)%*%H  
}
#print(sum(abs(H%*%f(X))))
xx=matrix(X,nrow=J)
alpha=xx[1:J]
beta=matrix(xx[-c(1:J)],nrow=J)
loglikvec=c(loglikvec,loglik(alpha,beta))
paramdiffvec=c(paramdiffvec,sum(abs(f(X))))
print(loglik(alpha,beta))
})

X=matrix(X,nrow=J)
alpha=X[1:J]
beta=matrix(X[-c(1:J)],nrow=J)
plot(c(1:length(loglikvec)),loglikvec,main="loglik")
plot(c(1:length(paramdiffvec),paramdiffvec,main="paramdiff"))

```