Universidade Federal do Paraná
Curso de Estatística
CE 089 -
Estatística Computacional II - 2014/2
Prof. Dr. Walmes Marques Zeviani
Tabela de conteúdo
##-----------------------------------------------------------------------------
## Obtém o valor da estatística t do teste de Student para a média de
## uma população. Assume que a distribuição de X seja normal.
simula0 <- function(n, mu0, sig0){
X <- rnorm(n, mean=mu0, sd=sig0)
T <- (mean(X)-mu0)/(sqrt(var(X)/n))
return(T)
}
simula0(n=10, mu0=0, sig0=1)
## [1] 1.522
t <- replicate(10000, simula0(n=10, mu0=0, sig0=1))
## Comparação por distribuições acumuladas.
plot(ecdf(t), xlim=c(-5, 5))
curve(pt(x, df=10-1), add=TRUE, col=2)
curve(pnorm(x), add=TRUE, col=3)
## Comparação pela densidade.
plot(density(t), xlim=c(-5, 5))
curve(dt(x, df=10-1), add=TRUE, col=2)
curve(dnorm(x), add=TRUE, col=3)
## p-valor da simulação.
sum(abs(t)>=qt(0.025, df=10-1, lower=FALSE))/length(t)
## [1] 0.0507
##-----------------------------------------------------------------------------
## Distribuição da estatística com afastamento dos pressupostos sobre a
## distribuição da população (X) que não tem distribuição normal.
simula1 <- function(n, mu0=1){
## X <- runif(n, -0.5, 0.5)+mu0
X <- rexp(n, 1)
T <- (mean(X)-mu0)/(sqrt(var(X)/n))
return(T)
}
n <- 5
t <- replicate(10000, simula1(n=n))
plot(ecdf(t), xlim=c(-5, 5))
curve(pt(x, df=n-1), add=TRUE, col=2)
curve(pnorm(x), add=TRUE, col=3)
## p-valor real vs nível se significância nominal.
sum(abs(t)>=qt(0.025, df=n-1, lower=FALSE))/length(t)
## [1] 0.1143
##-----------------------------------------------------------------------------
## Função que obtém o valor da estatística t para H0: mu==mu0 simulando
## de uma população normal com média mu.
simula2 <- function(n, sig0, mu0, mu){
X <- rnorm(n, mu, sig0)
t <- t.test(X, mu=mu0)
return(t$p.value)
}
## Sequência de valores para mu.
mu <- seq(0, 2, length.out=30)
n <- 10
N <- 1000
## Probabilidade de rejeitar H0 quando ela é falsa.
rej <- sapply(mu,
function(mu){
r <- replicate(N, simula2(n=n, sig0=1, mu0=0, mu))
sum(r<0.05)/N
})
## plot(rej~mu, type="o")
##-----------------------------------------------------------------------------
## Consultando com a função analítica para cálculo do poder do teste t.
d <- 1
pw <- power.t.test(n=n, delta=d, sd=1, sig.level=0.05,
type="one.sample")
str(pw)
## List of 8
## $ n : num 10
## $ delta : num 1
## $ sd : num 1
## $ sig.level : num 0.05
## $ power : num 0.803
## $ alternative: chr "two.sided"
## $ note : NULL
## $ method : chr "One-sample t test power calculation"
## - attr(*, "class")= chr "power.htest"
plot(rej~mu, type="o")
abline(v=d, h=pw$power, lty=2)
##-----------------------------------------------------------------------------
## Obtendo a curva de poder analítica.
my.ptt <- function(delta){
sapply(delta,
function(d){
pw <- power.t.test(n=n, delta=d, sd=1,
sig.level=0.05,
type="one.sample",
alternative="two.sided")
pw$power
})
}
muu <- c(-rev(mu), mu)
pw <- my.ptt(muu)
## Curva de poder teórica e empírica (vinda da simulação).
plot(pw~muu, type="l", ylab="Poder", xlab=expression(mu))
lines(rej~mu, type="o", col=2)
## Possivelmente tem um problema com a função analítica pois no ponto
## zero retorna 0.025 de poder quando deveria ser 0.05 por definição.
## Nem todo teste de hipótese tem uma função para cálculo de poder. Além
## do mais, quandos os pressupostos não são atendidos, não se tem como
## obter a curva de poder se não por simulação de dados sob cada cenário
## de afastamento de pressuposto.
apropos("power")
## [1] "power" "power.anova.test" "power.prop.test" "power.t.test"