#CONJUGATE GRADIENT FOR THE BINARY REGRESSION PROBLEM 2D

#Simulation
library(MASS,MCMCpack)
set.seed(23432)
n<-100
p<-2
X<-rnorm(n)
be<-c(-2,4)
X<-cbind(rep(1,n),X)
eta<-X %*% be
muvec<-1/(1+exp(-eta))
Y<-rbinom(n,1,muvec)

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

fBD2<-glm(Y~X1+X2-1,family=binomial,data=BinaryData2D)
summary(fBD2)

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


log.like.br<-function(bev,Xm,Yv){

	ev<-Xm %*% bev
	pv<-1/(1+exp(-ev))
	tvec<-sum(dbinom(Yv,1,pv,log=T))
	return(tvec)

}


be0<-be1<-seq(-10,10,by=0.1)
lmat<-matrix(0,nrow=length(be0),ncol=length(be1))
for(i in 1:length(be0)){
	for(j in 1:length(be1)){
	
		lmat[i,j]<-log.like.br(c(be0[i],be1[j]),X,Y)

	}
	
}

image(be0,be1,lmat,col=heat.colors(50),xlab=expression(beta[0]),ylab=expression(beta[1]))
contour(be0,be1,lmat,nlevels=50,add=T,col='gray')
abline(v=coef(fBD2)[1],h=coef(fBD2)[2],lty=3)


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

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),Y,X)
dot.ell(coef(fBD2),Y,X)


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

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

#Initialize
old.theta<-c(7.5,-5)
#old.theta<-c(-5,-10)

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

contour(be0,be1,lmat,nlevels=50,xlab=expression(beta[0]),ylab=expression(beta[1]))
abline(v=coef(fBD2)[1],h=coef(fBD2)[2],col="red",lty=2)


points(old.theta[1],old.theta[2],col="red",pch=3,lwd=2)


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


lines(c(old.theta[1],new.theta[1]),c(old.theta[2],new.theta[2]),col="red")

old.theta<-new.theta

points(old.theta[1],old.theta[2],col="red",pch=3,lwd=2)


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

	lines(c(old.theta[1],new.theta[1]),c(old.theta[2],new.theta[2]),col="red")

	old.theta<-new.theta

	points(old.theta[1],old.theta[2],col="red",pch=3,lwd=2)

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

	iter<-iter+1

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

}
summary(fBD2)