#===============================================================================
#. Título: Inferência Estatística (IE) - Teste de precisão
#. Curso : Introdução à estatística
#. Autor : José Cláudio Faria/UESC/DCET
# Data  : 01/09/2021 12:34:32
#===============================================================================
# Objetivos gerais
#===============================================================================
# a) Apresentar os recursos básicos do R de apoio a IE;
# b) Fundamentar o estudante em inferência básica usando a distribuição F
#    no teste de precisão a partir de amostras de dois processos: A e B.
#===============================================================================
# O script permite que se altere, livremente:
#   - os valores das amostras
#   - os número de elementos das amostras
#   - o erro adotado na inferência
#===============================================================================
#. Hipóteses para teste bilateral
#.. H0: Var(I)  = Var(J)
#.. H1: Var(I) != Var(J)
#===============================================================================

#.. ___Ini opções___
#.. Amostras
#... Informadas
#A <- c(10.2, 8.7, 9.5, 12.0, 9.0, 11.2, 12.5, 10.9, 8.9, 10.6)
#B <- c(9.9, 9.2, 10.4, 10.5, 11.0, 11.3, 9.6, 9.4, 10.0, 10.4)

#... Geradas
A <- rnorm(15,
           mean=10,
           sd=2)

B <- rnorm(15,
           mean=10,
           sd=2)

#... Inferência: erro tipo I (alfa) = 5%
alfa <- 5/100
#.. ___Fim opções___

#... O teste é bilateral, não importa a ordem das estimativas das variâncias
#... H1: Var(I) != Var(J)
n1 <- length(A)-1
n2 <- length(B)-1

(Fcal <- var(A)/var(B))

# Valores limite da distribuição F para aceitação de Ho
(Flim_i <- qf(p=alfa/2,
              df1=n1,
              df2=n2,
              lower=TRUE))

(Flim_s <- qf(p=alfa/2,
              df1=n1,
              df2=n2,
              lower=FALSE))

# Para imprimir na tela do console
cat('var(A) = ', var(A), '\n')
cat('var(B) = ', var(B), '\n\n')
cat('Fcal = var(A)/var(B) = ', round(var(A)/var(B), 2), '\n')

# Conclusão da inferência
ifelse(Fcal > Flim_i & Fcal < Flim_s, 
       res <- ('Aceita H0\n'), 
       res <- ('Rejeita H0\n'))  
cat(res)

# Generalização do teste usando a função var.test
(var.test(A,
          B,
          alt='two.side'))

#-----------------------------------------------------------------------------
# Adorno (opcional)
#-----------------------------------------------------------------------------
# Abrindo um novo dispositivo gráfico e preparando para receber a curva da
# função de densidade de probabiliddes F e as condições da inferência.
# Prepara o dispositivo gráfico
max_x <- min(Fcal, Flim_s)                

max_y <- max(df(seq(0, 
                    2, 
                    length=20),
                n1, 
                n2))

plot(1,
     main= 'Teste bilateral',
     xlim=c(-3, 10),
     ylim=c(0, 1.2*max_y), 
     type='n',
     xlab='F',
     ylab='f(F)')

# Plota a curva específica f(F) relativa ao tamanho das amostras
curve(df(x,
         n1,
         n2),
      n=1e3,
      col='darkblue',
      add=TRUE,
      lwd=2)

# Linha decisão
segments(x0=c(Flim_i, Flim_s),
         y0=0,
         x1=c(Flim_i, Flim_s),
         y1=max_y,
         col='red',
         lty=3)

# Valor do FLim
text(x=c(Flim_i, Flim_s),
     y=.3*max_y,
     pos=c(2, 4),
     col='darkblue',
     label=paste('Ftab',
                 '(',
                 n1,
                 ', ',
                 n2,
                 ')',
                 ' = ',      
                 round(c(Flim_i, Flim_s), 2),
                 sep=''))

# Texto das hipóteses
text(x=c(Flim_i, (Flim_i+Flim_s)/2, Flim_s),
     y=.6*max_y,
     labels=c('RRH0',
              'RAH0',
              'RRH0'),
     col=c('red',
           'darkgreen', 
           'red'),
     pos=c(2, 3, 4),
     offset=3)

# Hipóteses
text(x=Flim_s,
     y=c(max_y, .95*max_y),
     pos=4,
     offset=2,
     col='darkblue',
     label=c('H0:',
             'H1:'))
text(x=Flim_s,
     y=c(max_y, .95*max_y),
     pos=4,
     offset=3.8,
     col='darkblue',
     label=c(expression(sigma^2*(A) == sigma^2*(B)),
             expression(sigma^2*(A) != sigma^2*(B))))

# Texto do valor de prova
text(x=Flim_s,
     y=c(.8*max_y, .75*max_y),
     pos=4,
     offset=2,
     col='darkblue',
     label= c(expression(Fcal == s^2*(A) / s^2*(B)),
              substitute(Fcal == val, list(val = round(Fcal, 2)))))

# Texto do erro
text(x=c(Flim_i, Flim_s),
     y=.1*max_y,
     pos=c(2, 4),
     offset=1,
     col='red',
     substitute(alpha == paste(val, '%'),
                list(val = round(100*alfa/2, 2))))

# Achuras sob a curva
# Erro Superior
xa1 <- seq(0,
           Flim_i,
           len=1e3)

ya1 <- df(xa1,
          df1=n1,
          df2=n2)

# Para evitar Inf no caso da amostra do numerador ser de tamanho 2
ya1 <- sub(Inf,
           1e2,
           ya1)

polygon(x=c(0,
            xa1,
            Flim_i),
        y=c(0,
            ya1,
            0),
        col=adjustcolor('red', alpha.f=.1),
        border='red')

# Aceitação
xa2 <- seq(Flim_i,
           Flim_s,
           len=1e3)

ya2 <- df(xa2,
          df1=n1,
          df2=n2)

polygon(x=c(Flim_i,
            xa2,
            Flim_s),
        y=c(0,
            ya2,
            0),
        col=adjustcolor('darkgreen', alpha.f=.2),
        border='darkgreen')

# Erro Superior
xa3 <- seq(Flim_s,
           20,
           len=1e3)

ya3 <- df(xa3,
          df1=n1,
          df2=n2)

polygon(x=c(Flim_s,
            xa3),
        y=c(0,
            ya3),
        col=adjustcolor('red', alpha.f=.2),
        border='red')

# Ponto Fcal
points(x=Fcal,
       y=0,
       pch=18,
       col='red',
       bg='red',
       cex=1.2)

# Texto Fcal
text(x=Fcal,
     y=.01,
     col='red',
     pos=1,
     cex=.8,
     label=paste('Fcal = ',
                 round(Fcal, 2),
                 sep=''))