混合ロジスティックモデル

#mixed-logistic model 

#multi dimension

#見直し中

library(dplyr)

data=data.frame(num=1:20,y=c(0.20, 0.10, 0.49, 0.26, 0.92, 0.95, 0.18, 0.19, 0.44, 0.79, 0.61, 0.41, 0.49, 0.34, 0.62, 0.99, 0.38, 0.22,0.71,0.40),x1=c(84,87,86,85,82,87,92,94,88,84.9,78,90,88,87,82,84,86,83,85.2,82),x2=c(61,55.5,57,57,50,50,66.5,65,60.5,49.5,49.5,61,59.5,58.4,53.5,54,60,58.8,54,56),x3=c(24.5,24,23,22.5,22,24,25,26,25,24,23,21,26,25.5,24,23,24,24,24,22),sign=0)

for(j in 1:nrow(data)){

data$sign[j]=ifelse(data$y[j]>0.49,1,0)

}

X=as.matrix(cbind(rep(1,nrow(data)),data[,colnames(data) %in% c("x1","x2")]))

Y=as.numeric(data$sign)




#Broyden Quasi-Newton Algorithm p.213

f=function(x){

pis=x[1:K];x=x[-c(1:K)];alpha=matrix(c(x),ncol=K)

pt=array(0,dim=c(K,nrow(X)));mat=array(0,dim=c(K,nrow(X)))

for(i in 1:K){
for(j in 1:nrow(X)){

pt[i,j]=pis[i]/(1+exp(-sum(alpha[,i]*X[j,])))

mat[i,j]=exp(-sum(alpha[,i]*X[j,]))/(1+exp(-sum(alpha[,i]*X[j,])))

}}

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

vectors=array(0,dim=c(K,ncol(X)))

for(i in 1:K){

vectors[i,]=apply(X*c((pt*mat)[i,]),2,sum)

}

return(c(t(vectors)))

}

K=2

param=c(sample(seq(0.01,1,0.01),K),rep(0.01,K*ncol(X)))

XX=param

ite=10^7;

eta=10^(-5)

H=diag(f(XX))

for(l in 1:ite){

if(sum(abs(f(XX)))>10^(-9)){  

X_pre=XX  

XX[-c(1:K)]=XX[-c(1:K)]-eta*H%*%f(XX)  

s=XX[-c(1:K)]-X_pre[-c(1:K)]

y=f(XX)-f(X_pre)

H=H+((s-H%*%y)/as.numeric(t(s)%*%H%*%y))%*%t(s)%*%H

#print(sum(f(XX)))

ps=XX[1:K];XX=XX[-c(1:K)];a=matrix(XX,ncol=K)

Pt=array(0,dim=c(K,nrow(X)))

log_f=array(0,dim=c(K,nrow(X)))

for(i in 1:K){
for(j in 1:nrow(X)){

Pt[i,j]=ps[i]/(1+exp(-sum(X[j,]*a[,i])))

log_f[i,j]=-log(1+exp(-sum(X[j,]*a[,i])))

}}

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

ps=apply(Pt,1,mean)

Q=sum(Pt*log_f)+sum(Pt*c(log(ifelse(ps>0,ps,1))))

XX=c(ps,XX)

print(Q)

}

}