Adiciona scripts de Monte Carlo.

scripts/ce089-07.R

+#-----------------------------------------------------------------------
+# 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
+    # wilcox.test(x, y)$p.value
+}
+simul()
+
+# P-valores do teste sob H0.
+p <- replicate(1000, simul())
+
+# Distribuição do teste é uniforme quando é exato.
+plot(ecdf(p))
+abline(a = 0, b = 1, col = 2)
+
+#-----------------------------------------------------------------------
+# O teste F tem distribuição exata?
+
+g <- gl(5, 4)
+
+simul <- 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())
+
+# Distribuição do teste é uniforme quando é exato.
+plot(ecdf(p))
+abline(a = 0, b = 1, col = 2)
+
+#-----------------------------------------------------------------------
+# O teste da deviance tem distribuição exata chi-quadrado?
+
+g <- gl(2, 5)
+
+simul <- 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())
+
+# Distribuição do teste é uniforme quando é exato.
+plot(ecdf(p))
+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)
+y <- X %*% beta + rnorm(x, 0, 0.3)
+
+plot(y ~ x)
+
+m0 <- lm(y ~ x)
+confint(m0)
+
+simul <- function() {
+    y <- X %*% beta + rnorm(x, 0, 0.3)
+    m0 <- lm(y ~ x)
+    confint(m0)
+}
+
+ic <- replicate(100, simul())
+str(ic)
+
+# Com i = 1 é o intercepto e i = 2 é a taxa.
+i <- 1
+plot(NULL, NULL, xlim = c(0, dim(ic)[3]), ylim = range(ic[i, , ]))
+segments(1:ncol(ic[i, , ]),
+         ic[i, 1, ],
+         1:ncol(ic[i, , ]),
+         ic[i, 2, ])
+abline(h = beta[i])
+
+# Aumentando o número de simulações
+ic <- replicate(1000, simul())
+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)
+
+#-----------------------------------------------------------------------
+# Taxa de cobertura em modelo não linear.
+
+A <- 10
+V <- 2
+x <- 0:10
+y <- A * x/(V + x) + rnorm(x, 0, 0.3)
+
+plot(y ~ x)
+curve(A * x/(V + x), add = TRUE)
+
+n0 <- nls(y ~ A * x/(V + x), start = list(A = 10, V = 2))
+
+# Intervalo de Wald.
+confint.default(n0)
+
+# Limite superior. Note que IC de Wald é est + z * stderror.
+summary(n0)$coefficients[, 1] +
+              qnorm(0.975) * summary(n0)$coefficients[, 2]
+
+# Intervalo de perfil de verossimilhança.
+confint(n0)
+
+par(mfrow = c(1, 2))
+plot(profile(n0))
+layout(1)
+
+# x <- 0:10
+# X <- cbind(1, x)
+# beta <- c(3, 0.5)
+# y <- X %*% beta + rnorm(x, 0, 0.3)
+#
+# n0 <- nls(y ~ b0 + b1 * x, start = list(b0 = 3, b1 = 0.5))
+# confint.default(n0)
+# confint(n0)
+#
+# summary(n0)$coefficients[, 1] +
+#               qnorm(0.975) * summary(n0)$coefficients[, 2]
+#
+# par(mfrow = c(1, 2))
+# plot(profile(n0))
+# layout(1)
+
+simul <- function() {
+    y <- A * x/(V + x) + rnorm(x, 0, 0.3)
+    n0 <- nls(y ~ A * x/(V + x),
+              start = list(A = 10, V = 2))
+    confint.default(n0)
+    # confint(n0)
+}
+
+ic <- replicate(100, simul())
+str(ic)
+
+# Com i = 1 é o intercepto e i = 2 é a taxa.
+i <- 1
+plot(NULL, NULL, xlim = c(0, dim(ic)[3]), ylim = range(ic[i, , ]))
+segments(1:ncol(ic[i, , ]),
+         ic[i, 1, ],
+         1:ncol(ic[i, , ]),
+         ic[i, 2, ])
+abline(h = c(A, V)[i])
+
+# Aumentando para ter uma estimativa mais precisa.
+n <- 1000
+ic <- replicate(n, simul())
+str(ic)
+
+# Taxas de cobertura.
+sum(ic[1, 1, ] < A & ic[1, 2, ] > A)/n
+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.
+
+#-----------------------------------------------------------------------

scripts/ce089-08.R

+#-----------------------------------------------------------------------
+# 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))
+xy <- locator(n = 20, type = "p", pch = 19)
+# dput(lapply(xy, round, digits = 3))
+
+xy <- structure(list(x = c(0.204, 0.186, 0.529, 0.529, 0.385, 0.579,
+                           0.918, 0.798, 0.634, 0.761, 0.704, 0.485,
+                           0.291, 0.341, 0.402, 0.833, 0.972, 0.625,
+                           0.732, 0.315),
+                     y = c(0.829, 0.545, 0.526, 0.752, 0.674, 0.648,
+                           0.792, 0.33, 0.121, 0.127, 0.332, 0.352,
+                           0.188, 0.452, 0.221, 0.524, 0.957, 0.964,
+                           0.755, 0.332)),
+                .Names = c("x", "y"))
+
+#-----------------------------------------------------------------------
+
+# 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).
+d <- dist(xy)
+d
+
+# Obtém a menor distância.
+m <- min(d)
+
+# 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, {
+    min(dist(cbind(x = runif(n), y = runif(n))))
+})
+
+# 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.
+
+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.863, 0.64, 0.337, 0.088, 0.077, 0.346,
+                           0.654, 0.619),
+                     y = c(0.92, 0.922, 0.916, 0.935, 0.737, 0.674,
+                           0.67, 0.665, 0.452, 0.502, 0.461, 0.454,
+                           0.256, 0.26, 0.26, 0.219, 0.045, 0.06, 0.058,
+                           0.439)),
+                .Names = c("x", "y"))
+
+# 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))
+
+# Todas as distância entre os pontos (a estatística de teste é um vetor
+# e não um escalar).
+d <- c(dist(xy))
+d
+
+# Gerando estatísticas de teste por simulação.
+D <- replicate(9999, {
+    dist(cbind(x = runif(n), y = runif(n)))
+})
+
+# Juntanto observado com simulado.
+D <- cbind(d, D)
+str(D)
+
+plot(ecdf(D[, 1]), cex = 0)
+for (i in 2:ncol(D)) {
+    lines(ecdf(D[, i]), cex = 0)
+}
+lines(ecdf(D[, 1]), cex = 0, col = 2, lwd = 2)
+
+#-----------------------------------------------------------------------
+# Amostragem por Conjuntos Ordenados.
+
+# Extraí valores de uma população normal padrão.
+rand <- function(n) {
+    rnorm(n, 0, 1)
+}
+rand(10)
+
+# Ranked Set Sampling (Perfect Ranking). Calcula a média.
+rss <- function(m = 3) {
+    x <- matrix(rand(m * m), nrow = m)
+    x <- apply(x, MARGIN = 1, sort)
+    mean(diag(x))
+}
+rss(m = 5)
+
+# Simple Random Sampling. Calcula a média.
+srs <- function(m = 3) {
+    mean(rand(m))
+}
+rss(m = 5)
+
+# Ranked Set Sampling with Unperfect Ranking. Calcula a média.
+rssur <- function(m = 3, rho) {
+    x <- MASS::mvrnorm(m * m,
+                       mu = c(0, 0),
+                       Sigma = matrix(c(1, rho, rho, 1), nrow = 2))
+    g <- gl(m, m)
+    b <- by(data = x,
+            INDICES = g,
+            FUN = function(y) {
+                o <- order(y[, 2])
+                # cbind(y[o, 1], y[o, 2])
+                y[o, 1]
+            }, simplify = TRUE)
+    mean(diag(do.call(rbind, b)))
+}
+rssur(m = 5, rho = 0.95)
+
+#-----------------------------------------------------------------------
+# Simulação.
+
+# Tamanho da amostra.
+m <- 10
+
+# Distribuição da média na amostra aleatória simples.
+system.time(m1 <- replicate(1000, srs(m = m)))
+
+# Distribuição da média na amostra por conjuntos ordenados perfeito.
+system.time(m2 <- replicate(1000, rss(m = m)))
+
+# Distribuição da média na amostra por conjuntos ordenados imperfeito.
+system.time(m3 <- replicate(1000, rssur(m = m, rho = 0.75)))
+
+library(lattice)
+library(latticeExtra)
+
+densityplot(~m1 + m2 + m3,
+            auto.key = TRUE,
+            layout = c(NA, 1)) +
+    layer(panel.abline(v = 0, lty = 2))
+
+ecdfplot(~m1 + m2 + m3,
+         auto.key = TRUE,
+         layout = c(NA, 1))
+
+qqmath(~m1 + m2 + m3,
+       auto.key = TRUE,
+       layout = c(NA, 1))
+
+#-----------------------------------------------------------------------