推薦システムーマトリクス分解の多彩なすがたー 共立出版 廣瀬英雄先生

#p.19 IRT(項目反応理論) 確率勾配法

n=10
N=100

mat=array(0,dim=c(N,n))
for(j in 1:n){
vec=sample(0:1,N,replace=T,prob = prob)
mat[,j]=vec
}

#モデル選択(どれかを１つ選択)
logit=0　#ロジスティック
probit=0　#プロビット
comple_tran=1　#補対数対数モデル


loglik=function(p,no){
alpha=p[1:n];p=p[-c(1:n)]
beta=p[1:n];p=p[-c(1:n)]
loglik1=mat;loglik2=1-mat
for(i in no){

if(logit==1){  
loglik1[i,]=loglik1[i,]*log(1/(1+exp(-1.7*alpha*(p[i]-beta))))
loglik2[i,]=loglik2[i,]*log(1-1/(1+exp(-1.7*alpha*(p[i]-beta))))
}
  
if(probit==1){  
loglik1[i,]=loglik1[i,]*log(pnorm(alpha*(p[i]-beta)))
loglik2[i,]=loglik2[i,]*log(1-pnorm(alpha*(p[i]-beta)))
}
  
if(comple_tran==1){  
loglik1[i,]=loglik1[i,]*log(exp(-exp(alpha*(p[i]-beta))))
loglik2[i,]=loglik2[i,]*log(1-exp(-exp(alpha*(p[i]-beta))))
}  

}    
return(sum(loglik1)+sum(loglik2))    
}


param=rep(0.001,2*n+N)

ite=10^5
eta=0.001
h=0.01
p=0.3

Numbers=sample(1:N,ite,replace=T)  

for(l in 1:ite){

dL=c()
for(j in 1:length(param)){
param_sub=param;param_sub[j]=param_sub[j]+h  
dL=c(dL,(loglik(param_sub,Numbers[l])-loglik(param,Numbers[l]))/h)
}

param=param+eta*dL
print(loglik(param,c(1:N)))    
}

a=param[1:n];param=param[-c(1:n)]
b=param[1:n];param=param[-c(1:n)]
sita=param

#精度の確認
mat1=mat;mat2=1-mat
for(i in 1:N){
  
if(logit==1){  
mat1[i,]=mat1[i,]*ifelse(1/(1+exp(-1.7*a*(sita[i]-b)))>0.5,1,0)
mat2[i,]=mat2[i,]*ifelse(1/(1+exp(-1.7*a*(sita[i]-b)))>0.5,0,1)
}
  
if(probit==1){  
mat1[i,]=mat1[i,]*ifelse(pnorm(a*(sita[i]-b))>0.5,1,0)
mat2[i,]=mat2[i,]*ifelse(pnorm(a*(sita[i]-b))>0.5,0,1)
}
  
if(comple_tran==1){  
mat1[i,]=mat1[i,]*ifelse(exp(-exp(a*(sita[i]-b)))>0.5,1,0)
mat2[i,]=mat2[i,]*ifelse(exp(-exp(a*(sita[i]-b)))>0.5,0,1)
}  
  
}
(sum(mat1)+sum(mat2))/prod(dim(mat))

#p.52 マトリックス分解
#Ridge

R=matrix(c(3,3,3,1,NA,3,3,2,3,3,4,5),nrow=4)

k=5
U=matrix(sample(1:100,k*nrow(R),replace = T)/100,nrow=nrow(R))
V=matrix(sample(1:100,k*ncol(R),replace = T)/100,nrow=ncol(R))

ite=500
ku=1
kv=1
eta=0.001

for(l in 1:ite){
mat=R-U%*%t(V)  
mat[is.na(mat)==T]=0  
U=U-eta*(-2*(mat)%*%V+2*ku*U)
V=V-eta*(-2*t(mat)%*%U+2*ku*V)
S=sum((mat)^2)+ku*sum(U^2)+kv*sum(V^2)
RMSE=sqrt(sum(mat^2)/length(c(R[is.na(R)!=T])))
print(S)
}

#p.52 マトリックス分解
#lasso

R=matrix(c(3,3,3,1,NA,3,3,2,3,3,4,5),nrow=4)

k=5
U=matrix(sample(1:100,k*nrow(R),replace = T)/100,nrow=nrow(R))
V=matrix(sample(1:100,k*ncol(R),replace = T)/100,nrow=ncol(R))

ite=10000
ku=0.1
kv=0.1
eta=0.001

for(l in 1:ite){
mat=R-U%*%t(V)  
mat[is.na(mat)==T]=0  

valU=U-eta*(-2*(mat)%*%V)
valV=V-eta*(-2*t(mat)%*%U)

U=ifelse(abs(valU)<=ku*eta,0,ifelse(valU>ku*eta,valU-ku*eta,valU+ku*eta))
V=ifelse(abs(valV)<=kv*eta,0,ifelse(valV>kv*eta,valV-kv*eta,valV+kv*eta))

S=sum((mat)^2)+ku*sum(abs(U))+kv*sum(abs(V))
RMSE=sqrt(sum(mat^2)/length(c(R[is.na(R)!=T])))
print(S)
}

#p.58 非負マトリックス分解

R=matrix(c(3,3,3,1,3,3,3,2,3),nrow=3)

k=4

U=matrix(rep(0.01,k*nrow(R)),nrow=nrow(R))
V=matrix(rep(0.01,k*nrow(R)),nrow=nrow(R))

ite=100
for(l in 1:ite){
U=(R%*%V)*U/(U%*%t(V)%*%V)  
V=t(R)%*%U*V/(V%*%t(U)%*%U)
print(sum((R-U%*%t(V))^2))
}

#p.89 ハイブリッド法(lassoと非負マトリクス分解)

R=matrix(c(3,3,3,1,2,3,3,2,3,3,4,5),nrow=4)

k=5
U=matrix(sample(1:100,k*nrow(R),replace = T)/100,nrow=nrow(R))
V=matrix(sample(1:100,k*ncol(R),replace = T)/100,nrow=ncol(R))
UU=matrix(sample(1:100,k*nrow(R),replace = T)/100,nrow=nrow(R))
VV=matrix(sample(1:100,k*ncol(R),replace = T)/100,nrow=ncol(R))
alasso=0.01
adecom=0.01

ite=10000
ku=0.1
kv=0.1
eta=0.001

for(l in 1:ite){
#lasso
mat=R-U%*%t(V)  
valU=U-eta*(-2*(mat)%*%V)
valV=V-eta*(-2*t(mat)%*%U)
U=ifelse(abs(valU)<=ku*eta,0,ifelse(valU>ku*eta,valU-ku*eta,valU+ku*eta))
V=ifelse(abs(valV)<=kv*eta,0,ifelse(valV>kv*eta,valV-kv*eta,valV+kv*eta))
Slasso=sum((R-U%*%t(V))^2)

#非負マトリクス分解法
UU=(R%*%VV)*UU/(UU%*%t(VV)%*%VV)  
VV=t(R)%*%UU*VV/(VV%*%t(UU)%*%UU)
Sdecom=sum((R-UU%*%t(VV))^2)

alasso=alasso-2*eta*exp(alasso)*sum(R-U%*%t(V))
adecom=adecom-2*eta*exp(adecom)*sum(R-UU%*%t(VV))

cost=sum((exp(alasso)*(R-U%*%t(V))+exp(adecom)*(R-UU%*%t(VV)))^2)
print(cost)
}