カテゴリカルデータのモデル分析　坂元 慶行先生　共立出版




```r

#p.78

#年齢別データ
data=data.frame(age=c(20:94),ans1=c(3,5,10,5,6,6,5,10,9,6,12,10,11,14,12,6,6,7,14,13,7,6,5,9,10,7,13,11,8,7,4,12,9,6,7,6,5,5,9,5,4,7,5,5,7,1,1,3,0,1,0,3,2,3,0,1,1,1,2,0,2,0,1,rep(0,12)),ans2=c(4,5,3,6,4,2,7,3,5,6,4,2,6,5,7,4,1,5,3,3,4,3,2,3,3,5,1,3,0,1,5,3,7,4,5,2,5,4,3,3,2,3,5,1,2,2,3,1,3,4,1,6,3,1,0,2,0,1,2,0,1,0,1,0,0,1,rep(0,8),1))

#20代、30代、として計算
nj1=c()
nj2=c()
ages=sort(unique(round(data$age/10)*10))
C=length(ages)
for(j in 1:C){
nj1=c(nj1,sum(data$ans1[data$age>=ages[j] & data$age<=(ages[j]+9)]))
nj2=c(nj2,sum(data$ans2[data$age>=ages[j] & data$age<=(ages[j]+9)]))
}
nj=nj1+nj2

D=diag(1,C)
for(j in 1:C){
for(i in j:(j+2)){  
if(i==j){D[i,j]=1}
if((i==(j+1))&(i<=C)){D[i,j]=-2} 
if((i==(j+2))&(i<=C)){D[i,j]=1}
}}

#対数尤度の最大化
#newton-raphson
qj=rep(1,C)
ite=10

for(l in 1:ite){
pj=exp(qj)/(1+exp(qj))  
H=diag(-nj*pj*(1-pj))
df=nj1-nj*pj
qj=c(qj-solve(H)%*%df)
loglik=sum(nj1*qj-nj*log(1+exp(qj)))
print(loglik)
}



#ABIC
ABIC=function(q){
v=q[1]
qvec=q[-1]
q01=qvec[-c(1:2)]
val=sum(nj1*q01-nj*log(1+exp(q01)))
for(j in 1:C){
val=val-sum(qvec[(j+2)]-2*qvec[(j+1)]+qvec[j])^2/(2*v^2)  
}
val=-2*val-C*log(2*pi)
p=exp(q01)/(1+exp(q01))
L=diag(sqrt(nj*p*(1-p)))
val=val+log(det(t(L)%*%L+t(D)%*%D/(v^2)))
return(val)
}

param=rep(1,(C+3))
ite=10000
eta=0.01
h=0.01

for(l in 1:ite){

dparam=c()
for(j in 1:length(param)){
param_sub=param
param_sub[j]=param_sub[j]+h
dparam=c(dparam,(ABIC(param_sub)-ABIC(param))/h)
}  
param=param-eta*dparam
print(ABIC(param))
}

exp(param[-1])/(1+exp(param[-1]))

```

```r

#カテゴリかるデータのモデル分析　p.106

#ベイズ型一次元二値回帰モデル

library(dplyr)

data=data.frame(I0_1=c(2,3,2,1,3,3,0,2,1,2,3,1,2,1,0,2,2,2,4,0,1,1,1,0,0,1,2,1,1,1,0,1,1,0,1,1,0,2,0,0,2,1,2,0,1,1,2,2,2,1,0,0,0,0,0,1,3,0,1,2,0,0,1,1,0,0,3,1,0,1,2,0,2,0,3,2,3,1,1,3,2,1,0,2,3,1,3,1,2,1,2,2,2,2,3,3,2,0,2,1,1),I0_2=c(2,1,2,3,1,1,4,2,3,2,1,3,2,3,4,2,2,2,0,4,3,3,3,4,4,3,2,3,3,3,4,3,3,4,3,3,4,2,4,4,2,3,2,4,3,3,2,2,2,3,4,4,4,4,4,3,1,4,3,2,4,4,3,3,4,4,1,3,4,3,2,4,2,4,1,2,1,3,3,1,2,3,4,2,1,3,1,3,2,3,2,2,2,2,1,1,2,4,2,3,3))



summary1=data %>% group_by(I0_1) %>% summarise(n1=n())

summary2=data %>% group_by(I0_2) %>% summarise(n2=n())

summary=cbind(summary1,summary2)

n1_j=summary1$n1;n2_j=summary2$n2;n_j=summary1$n1+summary2$n2

class=sort(unique(c(data$I0_1,data$I0_2)))


D=array(0,dim=c(length(class),length(class)))

for(j in 1:(length(class)-2)){
  
D[j:(j+2),j]=c(1,-2,1)  
  
}

D[(length(class)-1),(length(class)-1)]=1

D[length(class),(length(class)-2):length(class)]=c(1,-2,1)

ABIC=function(q,q0,q1,v){
  
log_liklihood=sum(n1_j*q-n_j*log(1+exp(q)))  
  
val=c(length(class)*log(1/sqrt(2*pi*v^2)),-(q[1]-2*q0+q1)^2/(2*v^2),-(q[2]-2*q[1]+q0)^2/(2*v^2))

for(s in 3:length(class)){
  
val=c(val,-(q[j]-2*q[j-1]+q[j-2])^2/(2*v^2))  
  
}

p_j=exp(q)/(1+exp(q))  

L=diag(sqrt(c(n_j*p_j*(1-p_j))))

return(-2*(log_liklihood+(sum(val)))-length(class)*log(2*pi)+log(det(L%*%L+t(D)%*%D/(v^2))))

}


log_likli=function(a,b){
  
q=a+b*class

log_like=n1_j*q-n_j*log(1+exp(q))

log_like=sum(log_like)

return(log_like)
  
}


q=rep(1,length(class));

a=rep(1,length(class));b=rep(1,length(class))

ite=10000;h=0.01;eta=0.0001


for(l in 1:ite){
  
q=a+b*class  
  
p_j=exp(q)/(1+exp(q))  

vec=a

for(i in 1:length(a)){
  
vec_sub=vec;vec_sub[i]=vec_sub[i]+h  
  
a[i]=a[i]+eta*(log_likli(vec_sub,b)-log_likli(vec,b))/h

}

vec=b

for(i in 1:length(b)){
  
vec_sub=vec;vec_sub[i]=vec_sub[i]+h  
  
b[i]=b[i]+eta*(log_likli(a,vec_sub)-log_likli(a,vec))/h

}

print(log_likli(a,b))

}



ite=10000;

q0=1;q1=1;v=1;q=a+b*class  

cost=1000

for(l in 1:ite){

if(cost>0){  
  
#vec=q

#for(s in 1:length(q)){
  
#vec_sub=vec;vec_sub[s]=vec_sub[s]+h

#q[s]=q[s]+eta*(ABIC(vec_sub,q0,q1,v)-ABIC(vec,q0,q1,v))/h
  
#}
  
q0=q0-eta*(ABIC(q,q0+h,q1,v)-ABIC(q,q0,q1,v))/h

q1=q1-eta*(ABIC(q,q0,q1+h,v)-ABIC(q,q0,q1,v))/h

v=v-eta*(ABIC(q,q0,q1,v+h)-ABIC(q,q0,q1,v))/h

cost=ABIC(q,q0,q1,v)

}

print(ABIC(q,q0,q1,v))
  
}

ABIC(q,q0,q1,v)


```


```r

#p.119 密度関数の推定法

library(dplyr)

C=1000

d=seq(0,1,1/C)

n=100000

val=rpois(n,sample(10:20,1));val=val/max(val)

n_vec=c()

for(j in 2:(length(d)-1)){
  
n_vec=c(n_vec,length(val[val>=d[j-1] & val<d[j]]))  
  
}

n_vec=c(n_vec,length(val[val>=d[length(d)-1] & val<=d[length(d)]]))

h_vec=rep(1,(C-1));h_vec=c(h_vec,-sum(h_vec))



ite=100000



for(l in 1:ite){
  
dl=(n_vec-n_vec[length(n_vec)])/C-n*(exp(h_vec/C)-exp(h_vec[length(h_vec)]/C))/(C*sum(exp(h_vec/C))) 

dl=dl[1:(C-1)]
  
p_j=exp(h_vec/C)/sum(exp(h_vec/C));p_j[is.nan(p_j)>0]=1
  
mat=cbind(diag(sqrt(p_j[1:(C-1)])),rep(sqrt(p_j[C]),C-1))

p_vec=c(p_j-p_j[length(p_j)])[1:(length(p_j)-1)]
  
#G0=mat%*%t(mat)-t(t(p_vec))%*%t(p_vec)

#dl2=-n*G0/(C^2);

#H=svd(dl2)$v%*%(diag(1/svd(dl2)$d))%*%t(svd(dl2)$u)

#diag(dl2)=diag(dl2)+(10^(-2))

h_vec2=h_vec[1:(C-1)]+(10^(-1))*dl

h_vec[1:(C-1)]=h_vec2;h_vec[C]=-sum(h_vec2)

log_lik=sum(n_vec*log(exp(h_vec/C)/sum(exp(h_vec/C))))

print(log_lik)

}

cor(C*p_j,n_vec)

ABIC=function(q,q0,q1,v){
  
log_liklihood=sum(n_vec*log(exp(q/C)/sum(exp(q/C))))
  
val=c((C)*log(1/sqrt(2*pi*v^2)),-(q[1]-2*q0+q1)^2/(2*v^2),-(q[2]-2*q[1]+q0)^2/(2*v^2))

for(s in 3:(C)){
  
val=c(val,-(q[j]-2*q[j-1]+q[j-2])^2/(2*v^2))  
  
}

#D=array(0,dim=c(C-1,C-1))

#for(j in 1:(C-3)){
  
#D[j:(j+2),j]=c(1,-2,1)  
  
#}

#D[(C-2),(C-2)]=1

#D[(C-1),(C-3):(C-1)]=c(1,-2,1);D=rbind(D,c(rep(1,C-3),c(0,3)))

p_j=exp(q/C)/(1+exp(q/C))  

#mat=cbind(diag(sqrt(p_j[1:(C-1)])),rep(sqrt(p_j[C]),C-1))

#G_hat=mat%*%t(mat)*(n/(C^2));H=G_hat+t(D)%*%D/(v^2)

return(-2*(log_liklihood+(sum(val))))

}


ite=10000;h=0.01;eta=0.0001

q0=1;q1=1;v=1

cost=1000;q=h_vec

for(l in 1:ite){

if(cost>0){  
  
#vec=q

#for(s in 1:length(q)){
  
#vec_sub=vec;vec_sub[s]=vec_sub[s]+h

#q[s]=q[s]+eta*(ABIC(vec_sub,q0,q1,v)-ABIC(vec,q0,q1,v))/h
  
#}
  
cost=ABIC(q,q0,q1,v)  
  
q0=q0-eta*(ABIC(q,q0+h,q1,v)-cost)/h

q1=q1-eta*(ABIC(q,q0,q1+h,v)-cost)/h

v=v-eta*(ABIC(q,q0,q1,v+h)-cost)/h

#cost=ABIC(q,q0,q1,v)

}

print(cost)
  
}




```