#CONJUGATE GRADIENT FOR THE BINARY REGRESSION PROBLEM

#Simulation
library(MASS,MCMCpack)
set.seed(23432)
n<-50
p<-2
Sig<-rWishart(1,p+1,diag(1,p))[,,1]
X<-mvrnorm(n,mu=c(0,0),Sigma=0.1*Sig)
be<-c(1,-3,4)
X<-cbind(rep(1,n),X)
eta<-X %*% be
muvec<-1/(1+exp(-eta))
Y<-rbinom(n,1,muvec)

BinaryData<-data.frame(Y,X1=X[,1],X2=X[,2],X3=X[,3])

fBD<-glm(Y~X1+X2+X3-1,family=binomial,data=BinaryData)
summary(fBD)

#############################################################

link.fn<-function(tv){
	
	return(log(tv/(1-tv)))

}

h.fn<-function(tv){
	
	return(1/(1+exp(-tv)))

}

dot.h.fn<-function(tv){
	
	return(exp(tv)/((1+exp(tv))^2))

}

ddot.h.fn<-function(tv){
	
	return(exp(tv)*(1-exp(tv))/(1+exp(tv))^3)

}

V.fn<-function(muv){

	return(muv*(1-muv))

}

ell<-function(bev,Yv,Xm){
	etavec<-Xm %*% bev
	muvec<-h.fn(etavec)
	lterms<-Yv*log(muvec)+(1-Yv)*log(1-muvec)
	return(apply(lterms,2,sum))
}


dot.ell<-function(bev,Yv,Xm){
	etavec<-Xm %*% bev
	muvec<-h.fn(etavec)
	Vvec<-V.fn(muvec)
	mudot<-apply(Xm,2,function(x){return(x*dot.h.fn(etavec))})
	lterms<-apply(mudot,2,function(x){return(x*(Yv-muvec)/Vvec)})
	return(apply(lterms,2,sum))
}

dot.ell(c(0,0,0),Y,X)
dot.ell(coef(fBD),Y,X)


################################################
#Non-linear conjugate gradient method

#To do minimization: use g = -l, dot g = - dot l.

#Initialize
old.theta<-rep(0,p+1)
old.theta[1]<-log(mean(Y)/(1-mean(Y)))

#Line search

alpha.vec<-seq(-10,10,by=0.01)
Lmat<-outer(alpha.vec,dot.ell(old.theta,Y,X),"*")
thvec<-t(apply(Lmat,1,function(x){return(x+old.theta)}))
line.val<-apply(thvec,1,function(x){return(-ell(x,Y,X))})

#old alpha is the STEP SIZE.
old.alpha<-alpha.vec[which.min(line.val)]

#old d is the DIRECTION of the increment
old.d<-dot.ell(old.theta,Y,X)

#Update
new.theta<-old.theta+old.alpha*dot.ell(old.theta,Y,X)
old.beta<-sum(dot.ell(new.theta,Y,X)^2)/sum(dot.ell(old.theta,Y,X)^2)

old.theta<-new.theta
epsilon<-Inf
abstol<-1e-8
iter<-0
while(epsilon > abstol){

	old.d<-dot.ell(old.theta,Y,X)+old.beta*old.d

	Lmat<-outer(alpha.vec,dot.ell(old.theta,Y,X),"*")
	thvec<-t(apply(Lmat,1,function(x){return(x+old.theta)}))
	line.val<-apply(thvec,1,function(x){return(-ell(x,Y,X))})

	old.alpha<-alpha.vec[which.min(line.val)]
	new.theta<-old.theta+old.alpha*dot.ell(old.theta,Y,X)
	old.beta<-sum(dot.ell(new.theta,Y,X)^2)/sum(dot.ell(old.theta,Y,X)^2)

	epsilon<-sqrt(sum((new.theta-old.theta)^2))

	old.theta<-new.theta
	lval<-ell(old.theta,Y,X)

	iter<-iter+1

	print(c(iter,epsilon,old.theta,lval))

}