n <- c(50,30,40) ng <- length(n) N <- sum(n)
I <- diag(N) hacerJ <- function (d) matrix(1/d, d, d)
J <- hacerJ(N)
D <- as.matrix(Matrix::bdiag(lapply(n, hacerJ)))
W <- I-D B <- D-J mu <- runif(1, -15, 15); sigma <- runif(1, 1, 10)
x <- rnorm(N, mu, sigma)
SCE <- t(x) %*% W %*% x
SCF <- t(x) %*% B %*% x
c(SCE, SCF)
g <- factor(rep(1:ng, n))
sal <- aov(x ~ g)
summary(sal)
chi <- replicate(1000, {
x <- rnorm(N, mu, sigma)
SCE <- t(x) %*% W %*% x
SCF <- t(x) %*% B %*% x
N*log(1+SCF/SCE)}) c(mean(chi), sd(chi))
ks.test(chi, pchisq, ng-1)