#
# El problema de Behrens y Fisher.
#
rm(list=ls())
#
# Datos simulados.
#
n1 <- 5
n2 <- 5
mu1 <- 0
mu2 <- 5
delta <- mu1-mu2
sigma21 <- 100
sigma22 <- 100
alpha <- 0.05
#
# Usar simulaciones para calcular la verdadera cobertura de un intervalo frecuentista de 100(1-alpha)%. 
#
simuls <- 100000
nfuera <- 0
for (i in 1:simuls){
  y1 <- rnorm(n1,mu1,sqrt(sigma21))
  y2 <- rnorm(n2,mu2,sqrt(sigma22))
  ci <- t.test(y1,y2)$conf.int
  if (delta<ci[1]|delta>ci[2]){
    nfuera <- nfuera+1
  } 
}
trueconf <- 1-nfuera/simuls
#
# Si n1 y n2 son pequeñas, la aproximación es algo imprecisa.  Para n1, n2 más grandes, funciona algo mejor.
#
# Datos reales: alturas en pulgas de estudiantes de dos clases de estadística para biologia tomados del 
# Handbook of Biological Statistics por John Macdonald.
#
datos <- ("grupo valor
        2pm    69   
        2pm    70   
        2pm    66   
        2pm    63   
        2pm    68   
        2pm    70   
        2pm    69   
        2pm    67   
        2pm    62   
        2pm    63   
        2pm    76   
        2pm    59   
        2pm    62   
        2pm    62   
        2pm    75   
        2pm    62   
        2pm    72    
        2pm    63   
        5pm    68
        5pm    62
        5pm    67
        5pm    68
        5pm    69
        5pm    67
        5pm    61
        5pm    59
        5pm    62
        5pm    61
        5pm    69
        5pm    66
        5pm    62
        5pm    62
        5pm    61
        5pm    70
        ")
datos <- read.table(textConnection(datos),header=TRUE)
grupo <- datos$grupo
valor <- datos$valor
y1 <- valor[grupo=="2pm"]
y2 <- valor[grupo=="5pm"]
par(mfrow=c(1,2))
hist(y1,main="Alturas (en pulgas): grupo 2pm",xlab="altura",ylab="frecuencia")
hist(y2,main="Alturas (en pulgas): grupo 5pm",xlab="altura",ylab="frecuencia")
par(mfrow=c(1,1))
#
# Intervalo frecuentista.
#
alpha=0.05
ci <- t.test(y1,y2,conf.level=1-alpha)$conf.int
c(ci[1],ci[2])
#
# Calcular los estadísticos suficientes.
#
n1 <- length(y1)
media1 <- mean(y1)
dt1 <- sd(y1)
n2 <- length(y2)
media2 <- mean(y2)
dt2 <- sd(y2)
dtcomb <- sqrt((dt1^2)/n1+(dt2^2)/n2) # El estimador de la varianza
phi <- atan(dt1*sqrt(n2)/(dt2*sqrt(n1))) # El ángulo del BF.
#
# Inferencia bayesiana "exacta".  Usamos la biblioteca asht.
#
if (!require("asht")) install.packages("asht")
library(asht)
tq <- qbf(c(alpha/2,1-alpha/2),n1,n2,R=phi)
cibayes <- media1-media2+dtcomb*tq
cibayes
#
# Ahora hacer la inferencia a través del MonteCarlo.
#
simuls <- 100000
t1 <- rt(simuls,n1-1)
t2 <- rt(simuls,n2-1)
tdif <- t1*sin(phi)-t2*cos(phi)
tq <- quantile(tdif,c(alpha/2,1-alpha/2))
cimc <- media1-media2+dtcomb*tq
cimc