function(y, x)
{
	n <- length(y)
	k <- ncol(x)
	xs <- cbind(jay(n, 1), x)
	b <- inv(t(xs) %*% xs) %*% t(xs) %*% y
	yhat <- xs %*% b
	syx <- sqrt(sum((y - yhat)^2)/(n - k - 1))
	ds <- t(xs) %*% xs
	cookd <- c(1:n)
	a <- c(2:n)
	xs1 <- xs[a,  ]
	y1 <- y[a,  ]
	b1 <- inv(t(xs1) %*% xs1) %*% t(xs1) %*% y1
	cookd[1] <- t(b1 - b) %*% ds %*% (b1 - b)/(syx * (k +
		1))
	for(i in 2:(n - 1)) {
		a <- c(1:i - 1, i + 1, n)
		xs1 <- xs[a,  ]
		y1 <- y[a,  ]
		b1 <- inv(t(xs1) %*% xs1) %*% t(xs1) %*% y1
		cookd[i] <- t(b1 - b) %*% ds %*% (b1 - b)/(
			syx * (k + 1))
	}
	a <- c(1:n - 1)
	xs1 <- xs[a,  ]
	y1 <- y[a,  ]
	b1 <- inv(t(xs1) %*% xs1) %*% t(xs1) %*% y1
	cookd[n] <- t(b1 - b) %*% ds %*% (b1 - b)/(syx * (k +
		1))
	cookd
}
