Universidade Federal do Paraná
Curso de Estatística
CE 083 -
Estatística Computacional I - 2014/2
Prof. Dr. Walmes Marques Zeviani
##-----------------------------------------------------------------------------
## Introdução aos testes de hipótese. O curioso caso da moeda humana.
## x <- scan()
## dput(x)
## Sequência de caras de coroas de um voluntário (amostra).
x <- c(0, 0, 1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 1, 0, 1, 0, 0, 1, 1, 1,
1, 1, 0, 0, 1, 0, 1, 1, 1, 1, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0, 0,
1, 1, 0, 1, 0, 1, 1, 0, 1, 0, 0, 0, 1, 1, 0, 0, 1, 1, 0)
## Número de trocas (estatística).
t <- sum(abs(diff(x))); t
## Número de trocas para o lançamento de uma moeda verdadeira (H0).
sum(abs(diff(sample(0:1, size=60, repl=TRUE))))
## Comportamento de uma moeda verdadeira (distribuição sob H0).
empdis <- replicate(10000,
sum(abs(diff(sample(0:1, size=60, repl=TRUE)))))
tb <- table(empdis)
freq <- c(tb)
xval <- as.numeric(names(tb))
## Gráfico da distribuição empírica sob H0.
plot(freq~xval, type="h")
abline(v=t+0.2, col=2)
## Probabilidade de uma medida tão ou mais extrema que a atual (teste
## unilateral).
sum(freq[xval>=32])/length(empdis)
## Distribuição teórica existe para o número de trocas, é binomial.
xb <- 0:59
pr <- dbinom(xb, size=59, prob=0.5)
## Empírica e teórica ambas sob H0.
plot(freq/length(empdis)~xval, type="h")
points(xb-0.2, pr, type="h", col=4)
abline(v=t+0.2, col=2)
##-----------------------------------------------------------------------------
## Mais um exemplo: detectar um padrão não uniforme na distribuição de
## pontos.
plot(xlim=c(0,1), ylim=c(0,1), x=NULL, y=NULL, asp=1)
abline(v=0:1, h=0:1, lty=2)
## L <- locator(n=20, type="p")
## dput(lapply(L, round, digits=4))
L <-
structure(list(x = c(0.157, 0.4668, 0.2835, 0.2333, 0.8029, 0.661,
0.1744, 0.5017, 0.1788, 0.0806, 0.5563, 0.8771, 0.3315, 0.7898,
0.4101, 0.6283, 0.7527, 0.6763, 0.8858, 0.4843), y = c(0.141,
0.2588, 0.5423, 0.2501, 0.6273, 0.9457, 0.9523, 0.7277, 0.7582,
0.5248, 0.4616, 0.0756, 0.0538, 0.322, 0.4289, 0.6513, 0.4747,
0.2239, 0.8716, 0.8978)), .Names = c("x", "y"))
points(L)
## x <- runif(20)
## y <- runif(20)
## plot(xlim=c(0,1), ylim=c(0,1), x=x, y=y, asp=1)
## abline(v=0:1, h=0:1, lty=2)
##
## min(dist(cbind(x,y)))
## obs <- min(dist(cbind(L$x, L$y)))
r <- replicate(10000,
{
x <- runif(20)
y <- runif(20)
min(dist(cbind(x,y)))
})
## Estatística observada na amostra.
obs <- min(dist(cbind(L$x, L$y)))
## Distribuição empírica sob H0 e a estística calculada.
plot(density(r), xlim=c(0, 1.4*obs))
abline(v=obs, col=2)
## P-valor.
prop.table(table(r>=obs))
## Em tracejado, o valor para um teste com 5% de significância.
plot(density(r), xlim=c(0, 1.4*obs))
abline(v=obs, col=2)
abline(v=quantile(r, prob=0.95), lty=2)
##-----------------------------------------------------------------------------
## Curva de poder em função da diferença entre médias (d).
## x <- rnorm(10, 1.8, 0.2)
## t <- t.test(x=x, mu=1.8)
## str(t)
delta <- seq(0, 0.8, by=0.02)
res <- sapply(delta,
function(d){
replicate(1000, {
x <- rnorm(6, 1.8+d, 0.2)
t <- t.test(x=x, mu=1.8)
t$p.value
})
})
p <- apply(res, MARGIN=2, function(x) sum(x<0.05)/nrow(res))
plot(p~delta, type="l", ylim=c(0,1),
ylab="Probabilidade de rejeitar H0",
xlab="Difereça de médias")
##-----------------------------------------------------------------------------
## Em função do tamanho da amostra. Diferença fixa em 0.2.
N <- seq(3, 25, by=1)
res <- sapply(N,
function(n){
replicate(1000, {
x <- rnorm(n, 1.8+0.2, 0.2)
t <- t.test(x=x, mu=1.8)
t$p.value
})
})
p <- apply(res, MARGIN=2, function(x) sum(x<0.05)/nrow(res))
## plot(p~N, type="l", ylim=c(0,1))
## Qual o tamanho mínimo da amostra para ter uma probabilidade de
## rejeitar H0 superior à 95% quando a popução tiver uma diferença de
## médias com relação à H0 de no mínimo 0.2?
ap <- approx(x=p, y=N, xou=0.95); ap
plot(p~N, type="l", ylim=c(0,1))
abline(h=ap$x, v=ap$y, lty=2)
