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

write.csv(BinaryData,"BinaryData.csv",row.names=F)

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

#Start for 2015

Bdata<-read.csv("BinaryData.csv")
Y<-Bdata$Y
X1<-Bdata$X1
X2<-Bdata$X2
X3<-Bdata$X3
n<-length(Y)

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


X<-cbind(X1,X2,X3)

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

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)
	return(sum(Yv*log(muvec)+(1-Yv)*log(1-muvec)))
}


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

#Test
dot.ell(c(0,0,0),Y,X)

ddot.ell<-function(bev,Yv,Xm){

	etavec<-Xm %*% bev
	muvec<-h.fn(etavec)
	Vvec<-V.fn(muvec)
	Term1<-t(Xm) %*% apply(Xm,2,function(x){return(x*ddot.h.fn(etavec)*(Yv-muvec)/Vvec)})
	
	Term21<-((Yv-muvec)*(1-2*muvec)+Vvec)/Vvec^2

	mudot<-apply(Xm,2,function(x){return(x*dot.h.fn(etavec))})
	Term21<-apply(mudot,2,function(x){return(x*Term21)})
	
	Term2<-matrix(0,length(bev),length(bev))
   	#Could do this loop better !
	for(i in 1:n){
		Term2<-Term2 + Term21[i,] %*% t(mudot[i,]) 
	}
	return(Term1-Term2)

}

Fisher.n<-function(bev,Yv,Xm){

	etavec<-Xm %*% bev
	muvec<-h.fn(etavec)
	Vvec<-V.fn(muvec)
	Fisher<-matrix(0,length(bev),length(bev))
	mudot<-apply(Xm,2,function(x){return(x*dot.h.fn(etavec))})
	for(i in 1:n){
		Fisher<-Fisher + mudot[i,] %*% t(mudot[i,]) /Vvec[i]
	
	}
	return(Fisher)

}

ddot.ell(c(0,0,0),Y,X)
lddot<-ddot.ell(coef(fBD),Y,X)
Fisher.n(coef(fBD),Y,X)

#Check Asymptotic Variance
sqrt(diag(solve(-lddot/n))/n)

summary(fBD)$coef[,2]

################################################
#Newton's method

p<-2
old.beta<-rep(0,p+1)
epsilon<-Inf
abstol<-1e-8
iter<-0
while(abs(epsilon) > abstol){
	new.beta<-old.beta - solve(ddot.ell(old.beta,Y,X)) %*% dot.ell(old.beta,Y,X)
	epsilon<-sqrt(sum((new.beta-old.beta)^2))
	old.beta<-new.beta
	dot.ellval<-dot.ell(new.beta,Y,X)
	iter<-iter+1
	print(c(iter,new.beta,epsilon,dot.ellval))
	Hessian<-ddot.ell(new.beta,Y,X)
}
#Check
print(old.beta)
coef(fBD)

#Replace ddot.ell(old.beta,Y,X) by Fisher.n(old.beta,Y,X)
#for Fisher scoring.


################################################
#Quasi-Newton method

p<-2
epsilon<-Inf
abstol<-1e-8
iter<-0

#Initialize
old.beta<-rep(0,3)
old.dot<-dot.ell(old.beta,Y,X)
old.M<--ddot.ell(old.beta,Y,X)

#Recursion until tolerance is met.
while(abs(epsilon) > abstol){

	#Newton-like update, with M replacing the second derivative matrix

	new.beta<-old.beta+solve(old.M) %*% old.dot
	new.dot<-dot.ell(new.beta,Y,X)

	old.g<-new.dot-old.dot
	old.s<-new.beta-old.beta

	old.v<-old.g+old.M%*%old.s
	old.c<-1/(t(old.v)%*%old.s)[1,1]

	#Rank-1 update for M
	old.M<-old.M-old.c*old.v%*%t(old.v)

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

	old.beta<-new.beta
	old.dot<-new.dot

	dot.ellval<-new.dot

	iter<-iter+1

	#Check convergence: print out
	#Parameter values, SSE, function value, gradient

	print(c(iter,new.beta,epsilon,ell(old.beta,Y,X),new.dot))
	Hessian<-ddot.ell(new.beta,Y,X)

}
#Check
print(old.beta)
coef(fBD)

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