#========================================================================================== # Aula 16 da disciplina ce223 (13/05/2011) # Estatística básica, construção de funções com opções de controle # Professor Walmes M. Zeviani # www.leg.ufpr.br/ce223 #========================================================================================== #------------------------------------------------------------------------------------------ # aprimorando a função para teste de hipótese sobre a média incluindo o tipo de teste # bilateral, unilateral à esquerda, unilateral à direita my.test <- function(x, y, alpha, tipo){ ## x, é o vetor com dados ## y, é o vetor ou escalar ## alpha, é o nível de significância do teste ## é o tipo de teste: bilateral, à esquerda, à direita if(length(y)==1){ # vai fazer o teste para uma média h0 <- y m <- mean(x); v <- var(x); n <- length(x); e <- sqrt(v/n) tc <- (m-h0)/e if(tipo=="bilateral"){ # vai aplicar o teste bilateral, P(xbar!=h0) tt <- -pt(alpha/2, df=n-1) dec <- ifelse(abs(tc)>tt, "rejeito", "aceito") } if(tipo=="esquerda"){ # vai aplicar o teste unilateral à esquerda, P(xbartt, "rejeito", "aceito") } if(tipo=="direita"){ # vai aplicar o teste unilateral à direita, P(xbar>h0) tt <- -pt(alpha, df=n-1) dec <- ifelse(tc>tt, "rejeito", "aceito") } return(list(stat=c(media=m, variancia=v, erropad=e, tcalc=tc, ttab=tt), n=n, decisao=dec, tipo=tipo)) } if(length(y)>1){ # vai fazer o teste para diferença de duas médias a1 <- x; a2 <- y; nivel <- alpha m1 <- mean(a1); m2 <- mean(a2) n1 <- length(a1); n2 <- length(a2) v1 <- var(a1); v2 <- var(a2) v <- ((n1-1)*v1+(n2-1)*v2)/(n1+n2-2) dif <- abs(m2-m1) tt <- -qt((1-nivel)/2, df=n1+n2-2) ep <- sqrt(v*(1/n1+1/n2)) tc <- (dif-0)/ep return(list(medias=c(m1=m1, m2=m2, dif=dif), variancias=c(v1=v1, v2=v2, v=v), n=c(n1=n1, n2=n2, gl=n1+n2-2), statistica=c(tc=tc, tt=tt))) } } #------------------------------------------------------------------------------------------ y <- rbeta(100, 2, 3) z <- rbeta(80, 2.5, 3) my.test(y, h0=0.5, alpha=0.05, tipo="bilateral") my.test(y, h0=0.5, alpha=0.05, tipo="esquerda") my.test(y, h0=0.5, alpha=0.05, tipo="direita") my.test(y, z, alpha=0.05, tipo="direita") my.test(y, 0.5, alpha=0.05, tipo="direita") #------------------------------------------------------------------------------------------ help(t.test, help_type="html") #------------------------------------------------------------------------------------------ x1 <- rgamma(20, 1, 0.9) x2 <- rgamma(22, 1.1, 0.9) t.test(x1, x2, var.equal=TRUE) dados <- data.frame(amostra=rep(c(1,2), c(length(x1), length(x2))), valor=c(x1, x2)) dados dados <- dados[sample(1:nrow(dados)),] dados subset(dados, amostra==1) with(dados, t.test(valor~amostra, var=TRUE)) #------------------------------------------------------------------------------------------ dados <- data.frame(peso=rnorm(100, 80, 4), altura=rnorm(100, 1.78, 0.2), idade=rpois(100, lambda=22), renda=rexp(100, rate=0.01)) str(dados) aux <- apply(dados, 2, t.test) str(aux) aux$peso$conf lapply(aux, function(slot){ c(slot$conf) }) sapply(aux, function(slot){ c(slot$conf) }) #------------------------------------------------------------------------------------------