bigAhat <- - 1/n * t(X.int) %*% diag(Y/muhat^2) %*% X.int
bigBhat <- 1/n * t(X.int) %*% diag(1 - 2*Y/muhat + Y^2/muhat^2) %*% X.int
sandvar <- diag(solve(bigAhat) %*% bigBhat %*% solve(t(bigAhat)))
sand.lower.bound <- Bhat - qnorm(0.975)*sqrt(sandvar)
sand.upper.bound <- Bhat + qnorm(0.975)*sqrt(sandvar)
if ((sand.lower.bound < B) && (B < sand.upper.bound)) count.sand = count.sand + 1