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))

#-----------------------------------------------------------------------