S <- matrix(c(521, 477, 484, 510, 523, 528,
              477, 576, 536, 575, 580, 584,
              484, 536, 601, 593, 598, 613,
              510, 575, 593, 755, 718, 722,
              523, 580, 598, 718, 797, 751,
              528, 584, 613, 722, 751, 802), 6, 6)

N2 <- 152/2

Kprime <- matrix( c(1,0,0,0,0,0,
					1,1,0,0,0,0,
					1,1,1,0,0,0,
					1,1,1,1,0,0,
					1,1,1,1,1,0,
					1,1,1,1,1,1), 6, 6)

K <- t(Kprime)
Imatrix <- diag(1, 6, 6)

ll.BB <- function(theta, S) {
	Delta <- diag(theta[1:6])
	gamma <- theta[7]
	KDK <- K %*% Delta %*% Kprime
	Sigma <- KDK + gamma*Imatrix
	SigmaInverse <- solve(Sigma)
	SigmaS <- SigmaInverse %*% S
	l <- N2 * (log(det(Sigma)) + (sum(diag(SigmaS))))
	
	W <- SigmaS %*% SigmaInverse
	Psi1Delta <- N2 * diag(Kprime %*% (SigmaInverse - W) %*% K)
	Psi1gamma <- N2 * (sum(diag(SigmaInverse - W)))
	Psi1 <- c(Psi1Delta, Psi1gamma)

	Psi2Delta <- N2 * ((Kprime %*% SigmaInverse %*% K) 
					* (Kprime %*% (2*W - SigmaInverse) %*% K))
	Psi2Deltagamma <- N2 * (diag(Kprime %*% SigmaInverse 
									%*% (2*W - SigmaInverse) %*% K))
	Psi2gamma <- N2 * sum(diag(SigmaInverse %*% (2*W - SigmaInverse)))

	Psi2top <- cbind(Psi2Delta,Psi2Deltagamma)
	Psi2bottom <- cbind(t(Psi2Deltagamma),Psi2gamma)
	Psi2 <- rbind(Psi2top,Psi2bottom)
	
	attr(l, "gradient") <- Psi1
	attr(l, "hessian") <- Psi2
	return(l)
 }

theta <- c(500, 60, 20, 80, 20, 25, 50)

system.time( result <- nlm(f=ll.BB,p=theta,hessian=TRUE,S=S) )

SE <- sqrt(diag(solve(result$hessian)))