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

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