##===============================================================================
## Título: Inferência Estatística (IE) - Teste de precisão
## Curso : Introdução à estatística
## Autor : José Cláudio Faria/UESC/DCET
## Data  : 04/10/2019 14:41:17
## Versão: v3
##===============================================================================
## 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 unilateral a direita
#.. 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___

#... A condicional abaixo visa automatizar a IE no caso de alteração das amostras
#... garantindo que o teste é sempre unilateral a direita
#... H1: Var(I) > Var(J)
if(var(A) > var(B)) {
  n1 <- length(A)-1
  n2 <- length(B)-1
  (Fcal <- var(A)/var(B))
} else {
  n1 <- length(B)-1
  n2 <- length(A)-1
  (Fcal <- var(B)/var(A))
  }

# Valor limite da distribuição F para aceitação de Ho
(Flim <- qf(p=alfa,
            df1=n1,
            df2=n2,
            lower=FALSE))

# Para imprimir na tela do console
if(var(A) > var(B)) {
  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')
} else {
  cat('var(B) = ', var(B), '\n')
  cat('var(A) = ', var(A), '\n\n')
  cat('Fcal = var(B)/var(A) =', round(var(B)/var(A), 2), '\n')
}

# Conclusão da inferência
if(Fcal >= Flim) {
  cat('Rejeita H0\n')
} else {
  cat('Aceita H0\n')}

# Generalização do teste usando a função var.test
if(var(A) > var(B)) {
  (var.test(A,
            B,
            alt='greater'))
} else {
  (var.test(B,
            A,
            alt='greater'))}

#-----------------------------------------------------------------------------
# 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 <- max(Fcal, Flim)                

max_y <- max(df(seq(0, 
                    Flim, 
                    length=10),
                n1, 
                n2))

plot(1,
     main= 'Teste unilateral a direita',
     xlim=c(0, 3*max_x),
     ylim=c(0, 1.4*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=Flim,
         y0=0,
         x1=Flim,
         y1=max_y,
         col='red',
         lty=3)

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

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

# Hipóteses
text(x=Flim,
     y=c(max_y, .95*max_y),
     pos=4,
     offset=2,
     col='darkblue',
     label=c(ifelse(var(A) > var(B),
                    expression(sigma^2*(A) > sigma^2*(B)),
                    expression(sigma^2*(B) > sigma^2*(A))),
             ifelse(var(A) > var(B),
                    expression(sigma^2*(A) <= sigma^2*(B)),
                    expression(sigma^2*(B) <= sigma^2*(A)))))

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

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

# Achuras sob a curva
xa1 <- seq(0,
           Flim,
           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),
        y=c(0,
            ya1,
            0),
        col=adjustcolor('darkgreen', alpha.f=.1),
        border='darkgreen')

xa2 <- seq(Flim,
           20,
           len=1e3)

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

polygon(x=c(Flim,
            xa2),
        y=c(0,
            ya2),
        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=''))