n <- 4                                  #tamaño muestra 1
m <- 9                                  #tamaño muestra 2
mux    <- 0 ; muy    <- 0               #esperanzas
sigmax <- 1 ; sigmay <- 4               #desvíos
varianzas0 <- c (sigmax^2, sigmay^2)    #vector con varianzas
gl <- function (varianzas)              #cálculo de los grados de libertad
{
    rel   <- varianzas[2]/varianzas[1]  #relación entre las varianzas
    num_c <- (1/n + rel/m) ^ 2
    den_c <- 1/(n^3-n^2) + rel^2/(m^3-m^2)
    num_c/den_c
}
c0 <- gl(varianzas0)  # c teórico
## Simulación de la distribución del estadístico
nr <- 100000
M  <- replicate (nr,
{
    x <- rnorm (n, mux, sigmax) ; mx <- mean(x) ; vx <- var(x)
    y <- rnorm (m, muy, sigmay) ; my <- mean(y) ; vy <- var(y)
    t <- (mx - my) / sqrt (vx/n + vy/m)
})
M <- sort (M)
plot (M, pt(M,c0), type='l')            #comparación gráfica de distribuciones
lines (M, seq(0,1,length=nr), col='red')
## comparación de cuantiles
cuantiles <- c(.001, .01, .05, .1, .90, .95, .99, .999)
qt (cuantiles, c0)
quantile (M, cuantiles)
## Ejemplo para comparar "var.equal = TRUE" con "var.equal = FALSE"
## en el t.test de R
set.seed(1)
muy <- mux + 2*sigmax
x <- rnorm (n, mux, sigmax) ; y <- rnorm (m, muy, sigmay)
medias    <- c (mean(x), mean(y))
varianzas <- c ( var(x),  var(y))
esta      <- (medias[1] - medias[2]) / sqrt (1/m + 1/n) /
    sqrt ((varianzas[1]*(n-1) + varianzas[2]*(m-1)) / (n+m-2))
p_valor   <- 2 * (1 - pt (abs(esta), n+m-2)) #p-valor del contraste
c (estadístico=esta, gl=n+m-2, pv=p_valor)
test      <- t.test (x, y, var.equal = TRUE) #varianzas iguales
c (test$statistic, test$parameter, pv=test$p.value)
test      <- t.test (x, y)              #varianzas distintas, por omisión
c (test$statistic, test$parameter, pv=test$p.value)
esta      <- (medias[1] - medias[2]) / sqrt (varianzas[1]/n + varianzas[2]/m)
c         <- gl (varianzas)                  #grados de libertad corregidos
p_valor   <- 2 * (1 - pt (abs(esta), c))     #p-valor del contraste
c (estadístico=esta, gl=c, pv=p_valor)
## Corrección de los grados de libertad y relación entre las varianzas
var_x <- 1                            
var_y <- seq (0.1, 30, 0.1)
VAR   <- cbind (var_x, var_y)           #matriz de varianzas de X y Y
ces   <- apply(VAR, 1, gl)              #valores teóricos de c
razón <- var_y / var_x
plot (sqrt (razón), ces,
      type='l', xlab=expression(sigma[Y]/sigma[X]), ylab="grados de libertad")
abline (h = n+m-2, col='red')           #g.l. con varianzas iguales
abline (v = 1, col='red')               #igualdad de varianzas
## con varianzas muestrales parecidas la corrección penaliza mucho los g.l.

## (%i1) assume (n>1, m>1)$

## (%i2) gl(ratio) := (1/n+ratio/m)^2 / 
##     (1/(n^3-n^2)+ratio^2/(m^3-m^2))$

## (%i3) solve (diff (gl (ratio), ratio), ratio);

##                                    2
##                                   m  - m            m
## (%o3)                    [ratio = ------, ratio = - -]
##                                    2                n
##                                   n  - n

## (%i4) gl ((m^2-m) / (n^2-n)), ratsimp;

## (%o4)                              n + m - 2