項目反応理論

```{r}
#項目反応理論

#PROX 推定法　p.57

#ラッシュモデル

rm(list=ls())

N=100;n=10

u=array(0,dim=c(N,n))

for(j in 1:nrow(u)){
  
val=rpois(2,sample(6:10,1))  
  
u[j,]=sample(0:1,n,replace=T,prob=val/sum(val))  
  
}





log_lik=function(b,sita){
  
return(sum((u*c(sita)-t(t(u)*c(b)))-log(1+exp(matrix(rep(sita,n),nrow=N)-c(b)))))
  
}

b=rep(1,n);sita=rpois(N,sample(10:20,1))

eta=0.00001

ite=100000

for(j in 1:ite){

vec=b;h=0.01  
    
for(i in 1:length(b)){
  
b_sub=vec;b_sub[i]=b_sub[i]+h

dL=(log_lik(b_sub,sita)-log_lik(b,sita))/h

b[i]=b[i]+eta*dL
  
}

print(log_lik(b,sita))    

}


```

```{r}

#PROX 推定法　p.67

#ラッシュモデル

rm(list=ls())

N=100;n=10

u=array(0,dim=c(N,n))

for(j in 1:nrow(u)){
  
val=rpois(2,sample(6:10,1))  
  
u[j,]=sample(0:1,n,replace=T,prob=val/sum(val))  
  
}





log_lik=function(a,C,g){
  
P_box=array(0,dim=c(N,n))

Q_box=array(0,dim=c(N,n))

for(i in 1:N){
for(j in 1:n){
  
P_box[i,j]=g[j]+(1-g[j])*exp(a[j]*sita[i]+C[j])/(1+exp(a[j]*sita[i]+C[j]))

Q_box[i,j]=1-P_box[i,j]
  
}  
}

return(sum(u*log(P_box)+(1-u)*log(Q_box)))

}

a=rep(0.1,n);sita=rpois(N,sample(10:20,1))

C=rep(0.2,n);g=rep(0.5,n)

eta=0.00001

ite=100000

mesh=0.01

log_lik_val=-10^6;log_lik_val_pre=-10^6-1

for(j in 1:ite){
  
if(log_lik_val>log_lik_val_pre){  

log_lik_val_pre=log_lik_val  
  
vec=a;

for(i in 1:length(a)){
  
a_sub=vec;a_sub[i]=a_sub[i]+mesh

dL=(log_lik(a_sub,C,g)-log_lik(a,C,g))/mesh

a[i]=a[i]+eta*dL
  
}
  
vec=C;

for(i in 1:length(C)){
  
C_sub=vec;C_sub[i]=C_sub[i]+mesh

dL=(log_lik(a,C_sub,g)-log_lik(a,C,g))/mesh

C[i]=C[i]+eta*dL
  
}

P_box=array(0,dim=c(N,n))

Q_box=array(0,dim=c(N,n))

for(i in 1:N){
for(j in 1:n){
  
P_box[i,j]=g[j]+(1-g[j])*exp(a[j]*sita[i]+C[j])/(1+exp(a[j]*sita[i]+C[j]))

Q_box[i,j]=1-P_box[i,j]
  
}  
}

Q_star=Q_box/c(1-g);P_star=(P_box-c(g))/c(1-g)

dg=apply(u*(1-P_star)/P_box+(1-u)*(-Q_star)/(Q_box),1,sum)

g=g+eta*dg

log_lik_val=log_lik(a,C,g)

print(log_lik(a,C,g))

}
  
}  


```


```{r}


#同時最尤推定法 p.67

#Flecher and Powell 1963 多次元尺度法 p.238

library(dplyr)

n=100;m=15

u=array(0,dim=c(n,m))

for(j in 1:n){

val=rpois(2,sample(10:20,1))  

u[j,]=sample(0:1,m,replace=T,prob=val/sum(val))  

}

sita=rpois(n,sample(10:20,1))

sita=sita/max(sita)

a_vec=rep(1,m);c_vec=rep(1,m);g_vec=rep(0.99,m)

log_lik=function(a,c,g){

p_box=array(0,dim=c(n,m))

q_box=array(0,dim=c(n,m))

#sita_box=matrix(rep(sita,m),ncol=m)

for(i in 1:n){
for(j in 1:m){

p_box[i,j]=1-(abs(g[j])/sum(abs(g)))/(1+exp(a[j]*sita[i]+c[j]))

q_box[i,j]=(abs(g[j])/sum(abs(g)))/(1+exp(a[j]*sita[i]+c[j]))

}  
}

q_box[floor(p_box)==1]=0.001

p_box[floor(p_box)==1]=0.999;

return(sum(u*log(p_box)+(1-u)*log(q_box)))

}

log_lik(a_vec,c_vec,g_vec)





ite=100;h=10^(-3)

log_lik_vec=c(log_lik(a_vec,c_vec,g_vec))

log_lik_val_pre=-10^10;log_lik_val=-10^10+1

G=diag(1,length(a_vec)+length(c_vec)+length(g_vec))

alpha=0.1

for(l in 1:ite){
  
vec=a_vec;da=c()

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

a_sub=vec;a_sub[i]=a_sub[i]+h  

da=c(da,(log_lik(a_sub,c_vec,g_vec)-log_lik(a_vec,c_vec,g_vec))/h)

}  

vec=c_vec;dc=c() 

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

c_sub=vec;c_sub[i]=c_sub[i]+h  

dc=c(dc,(log_lik(a_vec,c_sub,g_vec)-log_lik(a_vec,c_vec,g_vec))/h)

}

vec=g_vec;dg=c() 

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

g_sub=vec;g_sub[i]=g_sub[i]+h  

dg=c(dg,(log_lik(a_vec,c_vec,g_sub)-log_lik(a_vec,c_vec,g_vec))/h)

}

da_pre=da;dc_pre=dc;dg_pre=dg

params=c(a_vec,c_vec,g_vec)-alpha*G%*%c(a_vec,c_vec,g_vec)

a_vec=params[1:length(a_vec)]

c_vec=params[(length(a_vec)+1):(length(a_vec)+length(c_vec))]

g_vec=params[(length(a_vec)+length(c_vec)+1):(length(a_vec)+length(c_vec)+length(g_vec))]

vec=a_vec;da=c()

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

a_sub=vec;a_sub[i]=a_sub[i]+h  

da=c(da,(log_lik(a_sub,c_vec,g_vec)-log_lik(a_vec,c_vec,g_vec))/h)

}  

vec=c_vec;dc=c() 

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

c_sub=vec;c_sub[i]=c_sub[i]+h  

dc=c(dc,(log_lik(a_vec,c_sub,g_vec)-log_lik(a_vec,c_vec,g_vec))/h)

}

vec=g_vec;dg=c() 

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

g_sub=vec;g_sub[i]=g_sub[i]+h  

dg=c(dg,(log_lik(a_vec,c_vec,g_sub)-log_lik(a_vec,c_vec,g_vec))/h)

}

y=c(da,dc,dg)-c(da_pre,dc_pre,dg_pre)

B=-G%*%(t(t(y))%*%t(y))%*%G/as.numeric(t(y)%*%G%*%t(t(y)))

d=-G%*%c(da_pre,dc_pre,dg_pre)

H=alpha*(d%*%t(d))/as.numeric(t(c(da_pre,dc_pre,dg_pre))%*%G%*%t(t(c(da_pre,dc_pre,dg_pre))))

G=G+H+B

print(log_lik(a_vec,c_vec,g_vec))

}


```

```{r}

#p.86 多値IRTモデル

data=data.frame(販売率=c(3.0,1.8,1.5,3.3,2.2,2.0,3.5,2.0,1.7,2.3),大都市=c(1,1,0,1,1,1,1,0,1,1),地方都市=c(0,1,1,0,1,1,0,1,0,0),温暖=c(0,1,0,1,1,0,1,1,0,1),普通=c(1,0,1,1,0,1,0,1,1,1),寒冷=c(1,1,0,0,1,1,0,0,1,1))
u=data[,-1]
colnames(u)=c("X1","X1","X2","X2","X2")
column_names=unique(colnames(u))

count=c()
for(j in 1:length(column_names)){
count=c(count,sum(colnames(u)==column_names[j]))  
}

loglik=function(x,n){
sita=x[1:nrow(u)];x=x[-c(1:nrow(u))]
a=x[1:length(column_names)];x=x[-c(1:length(column_names))]
b=x

Z=array(1,dim=dim(u))
Z=Z*sita
for(j in 1:length(column_names)){
for(i in 1:ncol(u)){
if(colnames(u)[i]==column_names[j]){
Z[,i]=exp(Z[,i]*a[j]-a[j]*b[i])  
}}}

P=Z
for(j in 1:length(column_names)){
matsub=Z[,colnames(u) %in% column_names[j]]
P[,colnames(u) %in% column_names[j]]=matsub/apply(matsub,1,sum)
}

logliksub=u*log(P)+(1-u)*log(1-P)
return(sum(logliksub[n,]))
}


param=rep(1,nrow(u)+length(column_names)+ncol(u))
param=param/length(param)
ite=1000
eta=10^(-2)
h=0.001
Num=sample(1:nrow(u),ite,replace = T)

for(l in 1:ite){
  dparam=c()
  for(i in 1:length(param)){
    param_sub=param
    param_sub[i]=param_sub[i]+h
    dparam=c(dparam,(loglik(param_sub,Num[l])-loglik(param,Num[l]))/h ) 
  }
  param=param+eta*dparam
  print(loglik(param,c(1:nrow(u))))
}



```