交通ネットワークの均衡分析

```{r}

#p.75

m=10;n=20

a=(array(rnorm(m*n),dim=c(m,n)))

beta=rep(0.01,m)

cost=function(sit,bs){
  
V2=t(a)%*%bs  

P_k2=exp(sit*V2)/sum(exp(sit*V2))

return(sum(log(exp(sit*V2)))/sit)
  
}

ite=100000

sita=0.01

h=0.001


for(l in 1:ite){
  
log_lik_pre=cost(sita,beta)
  
V=t(a)%*%beta  
 
P_k=exp(sita*V)/sum(exp(sita*V)) 

db=a%*%P_k;db2=((a^2)%*%P_k)-(a%*%P_k)^2

dsita=(cost(sita+h,beta)-cost(sita,beta))/h

dsita2=(cost(sita+h,beta)-2*cost(sita,beta)+cost(sita-h,beta))/(h^2)
 
H=diag(c(dsita2,db2));G=diag(1/c(db2));dg=c(db)

params=c(beta)+h*G%*%dg

beta=params

print(log(sum(exp(sita*V)))/sita)  

}

```

```{r}

#p.75 applied

m=10;n=20

a=(array(rnorm(m*n),dim=c(m,n)))

beta=rep(0.01,m)

cost=function(sit,bs){
  
V2=t(a)%*%bs  

P_k2=exp(sit*V2)/sum(exp(sit*V2))

return(sum(log((P_k2)))/sit)
  
}

ite=100000

sita=1

h=0.001

log_lik_pre=-10;log_lik=-9

for(l in 1:ite){
  
if(abs(log_lik_pre-log_lik)>10^(-5)){  
  
log_lik_pre=cost(sita,beta)
  
V=t(a)%*%beta  
 
P_k=exp(sita*V)/sum(exp(sita*V)) 

ave=a%*%P_k;var=((a^2)%*%P_k)-(a%*%P_k)^2

dsita=(cost(sita+h,beta)-cost(sita,beta))/h

dsita2=(cost(sita+h,beta)-2*cost(sita,beta)+cost(sita-h,beta))/(h^2)

db=c();db2=c()

vec=beta

for(j in 1:length(beta)){
  
vec_sub1=vec;vec_sub1[j]=vec_sub1[j]+h

vec_sub2=vec;vec_sub2[j]=vec_sub2[j]-h
  
db=c(db,(cost(sita,vec_sub1)-cost(sita,vec))/h)

db2=c(db2,(cost(sita,vec_sub1)-2*cost(sita,beta)+cost(sita,vec_sub2))/(h^2))
  
}
 
H=diag(c(dsita2,db2));G=diag(1/c(dsita2,db2));dg=c(dsita,db)

params=c(sita,beta)-h*G%*%dg

sita=params[1];beta=params[-1]

log_lik=cost(sita,beta)

}

print(cost(sita,beta))

}

```

```{r}

#p.79

library(dplyr)

rs=10;k=5

C_k_rs=array(sample(seq(0,1,0.01),k*rs,replace=T),dim=c(k,rs))

V_k_rs=array(sample(seq(0,1,0.01),k*rs,replace=T),dim=c(k,rs))


sita=10;E_hat=5

L=function(v,zf){
  
exp_mat=exp(-sita*v)

sum_exp_mat=apply(exp_mat,2,sum) 

p_k_rs=t(t(exp_mat)/c(sum_exp_mat))

q_rs=rep(1/rs,rs);f_k_rs=t(t(p_k_rs)*q_rs)

eta_rs=log(sum_exp_mat)-log(q_rs)-1

H_rs=-apply(p_k_rs,1,sum)*log(apply(p_k_rs,1,sum))

if(zf==0){
  
return(sita*(E_hat-sum(f_k_rs*C_k_rs))^2+sum(eta_rs*(q_rs-apply(f_k_rs,2,sum))^2))

}

if(zf==1){
  
z_f=(sum(q_rs*H_rs)-sum(q_rs*log(q_rs)))^2  

return(z_f+sita*(E_hat-sum(f_k_rs*C_k_rs))^2+sum(eta_rs*(q_rs-apply(f_k_rs,2,sum))^2))

}
  
if(zf==2){
  
return(sita*(E_hat-sum(f_k_rs*C_k_rs))^2)

} 

  
if(zf==3){

return(sum(f_k_rs*C_k_rs)-sum(f_k_rs*log(f_k_rs))/sita)  
  
}

}



op=3

ite=100000;h=0.01;eta=0.01

for(l in 1:ite){

exp_mat=exp(-sita*C_k_rs)

sum_exp_mat=apply(exp_mat,2,sum) 

p_k_rs=t(t(exp_mat)/c(sum_exp_mat))

q_rs=rep(1/rs,rs);f_k_rs=t(t(p_k_rs)*q_rs)

eta_rs=log(sum_exp_mat)-log(q_rs)-1

H_rs=-apply(p_k_rs,1,sum)*log(apply(p_k_rs,1,sum))
  
  
  
vec=C_k_rs  
    
for(i in 1:k){
for(j in 1:rs){  
  
vec_sub=vec;vec_sub[i,j]=vec_sub[i,j]+h

C_k_rs[i,j]=C_k_rs[i,j]-eta*(L(vec_sub,op)-L(vec,op))/h
  
}}  
  
print(L(C_k_rs,op))

#print(sum(f_k_rs*C_k_rs))
  
}


```

```{r}

#p.79 apply

library(dplyr)

rs=10;k=5

C_k_rs=array(sample(seq(0,1,0.01),k*rs,replace=T),dim=c(k,rs))

V_k_rs=array(sample(c(10:20)*10,k*rs,replace=T),dim=c(k,rs))


sita=1;E_hat=150

L=function(v,zf){
  
exp_mat=exp(-sita*v)

sum_exp_mat=apply(exp_mat,2,sum) 

p_k_rs=t(t(exp_mat)/c(sum_exp_mat))

q_rs=rep(1/rs,rs);f_k_rs=t(t(p_k_rs)*q_rs)

eta_rs=log(sum_exp_mat)-log(q_rs)-1

H_rs=-apply(p_k_rs,1,sum)*log(apply(p_k_rs,1,sum))

if(zf==0){
  
return(sita*(E_hat-sum(f_k_rs*V_k_rs))^2+sum(eta_rs*(q_rs-apply(f_k_rs,2,sum))^2))

}

if(zf==1){
  
z_f=(sum(q_rs*H_rs)-sum(q_rs*log(q_rs)))^2  

return(z_f+sita*(E_hat-sum(f_k_rs*V_k_rs))^2+sum(eta_rs*(q_rs-apply(f_k_rs,2,sum))^2))

}
  
if(zf==2){
  
return(sita*(E_hat-sum(f_k_rs*V_k_rs))^2)

} 

  
if(zf==3){

return(sum(f_k_rs*V_k_rs)+sum(f_k_rs*log(f_k_rs))/sita)  
  
  
}

}



op=3

ite=100000;h=0.01;eta=0.01

for(l in 1:ite){

exp_mat=exp(-sita*C_k_rs)

sum_exp_mat=apply(exp_mat,2,sum) 

p_k_rs=t(t(exp_mat)/c(sum_exp_mat))

q_rs=rep(1/rs,rs);f_k_rs=t(t(p_k_rs)*q_rs)

eta_rs=log(sum_exp_mat)-log(q_rs)-1

H_rs=-apply(p_k_rs,1,sum)*log(apply(p_k_rs,1,sum))
  
  
  
vec=C_k_rs  
    
for(i in 1:k){
for(j in 1:rs){  
  
vec_sub=vec;vec_sub[i,j]=vec_sub[i,j]+h

C_k_rs[i,j]=C_k_rs[i,j]-eta*(L(vec_sub,op)-L(vec,op))/h
  
}}  
  
print(L(C_k_rs,op))

#print(sum(f_k_rs*log(f_k_rs)))
  
}


```


```{r}

#p.197

library(dplyr)

s=5;r=10;a_s=5

q_rs=array(0.001,dim=c(r,s));v_rs=matrix(sample(1:10,r*s,replace = T),nrow=r)/1000

x_a=rep(0.1,a_s);y_a=sample(1:10,a_s,replace = T)

beta=0.01





ite=100000;eta=0.001

for(k in 1:ite){
  
cost=function(alpha){
  
return(abs(sum(pnorm(x_a+alpha*(y_a-x_a)))-sum(pnorm(q_rs+alpha*(v_rs-q_rs)))))  
  
}  

for(l in 1:ite){
  
bs=beta-eta*((cost(beta+h)-cost(beta))/h)  
  
if(bs>=0 & bs<=1){  
  
beta=beta-eta*((cost(beta+h)-cost(beta))/h)

}

}

x_a=x_a+beta*(y_a-x_a)

q_rs=q_rs+beta*(v_rs-q_rs)

print(cost(beta))

if(beta<10^(-10)){
  
print("stop")  
  
}

}

```

```{r}

#p.254

n=10;a=10

x_a=sample(1:10,a,replace=T)

p_a_rs=array(0,dim=c(1,n+1))

for(j in 1:a){
  
p_a_rs_sub=cbind(rep(j,n),array(sample(seq(0.01,1,0.01),n^2,replace=T),dim=c(n,n)) ) 
  
p_a_rs=rbind(p_a_rs,p_a_rs_sub) 
 
}

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

p_a_rs=data.frame(p_a_rs) %>% rename(a=X1)



cost=function(lambda,Q_rs_mat){
  
Q_rs_mat=abs(Q_rs_mat);lambda=abs(lambda)  
  
t=sum(Q_rs_mat)

dlam=c()

for(j in 1:a){
  
p_a_rs_sub=p_a_rs %>% filter(a==j) %>% select(-a)  
  
dlam=c(dlam,(sum(Q_rs_mat*p_a_rs_sub)-x_a[j])^2)  
  
}

L=t*log(t)-t-sum(Q_rs_mat*log(Q_rs_mat))+sum(Q_rs_mat*log(Q_rs_mat/sum(Q_rs_mat)))+sum((lambda^2)*dlam)

return(L)
  
}



ite=100

h=0.001

lam=rep(1,a);Q_rs=array(1,dim=c(n,n))

for(l in 1:ite){
  
vec=Q_rs

dQ=array(0,dim=c(n,n));ddQ=array(0,dim=c(n,n))

for(i in 1:n){
for(j in 1:n){  
  
vec_sub1=vec;vec_sub1[i,j]=vec_sub1[i,j]+h 

vec_sub2=vec;vec_sub2[i,j]=vec_sub2[i,j]-h
  
dQ[i,j]=(cost(lam,vec_sub1)-cost(lam,vec))/h

ddQ[i,j]=(cost(lam,vec_sub1)-2*cost(lam,vec)+cost(lam,vec_sub2))/(h^2)
  
}}

Q_rs=Q_rs-h*dQ/ddQ



vec=lam

for(j in 1:a){
  
vec_sub1=vec;vec_sub1[j]=vec_sub1[j]+h

vec_sub2=vec;vec_sub2[j]=vec_sub2[j]-h

dlam=(cost(vec_sub1,Q_rs)-cost(vec,Q_rs))/h

ddlam=(cost(vec_sub1,Q_rs)-2*cost(vec,Q_rs)+cost(vec_sub2,Q_rs))/(h^2)

lam[j]=vec[j]-h*dlam/ddlam
  
}

print(cost(lam,Q_rs))
  
}

```

```{r}


#p.256

m=10;n=20;K=10

U=array(sample(seq(0.01,1,0.01),m*n,replace=T),dim=c(m,n))

p=sample(seq(0.01,1,0.01),n)

x_a=sample(1:10,m,replace = T)

cost=function(o,lambda,vs,t_vec){
  
o=abs(o);lambda=abs(lambda);  
  
return(K*sum((t_vec*p-o)^2)+(vs^2)*(t_vec-sum(o))^2+sum((lambda^2)*(U%*%o-x_a)^2))
  
}


O_s=rep(1,n)

lam=rep(1,m)

v=1;t=10


ite=100000

h=0.01

cost_pre_val=1;cost_val=0.5

for(l in 1:ite){
  
if(abs(cost_pre_val-cost_val)>10^(-8)){  
  
cost_pre_val=cost_val  
  
vec=O_s

do_s=c();ddo_s=c()

for(j in 1:length(O_s)){
  
vec_sub1=vec;vec_sub1[j]=vec_sub1[j]+h

vec_sub2=vec;vec_sub2[j]=vec_sub2[j]-h

do_s=c(do_s,(cost(vec_sub1,lam,v,t)-cost(vec,lam,v,t))/h)

ddo_s=c(ddo_s,(cost(vec_sub1,lam,v,t)-2*cost(vec,lam,v,t)+cost(vec_sub2,lam,v,t))/(h^2))
  
}

O_s=O_s-h*do_s/ddo_s;O_s=abs(O_s)

vec=lam

dlam=c();ddlam=c()

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

vec_sub1=vec;vec_sub1[j]=vec_sub1[j]+h

vec_sub2=vec;vec_sub2[j]=vec_sub2[j]-h
    
dlam=c(dlam,(cost(O_s,vec_sub1,v,t)-cost(O_s,vec,v,t))/h)

ddlam=c(ddlam,(cost(O_s,vec_sub1,v,t)-2*cost(O_s,vec,v,t)+cost(O_s,vec_sub2,v,t))/(h^2))
  
}

lam=lam-h*dlam/ddlam

dv=(cost(O_s,lam,v+h,t)-cost(O_s,lam,v,t))/h

ddv=(cost(O_s,lam,v+h,t)-2*cost(O_s,lam,v,t)+cost(O_s,lam,v-h,t))/(h^2)

v=v-h*dv/ddv

dt=(cost(O_s,lam,v,t+h)-cost(O_s,lam,v,t))/h

ddt=(cost(O_s,lam,v,t+h)-2*cost(O_s,lam,v,t)+cost(O_s,lam,v,t-h))/(h^2)

t=t-h*dt/ddt

cost_val=cost(O_s,lam,v,t)

}

print(cost(O_s,lam,v,t))
  
}



```