##=============================================================================== ## 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=''))