質的情報の多変量解析　松田紀之先生　朝倉出版

```r

#p.45

#表2.1
n=c(17,13,19,11,18,12)
m11=rep(15,6)

G2=2*sum(n*log(n/m11))
df=length(n)-1
1-pchisq(G2,df)

#M1-1:AIC
-2*log(dmultinom(n,sum(n),prob=m11/sum(n)))+2*1
#M1-2:AIC
-2*log(dmultinom(n,sum(n),prob=n/sum(n)))+2*(df+1)


#表2.2
n=c(182,174,155,165)
m21=rep(169,4)

G2=2*sum(n*log(n/m21))
df=length(n)-1
1-pchisq(G2,df)

#M2-1:AIC
-2*log(dmultinom(n,sum(n),prob=m21/sum(n)))+2*1
#M2-2:AIC
-2*log(dmultinom(n,sum(n),prob=n/sum(n)))+2*(df+1)

```

```r
#p.66 質的情報の多変量解析
#表2.1の例

m12=c(17,13,19,11,18,12)
X=rbind(rep(-1,5),diag(1,5))
X=cbind(rep(1,nrow(X)),X)
n=c(17,13,19,11,18,12)
n=m12

ite=10^3
eta=0.01
params=rep(0.1,ncol(X))
log_lik=-100
log_lik_vec=c()

for(l in 1:ite){
log_lik_pre=log_lik
dL=apply(X*c(n)-X*c(exp(X%*%params)),2,sum)
dLL=-t(X)%*%diag(c(exp(X%*%params)))%*%X
params=params-eta*(solve(dLL)%*%dL)
log_lik=sum((X%*%params)*n-exp(X%*%params))
print(log_lik)
log_lik_vec=c(log_lik_vec,log_lik)
}
plot(log_lik_vec)

u1=-sum(params[-1])
u123456=c(u1,params[-1])
ubar=params[1]

#全平均の値とparams[1]の値が等しい
log(prod(exp(X%*%params)))/length(u123456)

#パラメータの分散
sw=sqrt(diag(-solve(dLL)))

sd_u123456=c(sqrt(sum(-solve(dLL)[-1,-1])),sw[-1])
sd_ubar=sw[1]

z_value=u123456/sd_u123456

#上側信頼区間値
u123456+qnorm(1-0.05/2)*sd_u123456

#下側信頼区間値
u123456-qnorm(1-0.05/2)*sd_u123456


```


```r


#p.69

test_static=function(nij,alpha){
n_ip=apply(nij,1,sum)
n_pj=apply(nij,2,sum)
n_pp=sum(nij)
mij_hat=t(t(n_ip))%*%t(n_pj)/n_pp
phat_ip=n_ip/n_pp
phat_pj=n_pj/n_pp
X2=sum((nij-(mij_hat))^2/(mij_hat))
I=nrow(nij)
J=ncol(nij)
return(c(X2,paste0("棄却値：",qchisq(1-alpha,df=(I-1)*(J-1)))))
}

#表3.1
nij1=matrix(c(52,61,73,921,902,851),nrow=3)

#表3.2
nij2=matrix(c(5,7,8,9,3,21,7,9,6),nrow=3)

#表3.3
nij3=matrix(c(181,160,35,56),nrow=2)

#表3.4
nij4=matrix(c(50,28,11,14,0,45,174,78,150,42,8,84,110,185,72,18,154,223,714,320,8,55,96,447,411),nrow=5)

test_static(nij1,0.05)
test_static(nij2,0.05)
test_static(nij3,0.05)
test_static(nij4,0.05)

```

```r

#p.83

mhat=function(nij,l){

ubar=mean(log(nij))
uA=apply(log(nij),1,mean)-ubar
uB=apply(log(nij),2,mean)-ubar
uAB=array(0,dim=dim(nij))

for(i in 1:nrow(uAB)){
for(j in 1:ncol(uAB)){
uAB[i,j]=log(nij)[i,j]-ubar-(uA[i]+uB[j])  
}}
if(l==0){
mij_hat=array(ubar,dim=c(nrow(nij),ncol(nij)))
for(i in 1:nrow(nij)){
for(j in 1:ncol(nij)){
mij_hat[i,j]=mij_hat[i,j]+(uA[i]+uB[j])  
}}
}
  
if(l==1){
mij_hat=matrix(rep(ubar+uA,length(uB)),nrow=length(uA))  
}
  
if(l==2){
mij_hat=t(matrix(rep(uB+ubar,length(uA)),nrow=length(uB)))  
}
  
if(!(l %in% c(0,1,2))){
mij_hat=array(ubar,dim=c(nrow(nij),ncol(nij)))
for(i in 1:nrow(nij)){
for(j in 1:ncol(nij)){
mij_hat[i,j]=mij_hat[i,j]+(uA[i]+uB[j])+uAB[i,j]  
}}
}
  
return(mij_hat)
}

AIC=function(m,n,p){
return(c(-2*log(dmultinom(c(n),sum(n),prob=c(m)/sum(n)))+2*(p)))
}


#表3.1
nij1=matrix(c(52,61,73,921,902,851),nrow=3)
I=nrow(nij1)
J=ncol(nij1)

#[AB]
mijABhat=exp(mhat(nij1,3))
AIC(mijABhat,nij1,(I*J))

#[A][B]
mijA_Bhat=exp(mhat(nij1,0))
AIC(mijA_Bhat,nij1,(I+J-1))
X2=sum((nij1-mijA_Bhat)^2/mijA_Bhat)
G2=2*sum(nij1*log(nij1/mijA_Bhat))


#[A]
mijAhat=exp(mhat(nij1,1)) 
AIC(mijAhat,nij1,(I-1+1))
X2=sum((nij1-mijAhat)^2/mijAhat)
G2=2*sum(nij1*log(nij1/mijAhat))


#[B]
mijBhat=exp(mhat(nij1,2)) 
AIC(mijBhat,nij1,(J-1+1))
X2=sum((nij1-mijBhat)^2/mijBhat)
G2=2*sum(nij1*log(nij1/mijBhat))


```

```r

#p.149

library(dplyr)

mat=matrix(c(50,28,11,14,0.5,45,174,78,150,42,8,84,110,185,72,18,154,223,714,320,8,55,96,447,411),ncol=5)

a=nrow(mat);b=ncol(mat)



log_lik=function(p,u_ab){

u=p[1];p=p[-1];u_a=p[1:a];p=p[-c(1:a)];u_b=p[1:b];p=p[-c(1:b)];

m_mat=array(0,dim=c(a,b))

n_mat=0

for(i in 1:a){
for(j in 1:b){

m_mat[i,j]=u+u_a[i]+u_b[j]+u_ab[i,j]

n_mat=n_mat+sum(log(c(1:mat[i,j])))

}}

val=sum(mat*log(m_mat)-m_mat)-n_mat

return(val)

}

param=rep(0.1,a+b+1)

uab=array(0.1,dim=c(a,b))


#iteration

ite=100000

eta=0.01

h=0.01

for(l in 1:ite){

vec=param

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

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

param[j]=param[j]+eta*(log_lik(vec_sub,uab)-log_lik(vec,uab))/h

}

vec=uab

for(i in 1:a){
for(j in 1:b){

vec_sub=vec;vec_sub[i,j]=vec_sub[i,j]+h

uab[i,j]=uab[i,j]+eta*(log_lik(param,vec_sub)-log_lik(param,vec))/h

}}

print(log_lik(param,uab))

}






```


```r

#p.173



mat=matrix(c(118,58,20,16,74,34,14,20,20,14,30,80,18,22,66,120),ncol=4)

a=nrow(mat);b=ncol(mat);t=2

pij=mat/sum(mat)

pi_A_B1=pij

mat=matrix(c(18,78,20,36,74,84,14,20,50,14,30,80,18,22,66,120),ncol=4)

a=nrow(mat);b=ncol(mat);t=2

pij=mat/sum(mat)

pi_A_B2=pij


pi_x=sample(1:9,t,replace=T);pi_x=pi_x/sum(pi_x)

pi_A=matrix(sample(1:9,a*t,replace=T),nrow=a,ncol=t)

pi_A=t(t(pi_A)/apply(pi_A,2,sum))

pi_B=matrix(sample(1:9,b*t,replace =T),nrow=b,ncol=t)

pi_B=t(t(pi_B)/apply(pi_B,2,sum))



ite=210

for(l in 1:ite){

pi_A_B_sum=array(0,dim=c(a,b))

for(s in 1:t){

pi_Xt=pi_x[s]

pi_A_Xt=pi_A[,s]

pi_B_Xt=pi_B[,s]

pi_A_B_sum=pi_A_B_sum+pi_Xt*t(t(c(pi_A_Xt)))%*%t(c(pi_B_Xt))  

}

df_vec=c()  

dp_vec=c()

for(s in 1:t){ 

if(s==1){

pi_A_B=pi_A_B1

}else{

pi_A_B=pi_A_B2

}

pi_Xt=pi_x[s]

pi_A_Xt=pi_A[,s]

pi_B_Xt=pi_B[,s]

pi_Xt_pre=pi_Xt

pi_A_Xt_pre=pi_A_Xt

pi_B_Xt_pre=pi_B_Xt

pi_A_B_Xt=pi_Xt*t(t(c(pi_A_Xt)))%*%t(c(pi_B_Xt))

pi_Xt_A_B=pi_A_B_Xt/pi_A_B_sum

pi_Xt=sum(pi_A_B*pi_Xt_A_B)

pi_A_Xt=apply(pi_A_B*pi_Xt_A_B,1,sum)/pi_x[s]

pi_B_Xt=apply(pi_A_B*pi_Xt_A_B,2,sum)/pi_x[s]

df=sum(abs(pi_Xt-pi_Xt_pre))+sum(abs(pi_A_Xt-pi_A_Xt_pre))+sum(abs(pi_B_Xt-pi_B_Xt_pre))

dp=sum(abs(pi_A_B-pi_A_B_Xt))

pi_x[s]=pi_Xt

pi_A[,s]=pi_A_Xt

pi_B[,s]=pi_B_Xt

df_vec=c(df_vec,df)

dp_vec=c(dp_vec,dp)

}

#print(sum(df_vec))

print(sum(dp_vec))

}

```