項目反応理論

#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))    

}

#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))

}

}  

#同時最尤推定法 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))

}