#=======================================================================
# Testes de hipótese Monte Carlo.

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

# x <- scan()
# dput(x)

# Amostra observada.
x <- c(c(1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 0),
       c(1, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 1, 1, 1, 0, 1, 0, 1, 0, 1, 1, 0))

# Caracterização.
n <- length(x)
p <- mean(x)
d <- sum(abs(diff(x)))

D <- replicate(9999, {
    # Amostras sob a H_0.
    r <- rbinom(n = n, size = 1, prob = 0.5)
    # Estatística calculada.
    sum(abs(diff(r)))
})

D <- c(d, D)

plot(ecdf(D))
abline(v = d, col = 2)
abline(h = c(0.025, 0.975), col = 2, lty = 2)

2 * sum(D > d)/length(D)

#-----------------------------------------------------------------------
# Teste para independência de processo pontual.

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

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

# 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) # Faz sentido esse teste ser bilateral?
sum(M > m)/length(M)

#-----------------------------------------------------------------------
# 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,
                           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,
     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)
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).
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)

#=======================================================================
# Avaliação de métodos de amostragem ou delineamentos.

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

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

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

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