setwd("C:\\Work\\Courses\\Courses\\Current\\666\\Rpackages")

expit<-function(x){return(1/(1+exp(-x)))}
logit<-function(x){return(log(x/(1-x)))}

gbd<-function(n=1000,p=1,intercept=TRUE,betas=rep(0,p+intercept),mu=rep(0,p),Sigma=diag(1,p),seed=1){
	require(MASS)
	set.seed(seed)
	X<-mvrnorm(n,mu,Sigma)
	Xm<-X
	if(intercept){Xm<-cbind(rep(1,n),X)}
	lp<-Xm %*% betas
	pv<-expit(lp)
	Y<-rbinom(n,1,pv)
	return(list(Y=Y,X=X))
}

sigmat<-rWishart(1,3,diag(rep(1,2)))[,,1]

d1<-gbd(n=100,p=2,betas=c(0.2,-0.5,0.5),mu=rep(0,2),Sigma=sigmat,seed=101)

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

lr.loglike<-function(betas,Y,X,intercept=TRUE){
	if(!intercept){
		lp<-X %*% betas
	}else{
		lp<-cbind(rep(1,length(Y)),X) %*% betas
	}
	pv<-expit(lp)
	return(sum(dbinom(Y,1,pv,log=T)))
}

est.bes<-function(Y,X,intercept=TRUE,start.vals=rep(0,ncol(X)+as.numeric(intercept))){
	
	x<-start.vals
	opt<-optim(x,fn=lr.loglike,gr=NULL,Y=Y,X=X,intercept=intercept,method="BFGS",control=list(fnscale=-1),hessian=T)
	beta.hat<-opt$par
	beta.cov<-solve(opt$hessian)/length(Y)
	if(intercept){
		Xm<-cbind(rep(1,length(Y)),X)
	}else{
		Xm<-X
	}
	fitted.lp<-Xm %*% beta.hat
	fitted.Y<-expit(fitted.lp)
	conv<-TRUE
	if(opt$convergence!=0){
		conv<-FALSE
		print('Optim did not converge')
	}
	
	return(list(ests=beta.hat,varcov=beta.cov,fitted.lp=fitted.lp,fitted.Y=fitted.Y,conv=conv))
}

result<-est.bes(d1$Y,d1$X)

sigmat<-rWishart(1,3,diag(rep(1,1)))[,,1]
d1<-gbd(n=100,p=1,betas=c(-0.5,0.5),mu=rep(0,1),Sigma=sigmat)
result<-est.bes(d1$Y,d1$X)


b0<-seq(-2,2,by=0.01)
b1<-seq(-2,2,by=0.01)

lrs<-outer(b0,b1,"+")*0
for(i in 1:length(d1$Y)){
	lpv<-outer(b0,b1*d1$X[i],"+")
	lrs<-lrs+lpv*d1$Y[i]-log(1+exp(lpv))
}
image(b0,b1,lrs)
contour(b0,b1,lrs,add=T)

plot.lrsurf<-function(Y,X,b0=seq(-2,2,by=0.01),b1=seq(-2,2,by=0.01),ncl=20){
	lrs<-outer(b0,b1,"+")*0
	for(i in 1:length(Y)){
		lpv<-outer(b0,b1*X[i],"+")
		lrs<-lrs+lpv*Y[i]-log(1+exp(lpv))
	}
	image(b0,b1,lrs,xlab=expression(beta[0]),ylab=expression(beta[1]))
	contour(b0,b1,lrs,add=T,nlevels=ncl)
}

plot.lrsurf(d1$Y,d1$X,b0=seq(-1,0,by=0.005),b1=seq(0.5,1.5,by=0.005),ncl=50)
abline(v=result$ests[1],h=result$ests[2],lty=2)

package.skeleton(name="LogReg",list=c("logit","expit","gbd","est.bes","lr.loglike","est.bes","plot.lrsurf","test.data"))