# THERE ARE 8 THINGS "TO DO" BELOW. PLEASE DO THEM, AND HAND IN 
# THE RESULTS. IF YOU DON'T FINISH IT IN CLASS, YOU MAY WORK ON
# IT AT HOME AND HAND IT IN LATER.


# clear the memory
rm(list=ls())

#DATA
x <- c(-0.86,-0.30,-0.05,0.73)
n <- c(5,5,5,5)
y <- c(0,1,3,5)


# Make ys into 0s and 1s and use logistic regression. 
xs <- rep(x,each=5)
ys <- c(rep(0,5),1,rep(0,4),rep(1,3),rep(0,2),rep(1,5))

# Fit logistic reression. This is one way to do the analysis, but
# we're using it to find out where to put our grid of alphas and betas.

logist.reg.fit <- glm(ys~xs,family=binomial(link="logit"))

# Look at the logistic regression fit.

summary(logist.reg.fit)

# Establish a grid of alphas and betas on which to evaluate
# the posterior. The grid end points are the mode +/- 6 standard errors.

alpha.grid <- 
	seq(from=summary(logist.reg.fit)$coef[1,1]-
				6*summary(logist.reg.fit)$coef[1,2],
		to=summary(logist.reg.fit)$coef[1,1]+
				6*summary(logist.reg.fit)$coef[1,2],
		length=200)

beta.grid <- 
	seq(from=summary(logist.reg.fit)$coef[2,1]-
				6*summary(logist.reg.fit)$coef[2,2],
		to=summary(logist.reg.fit)$coef[2,1]+
				6*summary(logist.reg.fit)$coef[2,2],
		length=200)

# Build a function to calculate the posterior.

posterior <- function(alpha,beta,y,n,x)
{
	# It's numerically more stable to compute the log
	# likelihood, and exponentiate
	
	ps <- exp(alpha+beta*x)/(1+exp(alpha+beta*x))
	log.liks <- y*log(ps/(1-ps))+n*log(1-ps) 
	log.lik <- sum(log.liks)
	log.lik <- log.lik+6
	# Adding 6 makes this the maximum of this likelihood close to 1. 
	# That makes sure we don't have numerical round-off errors. 
		
	lik <- exp(log.lik)	
	return(lik)	
}


# evaluate the posterior on the grids of alphas and betas
post.evals <- matrix(-999,200,200) # where to keep the posterior

for (i in 1:200)
{
	for (j in 1:200)
	{
		post.evals[i,j] <- 
			posterior(alpha.grid[i],beta.grid[j],y,n,x)
	}
}
# Note: each row of post.evals contains a different alpha 
# and each column has a different beta.



layout(matrix(c(1,3,2,4),2,2,byrow=T),widths=c(1,1),heights=c(1,1))


# draw the contour plot
contour(alpha.grid,beta.grid,post.evals,xlab="alpha",ylab="beta",
		levels=round(seq(from=.05,to=.95,length=5)*max(post.evals),1))


# apply() is a useful function. It sums over the rows or columns of a matrix.
# To see how it works:

mat <- matrix(rep(1,10),5,2)

mat
apply(mat,1,sum)
apply(mat,2,sum)


# TO DO 1: Please use the apply function to create p(alpha|y) and p(beta|y)
# from post.evals. (hint: p(alpha|y) = int p(alpha,beta|y) dbeta)
#

post.alpha <- 
post.beta <- 

# Plot the marginal posteriors of alpha and beta.

plot(alpha.grid,post.alpha,typ="l",xlab="alpha",ylab="p(alpha|y)")
plot(beta.grid,post.beta,typ="l",xlab="beta",ylab="p(beta|y)")


# Simulate from p(alpha,beta|y) and plot the simulations.

plot(alpha.grid,beta.grid,type="n",xlab="alpha",ylab="beta")
NN <- 200 # number of simulations
keep <- matrix(-999,NN,2) # to keep the simulations

for (i in 1:NN)
{
	# Generate random alpha from marginal posterior of alpha.
	# Note: post.alpha doesn't add up to 1. The rmultinom() function
	# automatically normalizes.
	alpha.star <- 
		alpha.grid[which.max(as.vector(rmultinom(1,1,post.alpha)))]
	# Please read the help files for rmultinom and which.max
	# TO DO 2: What does rmultinom(1,1,post.alpha) produce?
	# TO DO 3: What does which.max(as.vector(rmultinom(1,1,post.alpha))) produce?

	post.beta.given.y.alpha.star <- rep(-999,200)
	for (j in 1:200)
		post.beta.given.y.alpha.star[j] <- posterior(alpha.star,beta.grid[j],y,n,x)
	# TO DO 4: Why can we use the posterior for p(beta|alpha,y)? 
	
	# generate random beta given alpha.star and y
	beta.star <- 
		beta.grid[which.max(as.vector(rmultinom(1,1,post.beta.given.y.alpha.star)))]
	points(alpha.star,beta.star,pch=16,cex=.5)
	keep[i,] <- c(alpha.star,beta.star)

}
# TO DO 5: Why does the algorithm above sample from the posterior?

# TO DO 6: Please estimate the posterior mean of alpha and beta.
# TO DO 7: What are approximate 95% intervals for alpha and beta?
# TO DO 8: Modify the code to change the prior for alpha and beta to 
# independent normals with mean 0 and variance 0.001, and reproduce
# the figure. How does the posterior change? What is a 
# normal prior that produces a posterior that looks similar to
# the prior p(alpha,beta) = 1?