accidents = read.table("florida_acc.txt")

SSE = function(n = 100000, hyper = c(0.01, 0.01, 0, 1000, 0.01, 0.01), accidents = accidents){

n_m = length(unique(accidents[,1])) #Just the number of months.
n_s = dim(accidents)[2]-1   #Just the number of sites. 

a = hyper[1]
b = hyper[2]
d = hyper[3]
e = hyper[4]
f = hyper[5]
g = hyper[6]

start.gamma = rep(0, n_m)
start.epsilon = rep(d, n_s)

err.gamma = rep(0.05, n_m)
err.epsilon = rep(0.05, n_s)

X = matrix(ncol = n_m, nrow = n)   #This matrix will take the MCMC output for the seasonal effects (gamma's) for log(lambda) - one row per iteration, one column per season

Y = matrix(ncol = n_s, nrow = n)  #This matrix will take the MCMC output for the site effects (epsilon's) for log(lambda) - one row per iteration, one column per site

can.X = matrix(ncol = n_m, nrow = n)  #This matrix will take the candidates for the seasonal effects
can.X[1, ] = start.gamma

can.Y = matrix(ncol = n_s, nrow = n)  #This matrix will take the candidates for the site effects
can.Y[1, ] = start.epsilon

#Now we finalise the first row of the MCMC matrices at their initial states:
X[1, ] = can.X[1, ]
Y[1, ] = can.Y[1, ]

#Matrices for acceptance probabilities:
aprob.gamma = matrix(ncol = n_m, nrow = n)
aprob.epsilon = matrix(ncol = n_s, nrow = n)

#Log-likelihood function:
loglik = function(GAMMA, EPSILON, DATA){
	llh = -length(DATA)*exp(GAMMA+EPSILON)+length(DATA)*mean(DATA)*(GAMMA+EPSILON) - sum(log(factorial(DATA)))
	return(llh)}

#Initialise tau, c and zeta using the starting values for the random effects:
tau = vector("numeric", n)
c = vector("numeric", n)
zeta = vector("numeric", n)

tau[1] = rgamma(1, a+(n_m/2), b+0.5*sum(X[1, ]^2))
c[1] = rnorm(1, (d*e + (f/g)*sum(Y[1, ]))/(e+n_s*(f/g)), 1/sqrt(e+n_s*(f/g)))
zeta[1] = rgamma(1, f+(n_s/2), g+0.5*sum((Y[1, ]-c[1])^2))

#Now start the MCMC/Gibbs iterations:

for(i in 2:n){
	 
	#Seasonal effects::
	lik.X = vector("numeric", n_m)
	lik.X[] = 0
	for(m in 1:n_m){
		can.X[i, m] = X[(i-1), m] + rnorm(1, 0, err.gamma[m])

	for(k in 1:n_s){
		#Careful with you data - I'm using (k+1) here as the first column of my accident data was a month indicator (so no actual data)
		data.contribution = accidents[, (k+1)][accidents[, 1]==m]
		data.contribution = data.contribution[!is.na(data.contribution)]
		if(length(data.contribution)==0){lik.cont.X = 0}else{
		lik.cont.X = loglik(can.X[i, m], Y[(i-1), k], data.contribution) - loglik(X[(i-1), m], Y[(i-1), k], data.contribution)}
		if(is.na(lik.cont.X)){print(k)}
		lik.X[m] = lik.cont.X + lik.X[m]}

	aprob.gamma[i, m] = min(1, (exp(lik.X[m])*(dnorm(can.X[i, m], 0, 1/sqrt(tau[i-1])))/(dnorm(X[(i-1), m], 0, 1/sqrt(tau[i-1])))))

	u = runif(1)
	if(u < aprob.gamma[i, m]){
		X[i, m] = can.X[i, m]}else{X[i, m] = X[(i-1), m]}  }  

	#Site effects::
	lik.Y = vector("numeric", n_s)
	lik.Y[] = 0
	for(j in 1:n_s){
		can.Y[i, j] = Y[(i-1), j] + rnorm(1, 0, err.epsilon[j])

	for(k in 1:n_m){
		data.contribution = accidents[, (j+1)][accidents[, 1]==k]
		data.contribution = data.contribution[!is.na(data.contribution)]
		if(length(data.contribution)==0){lik.cont.Y = 0}else{
		lik.cont.Y = loglik(X[i, k], can.Y[i, j], data.contribution) - loglik(X[i, k], Y[(i-1), j], data.contribution)}
		lik.Y[j] = lik.cont.Y + lik.Y[j]}

	aprob.epsilon[i, j] = min(1, (exp(lik.Y[j])*(dnorm(can.Y[i, j], c[i-1], 1/sqrt(zeta[i-1])))/(dnorm(Y[(i-1), j], c[i-1], 1/sqrt(zeta[i-1])))))
	
	u = runif(1)
	if(u < aprob.epsilon[i, j]){
		Y[i, j] = can.Y[i, j]}else{Y[i, j] = Y[(i-1), j]}   }

	tau[i] = rgamma(1, a+(n_m/2), b+0.5*sum(X[i, ]^2))
	c[i] = rnorm(1, (d*e + zeta[i-1]*sum(Y[i, ]))/(e+n_s*zeta[i-1]), 1/sqrt(e+n_s*zeta[i-1]))
	zeta[i] = rgamma(1, f+(n_s/2), g+0.5*sum((Y[i, ]-c[i])^2))  #### All updates now complete!  

	print(i)
		}

	return(list(GAMMA = X, EPSILON = Y, TAU = tau, C = c, ZETA = zeta, APROB.GAMMA = aprob.gamma, APROB.EPSILON = aprob.epsilon))  
	}