スポーツの数理科学　－もっと楽しむための数字の読み方ー　竹内啓・藤野和建先生

```{r}

#スポーツの数理科学

library(dplyr)

data=data.frame(anscombe)

a=c();b=c();R2=c();ht=c();sigma=c()

u_vec=array(0,dim=c(4,nrow(data)))

for(j in 1:4){
  
X=data[,colnames(data) %in% paste0("x",j)]  

Y=data[,colnames(data) %in% paste0("y",j)]

tre=cov(X,Y)/var(X);inter=mean(Y)-tre*mean(X)

a=c(a,cov(X,Y)/cov(X,X))

b=c(b,mean(Y)-tre*mean(X))

Se=sum((Y-(inter+X*tre))^2)

SR=sum(((inter+tre*X)-mean(Y))^2)

S_yy=Se+SR

ht=1/length(Y)+(X-mean(X))^2/sum((X-mean(X))^2)

et=(Y-(inter+tre*X))^2

ht[ht==1]=0.99999

sigma=sqrt((Se-et*(1/(1-ht)))/(length(Y)-3))

u=et/(sigma*sqrt(1-ht))

u_vec[j,]=u

R2=c(R2,SR/S_yy)
  
}

#uは標準正規分布にしたがう

#データ直線に当てはめられるものをさがす

u_vec>qnorm(1-0.01)


```

```{r}

#p.19

data=data.frame(anscombe)



for(j in 1:4){
  
X=data[,colnames(data) %in% paste0("x",j)]

Y=data[,colnames(data) %in% paste0("y",j)]

#iterarions

m=0.1;A=0.1;B=0.1;ave=1

Q=10000

eta=0.0001

diff=1

while(diff>0.0001)({
  
Q_pre=Q  
  
m=m-eta*(-2*sum(Y-A-B*exp(-m*X)))  
  
A=A-eta*(-2*sum(Y-A-B*exp(-m*X)*exp(-m*X)))

B=B-eta*(-2*B*sum((Y-A-B*exp(-m*X))*X*exp(-m*X)))

Q=sum((Y-A-B*exp(-m*X))^2)

diff=abs(Q_pre-Q)

print(Q)

})
  
pre=A+B*exp(-m*X)

Se=sum((Y-pre)^2)

sigma=sqrt(Se/length(Y))

}

```

```{r}

#p.22

data=data.frame(anscombe)

j=1

X=data[,colnames(data) %in% paste0("x",j)]

X=ifelse(X-mean(X)>0,1,0)

t=c(1:length(X))

Y=data[,colnames(data) %in% paste0("y",j)]

exp_sig=1

if(exp_sig==0){

L=function(a,b,sig,ave){
  
u=(Y-a-b*(t-ave))/sig

return(sum(X*(log(dnorm(u))-log(sig)))+sum((1-X)*log(1-pnorm(u))))
  
}

}else{
  
L=function(a,b,sig,ave,m_val){
  
u=(Y-a-b*exp(-m_val*(t-ave)))/sig

return(sum(X*(log(dnorm(u))-log(sig)))+sum((1-X)*log(2-pnorm(u))))
  
}  
  
}


ite=10000

alpha=0.01;beta=0.01;sigma=1;average=1;m=1

eta=0.0001

h=0.01

if(exp_sig==0){

diff=1

while(diff>0.00001)({
  
log_lik_pre=L(alpha,beta,sigma,average)  
  
alpha=alpha+eta*(L(alpha+h,beta,sigma,average)-L(alpha,beta,sigma,average))/h

beta=beta+eta*(L(alpha,beta+h,sigma,average)-L(alpha,beta,sigma,average))/h

sigma=sigma+eta*(L(alpha,beta,sigma+h,average)-L(alpha,beta,sigma,average))/h

average=average+eta*(L(alpha,beta,sigma,average+h)-L(alpha,beta,sigma,average))/h
  
log_lik=L(alpha,beta,sigma,average)

#print(log_lik)

diff=abs(log_lik_pre-log_lik)

})

}else{
  
diff=1

while(diff>0.00001)({
  
log_lik_pre=L(alpha,beta,sigma,average,m)  
  
alpha=alpha+eta*(L(alpha+h,beta,sigma,average,m)-L(alpha,beta,sigma,average,m))/h

beta=beta+eta*(L(alpha,beta+h,sigma,average,m)-L(alpha,beta,sigma,average,m))/h

sigma=sigma+eta*(L(alpha,beta,sigma+h,average,m)-L(alpha,beta,sigma,average,m))/h

average=average+eta*(L(alpha,beta,sigma,average+h,m)-L(alpha,beta,sigma,average,m))/h
 
if(m+eta*(L(alpha,beta,sigma,average,m+h)-L(alpha,beta,sigma,average,m))/h>0){

m=m+eta*(L(alpha,beta,sigma,average,m+h)-L(alpha,beta,sigma,average,m))/h
 
}

log_lik=L(alpha,beta,sigma,average,m)

#print(log_lik)

diff=abs(log_lik_pre-log_lik)

})
  
}


```

```{r}

#Bradley-Terry model 

#サンプルデータの作成（野球の試合のデータ）
dat1=array(0,dim=c(6,6))

dat1[1,]=c(0,14,12,18,15,14)
dat1[2,]=c(11,0,13,18,17,16)
dat1[3,]=c(8,11,0,9,16,16)
dat1[4,]=c(7,7,16,0,11,15)
dat1[5,]=c(8,7,10,13,0,16)
dat1[6,]=c(12,9,9,11,8,0)

rownames(dat1)=c("広島","巨人","阪神","大洋","中日","ヤクルト")

colnames(dat1)=c("広島","巨人","阪神","大洋","中日","ヤクルト")


#総試合数のデータ
nij=array(0,dim=c(6,6))

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

nij[i,j]=dat1[i,j]+dat1[j,i]    

}
}

#勝数
t=apply(dat1,1,sum)

#各チームの強さの平均が50になるようにする
k=300

#パラメータの初期値
pi=k*rep(1/6,6)

#反復計算回数
ite=100

for(l in 1:ite){

#以前のパラメータを保存
pi_pre=pi  

r=c()

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

pi_sub=pi[-i]

ri=sum(nij[i,-i]/(pi[i]+pi_sub))

r=c(r,ri)

}

pi_sub=t/r

#パラメータの更新
pi=k*pi_sub/sum(pi_sub)

pi_mat=array(0,dim=c(length(pi),length(pi)))

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

pi_mat[i,]=pi[i]+pi  

}

diag(pi_mat)=1

l_sub=as.matrix(nij)*log(pi_mat)

l_sub=sum(l_sub)/2

#対数尤度
l=sum(t*log(pi))-l_sub

print(l)

}

pij=array(0,dim=c(length(pi),length(pi)))

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

pij[j,-j]=pi[j]/(pi[j]+pi[-j])  

}

#BTモデルの検定（尤度比検定）

l0=dat1*log(dat1);diag(l0)=0

l0=sum(l0)-sum((nij*log(nij))[upper.tri(nij*log(nij))])

#H:モデルの精度の検定-BTモデルがなりたっているか

kai=2*(l0-l)

v=(6-1)*(6-2)/2

qchisq(1-0.05,df=v)

#H:pi=一定

l2=-sum((nij[upper.tri(nij)])*log(2))

v=(6-1)

kai2=2*(l-l2)

qchisq(1-0.05,v)

#表2.4 p.42

aij=array(0,dim=c(6,6))

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

aij[i,j]=log2(((dat1[i,j]+0.5)*pi[j])/((dat1[j,i]+0.5)*pi[i]))  

}  
}

mu2=eigen(t(aij)%*%aij)$values

v=eigen(t(aij)%*%aij)$vectors

u=t(t(aij%*%v)/sqrt(mu2))


```

```{r}

#Bradley-Terry model p.27

dat1=array(0,dim=c(6,6))

dat1[1,]=c(0,14,12,18,15,14)
dat1[2,]=c(11,0,13,18,17,16)
dat1[3,]=c(8,11,0,9,16,16)
dat1[4,]=c(7,7,16,0,11,15)
dat1[5,]=c(8,7,10,13,0,16)
dat1[6,]=c(12,9,9,11,8,0)

nij=array(0,dim=c(6,6))

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

nij[i,j]=dat1[i,j]+dat1[j,i]    

}
}

t=apply(dat1,1,sum)


k=300


#levenberg(1944)

#関数定義
f=function(x){

pi=x

f1=c()

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

f1=c(f1,t[i]/pi[i]-sum(nij[i,-i]/(pi[i]+pi[-i])))

}

return(c(f1))

}

#初期値
X=k*rep(1/6,6)

#反復計算回数
ite=100000

#学習率
eta=10^(-2)

#差分メッシュ
h=0.01

#正の小さいスカラ
w=0.01

for(l in 1:ite){

#ヤコビアンを計算する箱を用意
df=array(0,dim=c(length(X),length(X)))

#1次差分を用いて数値微分し、ヤコビアンを計算する
for(j in 1:length(X)){

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

df[j,]=(f(vec)-f(X))/h  

}

#逆行列が特異にならないよう、対角成分をwで調整
A=df%*%t(df)+diag(w,length(X))  

y=df%*%f(X)

#座標点更新の際に必要なベクトルを計算（逆行列に関しては特異値分解で計算)
dx=svd(A)$u%*%diag(1/svd(A)$d)%*%t(svd(A)$v)%*%(-y)

#学習幅に従って座標点を更新
X=X+eta*dx

X=k*X/sum(X)

print(f(X))

}


```

```{r}

#Bradley-Terry model p.27

dat1=array(0,dim=c(6,6))

dat1[1,]=c(0,14,12,18,15,14)
dat1[2,]=c(11,0,13,18,17,16)
dat1[3,]=c(8,11,0,9,16,16)
dat1[4,]=c(7,7,16,0,11,15)
dat1[5,]=c(8,7,10,13,0,16)
dat1[6,]=c(12,9,9,11,8,0)

nij=array(0,dim=c(6,6))

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

nij[i,j]=dat1[i,j]+dat1[j,i]    

}
}

t=apply(dat1,1,sum)


k=300

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

#関数定義
f=function(x){

pi=x

f1=c()

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

f1=c(f1,t[i]/pi[i]-sum(nij[i,-i]/(pi[i]+pi[-i])))

}

return(c(f1))

}

#初期値
X=k*rep(1/6,6)

#ヘッシアンの近似行列の初期値
H=diag(1,length(X))

#反復計算回数
ite=10000

#最小点から非常にかけ離れた点になるのを防ぐためのスカラ
eta=10^(-2)

for(l in 1:ite){

#各代数方程式の値が入ったベクトル
f_val=f(X)

#座標点に対する更新量ベクトル
sigma=t(t(-H%*%f_val))

#etaによって学習幅を設定し座標点を更新
X=X+eta*sigma

X=k*X/sum(X) 

#以前の座標点と更新された座標点の関数値の差のベクトル(y)
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}


#Bradley-Terry model p.27

dat1=array(0,dim=c(6,6))

dat1[1,]=c(0,14,12,18,15,14)
dat1[2,]=c(11,0,13,18,17,16)
dat1[3,]=c(8,11,0,9,16,16)
dat1[4,]=c(7,7,16,0,11,15)
dat1[5,]=c(8,7,10,13,0,16)
dat1[6,]=c(12,9,9,11,8,0)

nij=array(0,dim=c(6,6))

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

nij[i,j]=dat1[i,j]+dat1[j,i]    

}
}

t=apply(dat1,1,sum)


k=10



#Broyden Quasi-Newton Algorithm p.213

f=function(x){

pi=x

f1=c()

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

f1=c(f1,t[i]/pi[i]-sum(nij[i,-i]/(pi[i]+pi[-i])))

}

return(c(f1))

}

X=k*rep(1/6,6)

ite=10000

eta=0.001

H=diag(f(X))

for(l in 1:ite){

if(sum(abs(f(X)))>10^(-9)){  

X_pre=X  

X=X-eta*H%*%f(X) 

X=k*X/sum(X) 

s=X-X_pre

y=f(X)-f(X_pre)

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

print(sum(abs(f(X))))

}

}  



```




```{r}

#p.40

#table 2.2

nij=array(0,dim=c(3,3))

nij[1,]=c(0,11,17)

nij[2,]=c(15,0,10)

nij[3,]=c(6,17,0)

X=nij

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

nij[i,j]=X[i,j]+X[j,i]  
    
}
}

t=c(28,25,20)

```

```{r}

#table 2.3

library(dplyr)

pr12_ratio=function(x){
  
return((26-x)*(27-x)*(29-x)/(x*(x-1)*(x-5)))  
  
}

const=c()

for(j in 6:25){
  
const=c(const,pr12_ratio(j))  
  
}

const=cumprod(const)

pr_data=data.frame(x12=c(6:25),const=const,pr=(const)/sum((const)))

#勝数が11以下になる確率

sum(pr_data$pr[pr_data$x12<=11])

#予想される阪神の広島に対する勝数

pre_x=pr_data$x12[pr_data$pr==max(pr_data$pr)]

X=array(0,dim=c(3,3))

X[1,]=c(0,pre_x,28-pre_x)

X[2,]=c(26-pre_x,0,pre_x-1)

X[3,]=c(pre_x-5,25-pre_x,0)


```

```{r}



#Bradley-Terry model p.148

dat1=array(0,dim=c(6,6))

dat1[1,]=c(0,9,8,11,9,8)
dat1[2,]=c(4,0,9,10,9,7)
dat1[3,]=c(3,4,0,5,7,9)
dat1[4,]=c(2,3,8,0,7,9)
dat1[5,]=c(4,3,6,6,0,8)
dat1[6,]=c(5,5,4,4,4,0)

dat2=array(0,dim=c(6,6))

dat2[1,]=c(0,7,5,5,4,7)
dat2[2,]=c(5,0,7,4,4,4)
dat2[3,]=c(4,4,0,8,4,5)
dat2[4,]=c(7,8,4,0,7,7)
dat2[5,]=c(6,8,9,4,0,4)
dat2[6,]=c(6,9,7,6,8,0)

Z=dat1+dat2

t=apply(Z,1,sum)

w=c()

for(i in 1:6){
  
w=c(w,sum(c(dat1[i,],dat1[,i]))-dat1[i,i])  
  
}

k=1

omega=rep(0.1,6)

pi=rep(1/6,6)*k

ite=100

for(l in 1:ite){
  
pi_pre=pi  

r=c();vec=c()

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

 

for(i in 1:6){
  
val=c();val2=c();   

for(j in 1:6){
  
if(j!=i){  
  
val=c(val,Z[i,j]*omega[i]*omega[j]/(omega[i]*omega[j]*pi[i]+pi[j])+Z[j,i]/(omega[i]*omega[j]*pi[j]+pi[i])) 

val2=c(val2,Z[i,j]*pi[i]*omega[j]/(omega[i]*omega[j]*pi[i]+pi[j])+Z[j,i]*pi[j]*omega[j]/(omega[i]*pi[j]*pi[j]+pi[i]))
  
}
  
mat[i,j]=Z[i,j]*log(omega[i]*omega[j]*pi[i]+pi[j])  
  
}  
  
r=c(r,sum(val));vec=c(vec,sum(val2))

}
    
omega=w/vec

pi_hat=t/r 

pi=k*pi_hat/sum(pi_hat)

l0=sum(w*log(omega))+sum(t*log(pi))-sum(mat)

#print(sum(abs(pi-pi_pre)))

print(l0)
  
}


```