reiser.2000 <- function(Y1,Y2){
m <- nrow(Y1)
n <- nrow(Y2)
p <- ncol(Y1)
y1_bar <- rowMeans(Y1)
y2_bar <- rowMeans(Y2)
mu1.hat <- mean(Y1)
mu2.hat <- mean(Y2)
SSB1 <- p * sum( (y1_bar-mu1.hat)^2)
SSB2 <- p * sum( (y2_bar-mu2.hat)^2)
tau1_hat <-sqrt(( 1 / ((p-1)*m)) * sum((Y1-do.call(cbind,replicate(p,y1_bar,simplify = FALSE)))^2))
tau2_hat <-sqrt(( 1 / ((p-1)*n)) * sum((Y2-do.call(cbind,replicate(p,y2_bar,simplify = FALSE)))^2))
# Reiser 2000, just above 3.2
sig1 <- sqrt(((1/(n-1)) * SSB1 - tau1_hat^2)/p)
sig2 <- sqrt (((1/(n-1)) * SSB2 - tau2_hat^2)/p)
delta.hat <- pnorm((mu2.hat-mu1.hat)/(sqrt(sig1^2+sig2^2)))
delta.hat
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.