Commit 64080ace authored by Walmes Marques Zeviani's avatar Walmes Marques Zeviani

Merge branch 'master' of gitlab.c3sl.ufpr.br:walmes/ce089

parents f477c7fa e80a0a18
#=======================================================================
# 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 do p-valor em um teste cuja estatística de teste tem
# distribuição exata. # distribuição exata.
simul <- function() { simul_diff_medias <- function(dist = "normal", n = 20) {
x <- rnorm(20) amos <- switch(dist,
y <- rnorm(20) "normal" = {
t.test(x, y, var.equal = TRUE)$p.value 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 # wilcox.test(x, y)$p.value
} }
simul()
simul_diff_medias()
# P-valores do teste sob H0. # P-valores do teste sob H0.
p <- replicate(1000, simul()) p <- replicate(1000, simul_diff_medias())
# Distribuição do teste é uniforme quando é exato. # Distribuição do teste é uniforme quando é exato.
plot(ecdf(p)) plot(ecdf(p))
...@@ -22,15 +37,14 @@ abline(a = 0, b = 1, col = 2) ...@@ -22,15 +37,14 @@ abline(a = 0, b = 1, col = 2)
g <- gl(5, 4) g <- gl(5, 4)
simul <- function() { simul_f_anova <- function() {
y <- rnorm(g) y <- rnorm(g)
# plot(y ~ as.integer(g))
m0 <- lm(y ~ g) m0 <- lm(y ~ g)
anova(m0)[1, 5] anova(m0)[1, 5]
} }
# P-valores do teste sob H0. # P-valores do teste sob H0.
p <- replicate(1000, simul()) p <- replicate(1000, simul_f_anova())
# Distribuição do teste é uniforme quando é exato. # Distribuição do teste é uniforme quando é exato.
plot(ecdf(p)) plot(ecdf(p))
...@@ -41,23 +55,25 @@ abline(a = 0, b = 1, col = 2) ...@@ -41,23 +55,25 @@ abline(a = 0, b = 1, col = 2)
g <- gl(2, 5) g <- gl(2, 5)
simul <- function() { simul_chi_deviance <- function() {
y <- rpois(g, lambda = 5) y <- rpois(g, lambda = 5)
m0 <- glm(y ~ g, family = poisson) m0 <- glm(y ~ g, family = poisson)
anova(m0, test = "Chi")[2, 5] anova(m0, test = "Chi")[2, 5]
} }
# P-valores do teste sob H0. # P-valores do teste sob H0.
p <- replicate(1000, simul()) p <- replicate(1000, simul_chi_deviance())
# Distribuição do teste é uniforme quando é exato. # Distribuição do teste é uniforme quando é exato.
plot(ecdf(p)) plot(ecdf(p), cex = 0.1)
abline(a = 0, b = 1, col = 2) abline(a = 0, b = 1, col = 2)
#----------------------------------------------------------------------- #=======================================================================
# Cobertura de intervalos de confiança. # Cobertura de intervalos de confiança.
#-----------------------------------------------------------------------
# Gerando dados de um modelo de regressão linear simples. # Gerando dados de um modelo de regressão linear simples.
x <- 0:10 x <- 0:10
X <- cbind(1, x) X <- cbind(1, x)
beta <- c(3, 0.5) beta <- c(3, 0.5)
...@@ -68,13 +84,13 @@ plot(y ~ x) ...@@ -68,13 +84,13 @@ plot(y ~ x)
m0 <- lm(y ~ x) m0 <- lm(y ~ x)
confint(m0) confint(m0)
simul <- function() { simul_confint <- function() {
y <- X %*% beta + rnorm(x, 0, 0.3) y <- X %*% beta + rnorm(x, 0, 0.3)
m0 <- lm(y ~ x) m0 <- lm(y ~ x)
confint(m0) confint(m0)
} }
ic <- replicate(100, simul()) ic <- replicate(100, simul_confint())
str(ic) str(ic)
# Com i = 1 é o intercepto e i = 2 é a taxa. # Com i = 1 é o intercepto e i = 2 é a taxa.
...@@ -87,12 +103,12 @@ segments(1:ncol(ic[i, , ]), ...@@ -87,12 +103,12 @@ segments(1:ncol(ic[i, , ]),
abline(h = beta[i]) abline(h = beta[i])
# Aumentando o número de simulações # Aumentando o número de simulações
ic <- replicate(1000, simul()) ic <- replicate(1000, simul_confint())
str(ic) str(ic)
# Conta quantos IC contem o valor paramétrico. # Conta quantos IC contem o valor paramétrico real.
sum(ic[1, 1, ] < 3 & ic[1, 2, ] > 3) sum(ic[1, 1, ] < 3 & ic[1, 2, ] > 3)/dim(ic)[3]
sum(ic[2, 1, ] < 0.5 & ic[2, 2, ] > 0.5) sum(ic[2, 1, ] < 0.5 & ic[2, 2, ] > 0.5)/dim(ic)[3]
#----------------------------------------------------------------------- #-----------------------------------------------------------------------
# Taxa de cobertura em modelo não linear. # Taxa de cobertura em modelo não linear.
...@@ -137,7 +153,7 @@ layout(1) ...@@ -137,7 +153,7 @@ layout(1)
# plot(profile(n0)) # plot(profile(n0))
# layout(1) # layout(1)
simul <- function() { simul_confint_nls <- function() {
y <- A * x/(V + x) + rnorm(x, 0, 0.3) y <- A * x/(V + x) + rnorm(x, 0, 0.3)
n0 <- nls(y ~ A * x/(V + x), n0 <- nls(y ~ A * x/(V + x),
start = list(A = 10, V = 2)) start = list(A = 10, V = 2))
...@@ -145,7 +161,7 @@ simul <- function() { ...@@ -145,7 +161,7 @@ simul <- function() {
# confint(n0) # confint(n0)
} }
ic <- replicate(100, simul()) ic <- replicate(100, simul_confint_nls())
str(ic) str(ic)
# Com i = 1 é o intercepto e i = 2 é a taxa. # Com i = 1 é o intercepto e i = 2 é a taxa.
...@@ -159,7 +175,7 @@ abline(h = c(A, V)[i]) ...@@ -159,7 +175,7 @@ abline(h = c(A, V)[i])
# Aumentando para ter uma estimativa mais precisa. # Aumentando para ter uma estimativa mais precisa.
n <- 1000 n <- 1000
ic <- replicate(n, simul()) ic <- replicate(n, simul_confint_nls())
str(ic) str(ic)
# Taxas de cobertura. # Taxas de cobertura.
...@@ -169,4 +185,103 @@ sum(ic[2, 1, ] < V & ic[2, 2, ] > V)/n ...@@ -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 # TODO: substitua o IC de Wald pelo IC de perfil de verossimilhança e
# repita tudo. Tire suas próprias conclusões. # 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()
#----------------------------------------------------------------------- #-----------------------------------------------------------------------
#=======================================================================
# Testes de hipótese Monte Carlo.
#----------------------------------------------------------------------- #-----------------------------------------------------------------------
# Teste para independência de processo pontual. # Teste para independência de processo pontual.
plot(NULL, NULL, xlim = c(0, 1), ylim = c(0, 1), asp = 1) plot(xy,
lines(x = c(0, 1, 1, 0, 0), y = c(0, 0, 1, 1, 0)) 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) # xy <- locator(n = 20, type = "p", pch = 19)
# dput(lapply(xy, round, digits = 3)) # 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, ...@@ -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. # Transforma a lista em matriz.
xy <- do.call(cbind, xy) xy <- do.call(cbind, xy)
plot(xy, NULL, xlim = c(0, 1), ylim = c(0, 1), asp = 1) plot(xy,
lines(x = c(0, 1, 1, 0, 0), y = c(0, 0, 1, 1, 0)) NULL,
xlim = c(0, 1),
# Calcula todas as distâncias (euclidianas). 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 <- dist(xy)
d d
...@@ -52,7 +69,7 @@ abline(v = m, col = 2) ...@@ -52,7 +69,7 @@ abline(v = m, col = 2)
2 * sum(M > m)/length(M) 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, 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, 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, ...@@ -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. # Transforma a lista em matriz.
xy <- do.call(cbind, xy) xy <- do.call(cbind, xy)
plot(xy, NULL, xlim = c(0, 1), ylim = c(0, 1), asp = 1) plot(xy,
lines(x = c(0, 1, 1, 0, 0), y = c(0, 0, 1, 1, 0)) 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 # Todas as distância entre os pontos (a estatística de teste é um vetor
# e não um escalar). # e não um escalar).
...@@ -90,6 +150,9 @@ for (i in 2:ncol(D)) { ...@@ -90,6 +150,9 @@ for (i in 2:ncol(D)) {
} }
lines(ecdf(D[, 1]), cex = 0, col = 2, lwd = 2) lines(ecdf(D[, 1]), cex = 0, col = 2, lwd = 2)
#=======================================================================
# Avaliação de métodos de amostragem ou delineamentos.
#----------------------------------------------------------------------- #-----------------------------------------------------------------------
# Amostragem por Conjuntos Ordenados. # Amostragem por Conjuntos Ordenados.
...@@ -162,3 +225,23 @@ qqmath(~m1 + m2 + m3, ...@@ -162,3 +225,23 @@ qqmath(~m1 + m2 + m3,
layout = c(NA, 1)) 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))
#-----------------------------------------------------------------------
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment