welch-t-test.R

## Welch t-test

### Simple Normal-Normal example
mu1 <- 1
mu2 <- 2
var1 <- 1
var2 <- 0.5
n <- 10^3
y1 <- rnorm(n, mu1, sqrt(var1))
y2 <- rnorm(n, mu2, sqrt(var2))

## base implementation
(tt <- t.test(y1, y2, var.equal=FALSE))

## by hand
mu1 <- mean(y1)
mu2 <- mean(y2)
var1 <- sd(y1)^2
var2 <- sd(y2)^2
(tval <- (mu1-mu2) / sqrt(var1/n + var2/n))
(df <- ((var1/n) + (var2/n))^2 / (((var1/n)^2 / (n-1)) + ((var2/n)^2 / (n-1))))
(pval <- 2*pt(-abs(tval), df))

## a function
compare_distr <- function(mu1, var1, mu2, var2) {
    tval <- (mu1-mu2) / sqrt(var1/n + var2/n)
    df <- ((var1/n) + (var2/n))^2 / (((var1/n)^2 / (n-1)) + ((var2/n)^2 / (n-1)))
    pval <- 2*pt(-abs(tval), df)
    c(t=tval, df=df, p=pval)
}
compare_distr(mu1, var1, mu2, var2)