diff --git a/scripts/ce089-07.R b/scripts/ce089-07.R
index d91dbf43b9ab191c721dd797b90ad83f24e4ca5b..a5534011a1fdd8d5d0fd65ba5ca488144d86d977 100644
--- a/scripts/ce089-07.R
+++ b/scripts/ce089-07.R
@@ -1,17 +1,32 @@
+#=======================================================================
+# Avaliação da distribuição sob H_0 de testes de hipótese.
+
#-----------------------------------------------------------------------
# Distribuição do p-valor em um teste cuja estatÃstica de teste tem
# distribuição exata.
-simul <- function() {
- x <- rnorm(20)
- y <- rnorm(20)
- t.test(x, y, var.equal = TRUE)$p.value
+simul_diff_medias <- function(dist = "normal", n = 20) {
+ amos <- switch(dist,
+ "normal" = {
+ list(x = rnorm(n),
+ y = rnorm(n))
+ },
+ "exponencial" = {
+ list(x = rexp(n, rate = 2),
+ y = rexp(n, rate = 2))
+ },
+ "poisson" = {
+ list(x = rpois(n, lambda = 4),
+ y = rpois(n, lambda = 4))
+ })
+ t.test(amos$x, amos$y, var.equal = TRUE)$p.value
# wilcox.test(x, y)$p.value
}
-simul()
+
+simul_diff_medias()
# P-valores do teste sob H0.
-p <- replicate(1000, simul())
+p <- replicate(1000, simul_diff_medias())
# Distribuição do teste é uniforme quando é exato.
plot(ecdf(p))
@@ -22,15 +37,14 @@ abline(a = 0, b = 1, col = 2)
g <- gl(5, 4)
-simul <- function() {
+simul_f_anova <- function() {
y <- rnorm(g)
- # plot(y ~ as.integer(g))
m0 <- lm(y ~ g)
anova(m0)[1, 5]
}
# P-valores do teste sob H0.
-p <- replicate(1000, simul())
+p <- replicate(1000, simul_f_anova())
# Distribuição do teste é uniforme quando é exato.
plot(ecdf(p))
@@ -41,23 +55,25 @@ abline(a = 0, b = 1, col = 2)
g <- gl(2, 5)
-simul <- function() {
+simul_chi_deviance <- function() {
y <- rpois(g, lambda = 5)
m0 <- glm(y ~ g, family = poisson)
anova(m0, test = "Chi")[2, 5]
}
# P-valores do teste sob H0.
-p <- replicate(1000, simul())
+p <- replicate(1000, simul_chi_deviance())
# Distribuição do teste é uniforme quando é exato.
-plot(ecdf(p))
+plot(ecdf(p), cex = 0.1)
abline(a = 0, b = 1, col = 2)
-#-----------------------------------------------------------------------
+#=======================================================================
# Cobertura de intervalos de confiança.
+#-----------------------------------------------------------------------
# Gerando dados de um modelo de regressão linear simples.
+
x <- 0:10
X <- cbind(1, x)
beta <- c(3, 0.5)
@@ -68,13 +84,13 @@ plot(y ~ x)
m0 <- lm(y ~ x)
confint(m0)
-simul <- function() {
+simul_confint <- function() {
y <- X %*% beta + rnorm(x, 0, 0.3)
m0 <- lm(y ~ x)
confint(m0)
}
-ic <- replicate(100, simul())
+ic <- replicate(100, simul_confint())
str(ic)
# Com i = 1 é o intercepto e i = 2 é a taxa.
@@ -87,12 +103,12 @@ segments(1:ncol(ic[i, , ]),
abline(h = beta[i])
# Aumentando o número de simulações
-ic <- replicate(1000, simul())
+ic <- replicate(1000, simul_confint())
str(ic)
-# Conta quantos IC contem o valor paramétrico.
-sum(ic[1, 1, ] < 3 & ic[1, 2, ] > 3)
-sum(ic[2, 1, ] < 0.5 & ic[2, 2, ] > 0.5)
+# Conta quantos IC contem o valor paramétrico real.
+sum(ic[1, 1, ] < 3 & ic[1, 2, ] > 3)/dim(ic)[3]
+sum(ic[2, 1, ] < 0.5 & ic[2, 2, ] > 0.5)/dim(ic)[3]
#-----------------------------------------------------------------------
# Taxa de cobertura em modelo não linear.
@@ -137,7 +153,7 @@ layout(1)
# plot(profile(n0))
# layout(1)
-simul <- function() {
+simul_confint_nls <- function() {
y <- A * x/(V + x) + rnorm(x, 0, 0.3)
n0 <- nls(y ~ A * x/(V + x),
start = list(A = 10, V = 2))
@@ -145,7 +161,7 @@ simul <- function() {
# confint(n0)
}
-ic <- replicate(100, simul())
+ic <- replicate(100, simul_confint_nls())
str(ic)
# Com i = 1 é o intercepto e i = 2 é a taxa.
@@ -159,7 +175,7 @@ abline(h = c(A, V)[i])
# Aumentando para ter uma estimativa mais precisa.
n <- 1000
-ic <- replicate(n, simul())
+ic <- replicate(n, simul_confint_nls())
str(ic)
# Taxas de cobertura.
@@ -169,4 +185,103 @@ sum(ic[2, 1, ] < V & ic[2, 2, ] > V)/n
# TODO: substitua o IC de Wald pelo IC de perfil de verossimilhança e
# repita tudo. Tire suas próprias conclusões.
+#=======================================================================
+# Curva de poder de teste.
+
+#-----------------------------------------------------------------------
+# Teste t versus teste de Wilcox.
+
+# Várias opções de distribuição para comparar com fuga das suposições.
+simul_dados <- function(mu1 = 5, mu2 = mu1, sigma = 1, n = 10, dist = "normal") {
+ amos <- switch(dist,
+ "normal" = {
+ list(y1 = rnorm(n, mu1, sigma),
+ y2 = rnorm(n, mu2, sigma))
+ },
+ "exponencial" = {
+ list(y1 = rexp(n, rate = mu1),
+ y2 = rexp(n, rate = mu2))
+ },
+ "uniforme" = {
+ list(y1 = runif(n, min = mu1, max = mu1 + 1),
+ y2 = runif(n, min = mu2, max = mu2 + 1))
+ },
+ "poisson" = {
+ list(y1 = rpois(n, lambda = mu1),
+ y2 = rpois(n, lambda = mu2))
+ })
+ return(amos)
+}
+
+simul_t <- function(y1, y2) {
+ t.test(y1, y2, var.equal = TRUE)$p.value
+}
+
+simul_wilcox <- function(y1, y2) {
+ wilcox.test(y1, y2)$p.value
+}
+
+d_grid <- seq(0, 3, length.out = 30)
+p_grid <- sapply(d_grid,
+ FUN = function(d) {
+ alpha <- 0.05
+ N <- 1000
+ r <- replicate(N, {
+ # Avaliação pareada do poder.
+ dts <- simul_dados(mu1 = 5, mu2 = 5 + d)
+ c(simul_t(dts$y1, dts$y2) < alpha,
+ simul_wilcox(dts$y1, dts$y2) < alpha)
+ }, simplify = TRUE)
+ rowSums(r)/N
+ })
+
+p_grid <- t(p_grid)
+
+matplot(d_grid,
+ p_grid,
+ xlab = expression("Diferença entre médias:" ~ mu[1] - mu[2]),
+ ylab = expression(Pr("Rejeitar" ~ H[0] ~ "|" ~ mu[1] - mu[2])),
+ type = "l",
+ ylim = c(0, 1))
+abline(h = 0.05, col = 2)
+legend("right",
+ legend = c("Teste t", "Teste de Wilcox"),
+ lty = 1:2,
+ col = 1:2,
+ bty = "n")
+grid()
+
+#--------------------------------------------
+# Usando distribuição uniforme (fuga de pressuposto pro teste t).
+
+d_grid <- seq(0, 1.2, length.out = 30)
+p_grid <- sapply(d_grid,
+ FUN = function(d) {
+ alpha <- 0.05
+ N <- 1000
+ r <- replicate(N, {
+ # Avaliação pareada do poder.
+ dts <- simul_dados(n = 6, mu1 = 5, mu2 = 5 + d, dist = "uniforme")
+ c(simul_t(dts$y1, dts$y2) < alpha,
+ simul_wilcox(dts$y1, dts$y2) < alpha)
+ }, simplify = TRUE)
+ rowSums(r)/N
+ })
+
+p_grid <- t(p_grid)
+
+matplot(d_grid,
+ p_grid,
+ xlab = expression("Diferença entre médias:" ~ mu[1] - mu[2]),
+ ylab = expression(Pr("Rejeitar" ~ H[0] ~ "|" ~ mu[1] - mu[2])),
+ type = "l",
+ ylim = c(0, 1))
+abline(h = 0.05, col = 2)
+legend("right",
+ legend = c("Teste t", "Teste de Wilcox"),
+ lty = 1:2,
+ col = 1:2,
+ bty = "n")
+grid()
+
#-----------------------------------------------------------------------
diff --git a/scripts/ce089-08.R b/scripts/ce089-08.R
index 8906759a8e44dfcc1e17cd7461be9cfacd0628b6..9f914386d934cecc1acf216a4c95386aae3aed32 100644
--- a/scripts/ce089-08.R
+++ b/scripts/ce089-08.R
@@ -1,8 +1,18 @@
+#=======================================================================
+# Testes de hipótese Monte Carlo.
+
#-----------------------------------------------------------------------
# Teste para independência de processo pontual.
-plot(NULL, NULL, xlim = c(0, 1), ylim = c(0, 1), asp = 1)
-lines(x = c(0, 1, 1, 0, 0), y = c(0, 0, 1, 1, 0))
+plot(xy,
+ NULL,
+ xlim = c(0, 1),
+ ylim = c(0, 1),
+ asp = 1,
+ axes = FALSE,
+ ann = FALSE)
+lines(x = c(0, 1, 1, 0, 0),
+ y = c(0, 0, 1, 1, 0))
# xy <- locator(n = 20, type = "p", pch = 19)
# dput(lapply(xy, round, digits = 3))
@@ -21,10 +31,17 @@ xy <- structure(list(x = c(0.204, 0.186, 0.529, 0.529, 0.385, 0.579,
# Transforma a lista em matriz.
xy <- do.call(cbind, xy)
-plot(xy, NULL, xlim = c(0, 1), ylim = c(0, 1), asp = 1)
-lines(x = c(0, 1, 1, 0, 0), y = c(0, 0, 1, 1, 0))
-
-# Calcula todas as distâncias (euclidianas).
+plot(xy,
+ NULL,
+ xlim = c(0, 1),
+ ylim = c(0, 1),
+ asp = 1,
+ axes = FALSE,
+ ann = FALSE)
+lines(x = c(0, 1, 1, 0, 0),
+ y = c(0, 0, 1, 1, 0))
+
+# Calcula todas as distâncias (euclidianas) entre pontos.
d <- dist(xy)
d
@@ -52,7 +69,7 @@ abline(v = m, col = 2)
2 * sum(M > m)/length(M)
#-----------------------------------------------------------------------
-# Moficando a estatÃstica de teste.
+# ATTENTION: dados um tanto patológicos.
xy <- structure(list(x = c(0.088, 0.326, 0.577, 0.846, 0.857, 0.568,
0.306, 0.077, 0.077, 0.328, 0.597, 0.883,
@@ -67,8 +84,51 @@ xy <- structure(list(x = c(0.088, 0.326, 0.577, 0.846, 0.857, 0.568,
# Transforma a lista em matriz.
xy <- do.call(cbind, xy)
-plot(xy, NULL, xlim = c(0, 1), ylim = c(0, 1), asp = 1)
-lines(x = c(0, 1, 1, 0, 0), y = c(0, 0, 1, 1, 0))
+plot(xy,
+ NULL,
+ xlim = c(0, 1),
+ ylim = c(0, 1),
+ asp = 1,
+ axes = FALSE,
+ ann = FALSE)
+lines(x = c(0, 1, 1, 0, 0),
+ y = c(0, 0, 1, 1, 0))
+
+#-----------------------------------------------------------------------
+# Modificando a estatÃstica de teste 1.
+
+# A soma das menores distâncias entre vizinhos mais próximos.
+d <- dist(xy)
+D <- as.matrix(d)
+diag(D) <- Inf
+apply(D, MARGIN = 2, FUN = min)
+m <- sum(apply(D, MARGIN = 2, FUN = min))
+
+# Número de pontos.
+n <- nrow(xy)
+
+# Geração de estatÃsticas por simulação do modelo assumido sob
+# hipótese nula.
+M <- replicate(9999, {
+ d <- dist(cbind(x = runif(n), y = runif(n)))
+ D <- as.matrix(d)
+ diag(D) <- Inf
+ sum(apply(D, MARGIN = 2, FUN = min))
+})
+
+# Concatena as estatÃsticas simuladas com a observada.
+M <- c(m, M)
+
+# Gráfico da distribuição acumulada empÃrica.
+plot(ecdf(M))
+abline(h = c(0.025, 0.975), lty = 2)
+abline(v = m, col = 2)
+
+# P-valor.
+2 * sum(M > m)/length(M)
+
+#-----------------------------------------------------------------------
+# Moficando a estatÃstica de teste 2.
# Todas as distância entre os pontos (a estatÃstica de teste é um vetor
# e não um escalar).
@@ -90,6 +150,9 @@ for (i in 2:ncol(D)) {
}
lines(ecdf(D[, 1]), cex = 0, col = 2, lwd = 2)
+#=======================================================================
+# Avaliação de métodos de amostragem ou delineamentos.
+
#-----------------------------------------------------------------------
# Amostragem por Conjuntos Ordenados.
@@ -162,3 +225,23 @@ qqmath(~m1 + m2 + m3,
layout = c(NA, 1))
#-----------------------------------------------------------------------
+
+# Criando um grid de rho.
+rho_grid <- seq(0, 0.95, by = 0.05)
+
+res <- lapply(rho_grid,
+ FUN = function(rho) {
+ replicate(500, rssur(m = 10, rho = rho))
+ })
+str(res)
+
+n <- length(res[[1]])
+res <- do.call(c, res)
+
+res <- data.frame(rho = rep(rho_grid, each = n), valor = res)
+str(res)
+
+densityplot(~valor | rho, data = res) +
+ layer(panel.abline(v = 0, lty = 2))
+
+#-----------------------------------------------------------------------