diff --git a/scripts/metodos-intensivos.Rmd b/scripts/metodos-intensivos.Rmd
new file mode 100644
index 0000000000000000000000000000000000000000..3b29cd036655b7cbc75a8354d9559976884ca581
--- /dev/null
+++ b/scripts/metodos-intensivos.Rmd
@@ -0,0 +1,853 @@
+---
+title: "Métodos computacionalmente intensivos para inferência"
+author: "[Walmes Marques Zeviani](http://lattes.cnpq.br/4410617539281650)"
+date: '`r Sys.Date()`'
+output:
+ html_document:
+ theme: yeti
+---
+
+```{r, include = FALSE}
+opts_chunk$set(cache = TRUE,
+ tidy = FALSE,
+ fig.width = 7,
+ fig.height = 6,
+ fig.align = "center",
+ eval.after= "fig.cap",
+ warning = FALSE,
+ error = FALSE,
+ message = FALSE)
+```
+
+# Testes de aleatorização
+
+## Uma senhora toma chá
+
+```{r}
+#-----------------------------------------------------------------------
+# Uma senhora toma chá.
+
+# Quantidade de maneiras de gerar dois grupos de 4 xícaras usando 8.
+n_poss <- choose(8, 4)
+n_poss
+
+# De 4 xícaras selecionadas para um grupo, X representa o número de
+# acertos.
+n_acertar <- sapply(4:0,
+ FUN = function(i) {
+ # (formas de acertar x em 4) * (formas de errar x em 4).
+ choose(4, i) * choose(4, 4 - i)
+ })
+n_acertar
+
+# Pr(Acerto total) = 1/70 < 5%.
+cumsum(n_acertar)/n_poss
+
+# Comprovando por simulação.
+v <- rep(0:1, each = 4)
+mean(replicate(n = 100000,
+ expr = all(sample(v) == v)))
+
+# Essa forma de fazer a conta é mais realista mas dá na mesma.
+mean(replicate(n = 100000,
+ expr = all(sample(v) == sample(v))))
+```
+
+## Teste para a diferença de médias
+
+```{r}
+#-----------------------------------------------------------------------
+# Exemplo com teste para a diferença de médias.
+
+# Comprimentos de crânios de cães pré-históricos.
+m <- c(120, 107, 110, 116, 114, 111, 113, 117, 114, 112)
+f <- c(110, 111, 107, 108, 110, 105, 107, 106, 111, 111)
+
+plot(ecdf(m), xlim = range(c(m, f)), col = "cyan")
+lines(ecdf(f), col = "magenta")
+rug(m, col = "cyan")
+rug(f, col = "magenta")
+
+# Diferença de média observada.
+d <- mean(m) - mean(f)
+d
+
+# Todos as combinações possíveis.
+choose(n = 20, k = 10)
+
+#--------------------------------------------
+# Com todas as combinações possíveis (exaustão).
+
+# Para construir todas as combinações possíveis.
+k <- combn(x = 1:20, m = 10)
+dim(k)
+
+# Vetor com os valores dos dois grupos concatenados.
+mf <- c(m, f)
+
+# Vetor cheio de zeros.
+g <- integer(20)
+
+# Calcula a diferença para todo arranjo possível.
+D <- apply(k,
+ MARGIN = 2,
+ FUN = function(i) {
+ # Alguns elementos do vetor convertidos para 1.
+ g[i] <- 1L
+ # Cálculo da estatística de teste.
+ -diff(tapply(mf, g, FUN = mean))
+ })
+
+# Histograma da distribuição da estatística sib H_0.
+hist(D, col = "gray50")
+rug(D)
+abline(v = d, col = 2)
+
+plot(ecdf(D), cex = 0)
+rug(D)
+abline(v = d, col = 2)
+
+# P-valor do teste.
+2 * sum(D >= d)/length(D)
+
+#--------------------------------------------
+# Com simulação (não necessáriamente exaustivo).
+
+# Variáveis que indentifica os grupos.
+g <- rep(1:2, each = 10)
+
+# Apenas para conferir.
+cbind(g, mf)
+
+# Replicando a diferença para grupos formados por aleatorização.
+D <- replicate(9999, {
+ gg <- sample(g)
+ -diff(tapply(mf, gg, FUN = mean))
+})
+
+# Tem que juntar a estatística observada com as simuladas.
+D <- c(D, d)
+
+hist(D, col = "gray50")
+rug(D)
+abline(v = d, col = 2)
+
+plot(ecdf(D), cex = 0)
+rug(D)
+abline(v = d, col = 2)
+
+# P-valor do teste.
+2 * sum(D >= d)/length(D)
+```
+
+## Teste para a correlação
+
+```{r}
+#-----------------------------------------------------------------------
+# Teste de aleatorização para a correlação.
+
+# N = 5 para um par de medidas.
+x <- c(4.1, 8.3, 2.9, 10.8, 9.5)
+y <- c(3.7, 5.1, 1.0, 7.7, 8.9)
+
+cbind(x, y)
+
+plot(x, y)
+
+# Estatística de teste na amostra original.
+r0 <- cor(x, y, method = "pearson")
+r0
+
+# Todas as permutações possiveis: 5! = 120 do vetor x.
+X <- gtools::permutations(n = length(x), r = length(x), v = x)
+str(X)
+head(X)
+
+# As 120 correlações obtidas para cada arranjo.
+r <- apply(X, MARGIN = 1, FUN = cor, y = y, method = "spearman")
+
+# P-valor do teste.
+2 * sum(r >= r0)/length(r)
+```
+
+## Índice de Moran
+
+ * <http://gis.stackexchange.com/questions/161887/significance-test-for-morans-i-using-monte-carlo-simulation>.
+ * <https://en.wikipedia.org/wiki/Moran's_I>.
+ * <http://rstudio-pubs-static.s3.amazonaws.com/268058_dd37b98f25a4496b9f0a7eb2fcf892cd.html>.
+ * <http://rspatial.org/analysis/rst/3-spauto.html>.
+
+```{r}
+#-----------------------------------------------------------------------
+# Índice de Moran para medir dependência espacial.
+
+# Coordenadas dos eventos em uma malha regular.
+x <- 1:8
+y <- 1:8
+
+# Para criar a matriz de pesos.
+ind <- expand.grid(i = 1:length(x), j = 1:length(y))
+f <- function(i, j) {
+ u <- min(3, sum(abs(ind[i, ] - ind[j, ])))
+ c(0, 1, sqrt(1/2), 0)[u + 1]
+}
+w <- matrix(0, nrow(ind), nrow(ind))
+for (i in 1:nrow(ind)) {
+ for (j in 1:nrow(ind)) {
+ w[i, j] <- f(i, j)
+ }
+}
+w <- w/sum(w)
+
+image(w)
+
+# Índice de Moran.
+moran <- function(x, weights) {
+ # Tamanho da amostra.
+ n <- length(x)
+ # Valores normalizados.
+ z <- as.vector((x - mean(x))/sd(x))
+ # Índice de Moran.
+ as.vector(z %*% weights %*% (z * sqrt(n/(n - 1))))
+}
+
+# Teste de permutação com gráfico.
+ppt <- function(z, w, N = 10000, stat, ...) {
+ # Índice de Moran por reamostragem.
+ sim <- replicate(N,
+ moran(sample(z), w))
+ # Determina o p-valor.
+ p.value <- mean((all <- c(stat, sim)) >= stat)
+ # Histograma da distribuição empírica sob H_0.
+ hist(sim,
+ sub = paste("p =", round(p.value, 4)),
+ xlim = range(all),
+ ...)
+ abline(v = stat, col = "#903030", lty = 3, lwd = 2)
+ return(p.value)
+}
+
+# Observações simuladas.
+set.seed(17)
+par(mfrow = c(2, 3))
+
+# Dados com dependência espacial ---------------------------------------
+# Indução de autocorrelação por meio de uma tendência.
+z <- matrix(rexp(length(x) * length(y),
+ outer(x, y^2)),
+ length(x))
+image(log(z), main = "Com dependência")
+
+# Índice de Moran com dados originais.
+stat <- moran(z, w)
+stat
+
+hist(z)
+ppt(z, w, stat = stat, main = "I de Moran", xlab = "I")
+
+# Dados sem dependência espacial ---------------------------------------
+# Geração de de um conjunto de dados sob hipótese nula.
+z <- matrix(rnorm(length(x) * length(y), 0, 1/2), length(x))
+image(z, main = "Sem dependência")
+
+# Índice de Moran com dados originais.
+stat <- moran(z, w)
+stat
+
+hist(z)
+ppt(z, w, stat = stat, main = "I de Moran", xlab = "I")
+```
+
+# Jackknife
+
+## Regressão não linear
+
+```{r}
+#-----------------------------------------------------------------------
+# Aplicação de Jackknife em modelos de regressão não linear.
+
+library(latticeExtra)
+library(car)
+library(alr3)
+
+# Curva "palpite".
+start <- list(th0 = 75, th1 = 0.5, th2 = 50)
+xyplot(C ~ Temp, data = segreg) +
+ layer(panel.curve(th0 + th1 * (x - th2) * (x >= th2) +
+ 0 * (x < th2), lty = 2),
+ data = start)
+
+# Ajuste.
+n0 <- nls(C ~ th0 + th1 * (Temp - th2) * (Temp >= th2) +
+ 0 * (Temp < th2),
+ data = segreg,
+ start = start)
+
+# Estimativas e medidas de ajuste.
+summary(n0)
+
+# Observados e preditos.
+xyplot(C ~ Temp, data = segreg) +
+ layer(panel.curve(th0 + th1 * (x - th2) * (x >= th2) +
+ 0 * (x < th2)),
+ data = as.list(coef(n0)))
+
+#-----------------------------------------------------------------------
+# Estimativas Jackknife para os parâmetros.
+
+
+theta <- coef(n0) # Estimativas com todos os dados.
+n <- nrow(segreg) # Número de observações.
+
+# Matriz vazia para armazenar os pseudo-valores.
+thetaj <- matrix(0, nrow = n, ncol = length(theta))
+colnames(thetaj) <- names(theta)
+
+# Ajustes deixando uma observação de fora (leave-one-out).
+for (i in 1:n) {
+ # Reajusta o modelo, i.e. estima com leave-one-out.
+ n1 <- nls(C ~ th0 + th1 * (Temp - th2) * (Temp >= th2) +
+ 0 * (Temp < th2),
+ data = segreg[-i, ],
+ start = coef(n0))
+ # Calcula os pseudo valores.
+ thetaj[i, ] <- n * theta - (n - 1) * coef(n1)
+}
+
+# Os primeiros pseudo valores.
+head(thetaj)
+
+# Estimativas pontuais.
+jk <- colMeans(thetaj)
+cbind(MLE = theta, Jack = jk, Vício = theta - jk)
+
+# Erros padrões.
+cbind(MLE = summary(n0)$coefficients[, "Std. Error"],
+ Jack = apply(thetaj, MARGIN = 2, sd)/sqrt(n))
+
+# Pacote com função pronta para Jackknife para modelos não lineares.
+library(nlstools)
+# ls("package:nlstools")
+
+# Aplica Jackknife sobre o modelo.
+j0 <- nlsJack(n0)
+summary(j0)
+```
+
+## Dose letal para 50% da população
+
+```{r}
+#-----------------------------------------------------------------------
+# Inferência para a DL 50.
+
+# library(labestData)
+# data(PaulaTb3.12)
+# str(PaulaTb3.12)
+# help(PaulaTb3.12, h = "html")
+# # Converte a resposta para binário.
+# PaulaTb3.12$resp <- as.integer(PaulaTb3.12$resp) - 1
+# dput(PaulaTb3.12)
+
+PaulaTb3.12 <-
+structure(list(vol = c(3.7, 3.5, 1.25, 0.75, 0.8, 0.7, 0.6, 1.1,
+0.9, 0.9, 0.8, 0.55, 0.6, 1.4, 0.75, 2.3, 3.2, 0.85, 1.7, 1.8,
+0.4, 0.95, 1.35, 1.5, 1.6, 0.6, 1.8, 0.95, 1.9, 1.6, 2.7, 2.35,
+1.1, 1.1, 1.2, 0.8, 0.95, 0.75, 1.3), razao = c(0.825, 1.09,
+2.5, 1.5, 3.2, 3.5, 0.75, 1.7, 0.75, 0.45, 0.57, 2.75, 3, 2.33,
+3.75, 1.64, 1.6, 1.415, 1.06, 1.8, 2, 1.36, 1.35, 1.36, 1.78,
+1.5, 1.5, 1.9, 0.95, 0.4, 0.75, 0.03, 1.83, 2.2, 2, 3.33, 1.9,
+1.9, 1.625), resp = c(1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0,
+1, 1, 1, 1, 1, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 1, 0, 0, 1,
+1, 1, 0, 0, 1)), .Names = c("vol", "razao", "resp"), row.names = c(NA,
+39L), class = "data.frame")
+
+layout(1)
+plot(resp ~ vol, data = PaulaTb3.12)
+
+m0 <- glm(resp ~ vol,
+ data = PaulaTb3.12,
+ family = binomial)
+summary(m0)
+
+# DL_50.
+dl <- -coef(m0)[1]/coef(m0)[2]
+
+# str(m0$family)$link
+
+plot(resp ~ vol, data = PaulaTb3.12)
+curve(m0$family$linkinv(coef(m0)[1] + coef(m0)[2] * x),
+ add = TRUE)
+abline(v = dl, h = 0.5, lty = 2)
+
+n <- nrow(PaulaTb3.12)
+i <- 1:n
+
+# Estimativas por Jackknife.
+DL <- sapply(i,
+ FUN = function(i) {
+ m0 <- glm(resp ~ vol,
+ data = PaulaTb3.12[-i, ],
+ family = binomial)
+ # Estimativa Parcial.
+ DL <- -coef(m0)[1]/coef(m0)[2]
+ # Pseudo-valor.
+ n * dl - (n - 1) * DL
+ })
+
+m <- mean(DL)
+e <- sd(DL)/sqrt(n)
+
+cbind(Est = rbind(MLE = dl,
+ Jack = m),
+ StdErr = rbind(MLE = car::deltaMethod(m0, g = "-(Intercept)/vol")["SE"],
+ Jack = e))
+
+plot(ecdf(DL))
+rug(DL)
+
+plot(density(DL))
+rug(DL)
+
+# NOTE: alguma pista sobre porque a distribuição fica bimodal?
+```
+
+## Uma situação problemática
+
+```{r}
+#-----------------------------------------------------------------------
+# Aplicação na Mediana (derivada não suave) -> problema.
+
+# Ordena o vetor de valores.
+x <- sort(precip)
+
+# Calcula a mediana com os dados originais.
+M <- median(x)
+M
+
+n <- length(x)
+i <- 1:length(x)
+
+# Estimativas da mediana por Jackknife.
+y <- sapply(i,
+ FUN = function(i) {
+ ep <- median(x[-i])
+ pv <- n * M - (n - 1) * ep
+ return(pv)
+ })
+
+stem(y)
+
+# NOTE: Explique porque o Jackknife para a mediana só retornou 2
+# valores? Considere que o tamanho da amostra é dado abaixo.
+length(x)
+```
+
+# Bootstrap
+
+## Regressão não linear
+
+```{r}
+#-----------------------------------------------------------------------
+# Aplicação de bootstrap. Estudo de caso com modelos não lineares.
+
+library(alr3)
+str(turk0)
+
+# Gráfico dos valores observados.
+plot(Gain ~ A, data = turk0)
+
+# Lista com os valores iniciais.
+start <- list(Int = 625, Ass = 180, Mei = 0.1)
+
+# Adiciona as curvas.
+with(start, {
+ curve(Int + Ass * A/(Mei + A),
+ xname = "A", add = TRUE, col = 2)
+ curve(Int + Ass * (1 - 2^(-A/Mei)),
+ xname = "A", add = TRUE, col = 4)})
+
+# Ajuste do primeiro modelo.
+n0 <- nls(Gain ~ Int + Ass * A/(Mei + A),
+ data = turk0,
+ start = start)
+
+# Ajuste do segundo modelo.
+n1 <- nls(Gain ~ Int + Ass * (1 - 2^(-A/Mei)),
+ data = turk0,
+ start = start)
+
+cbind(coef(n0), coef(n1))
+
+#-----------------------------------------------------------------------
+# Criando funções que fazem reamostragem dos dados e ajustam o modelo.
+
+n <- nrow(turk0)
+i <- 1:n
+
+myboot0 <- function(formula, start) {
+ n0 <- nls(formula = formula,
+ data = turk0[sample(i, n, replace = TRUE), ],
+ start = start)
+ coef(n0)
+}
+
+#-----------------------------------------------------------------------
+# Replicando.
+
+# B = 999 reamostragens.
+b0 <- replicate(999,
+ myboot0(Gain ~ Int + Ass * A/(Mei + A), start))
+b1 <- replicate(999,
+ myboot0(Gain ~ Int + Ass * (1 - 2^(-A/Mei)), start))
+
+# Conjunto de 1000 estimativas bootstrap.
+b0 <- cbind(b0, coef(n0))
+b1 <- cbind(b1, coef(n1))
+str(b0)
+
+# Calcula o vício absoluto (bootstrap - original).
+rbind(n0 = rowMeans(b0) - coef(n0),
+ n1 = rowMeans(b1) - coef(n1))
+
+# Variância do estimador.
+rbind(n0 = apply(b0, MARGIN = 1, var),
+ n1 = apply(b1, MARGIN = 1, var))
+
+# Função para calcular o erro quadrático médio.
+eqm <- function(x) {
+ sum((x - tail(x, 1))^2)/length(x)
+}
+
+# Erro quadrático médio.
+rbind(n0 = apply(b0, 1, FUN = eqm),
+ n1 = apply(b1, 1, FUN = eqm))
+```
+
+## Inferência para valores preditos
+
+```{r}
+#-----------------------------------------------------------------------
+# Curvas ajustadas de onde é possível determinar uma banda de confiança.
+
+bts <- cbind(rep(1:2, each = ncol(b0)),
+ as.data.frame(rbind(t(b0), t(b1))))
+splom(bts[, -1], groups = bts[, 1])
+
+par(mfrow = c(1, 2))
+plot(Gain ~ A, data = turk0, type = "n")
+apply(b0,
+ MARGIN = 2,
+ FUN = function(b) {
+ curve(b[1] + b[2] * x/(b[3] + x),
+ add = TRUE, col = rgb(0, 1, 1, 0.25), n = 51)
+ invisible()
+ })
+points(Gain ~ A, data = turk0)
+abline(v = 0.2, lty = 2)
+plot(Gain ~ A, data = turk0, type = "n")
+apply(b1,
+ MARGIN = 2,
+ FUN = function(b) {
+ curve(b[1] + b[2] * (1 - 2^(-x/b[3])),
+ add = TRUE, col = rgb(1, 1, 0, 0.25), n = 51)
+ invisible()
+ })
+points(Gain ~ A, data = turk0)
+abline(v = 0.2, lty = 2)
+
+#-----------------------------------------------------------------------
+# Estimativa do valor da função em x = 0.2.
+
+p0 <- apply(b0, MARGIN = 2,
+ FUN = function(b){
+ b[1] + b[2] * 0.2/(b[3] + 0.2)
+ })
+p1 <- apply(b1, MARGIN = 2,
+ FUN = function(b){
+ b[1] + b[2] * (1 - 2^(-0.2/b[3]))
+ })
+
+cbind(Est = rbind(mean(p0), mean(p1)),
+ StdErr = rbind(sd(p0), sd(p1)))
+```
+
+## Usando o pacote `boot`
+
+```{r}
+#-----------------------------------------------------------------------
+# Usando o pacote boot.
+
+library(boot)
+library(latticeExtra)
+
+# dput(as.list(coef(n1)))
+start <- list(Int = 622.958054146589,
+ Ass = 178.251908347242,
+ Mei = 0.097321927254495)
+
+# Criar uma função com dois argumentos: o data.frame original e um vetor
+# que vai representar o índice das linhas usado para reamostrar dos
+# dados.
+fitmodel <- function(dataset, index) {
+ n0 <- nls(Gain ~ Int + Ass * (1 - 2^(-A/Mei)),
+ data = dataset[index, ],
+ start = start)
+ c(coef(n0), s2 = deviance(n0)/df.residual(n0))
+}
+
+set.seed(321)
+b0 <- boot(data = turk0,
+ statistic = fitmodel,
+ R = 1999)
+
+class(b0)
+methods(class = class(b0))
+
+summary(b0)
+
+# Como são obtidos.
+theta <- b0$t
+
+# Estimativas como conjunto de dados original.
+b0$t0
+
+# Vício.
+apply(theta, MARGIN = 2, FUN = mean) - b0$t0
+
+# Desvio-padrão boostrap.
+apply(theta, MARGIN = 2, FUN = sd)
+
+# Mediana bootstrap.
+apply(theta, MARGIN = 2, FUN = median)
+
+st <- summary(b0)
+str(st)
+
+# Gráfico de pares.
+pairs(b0)
+
+vcov(b0)
+cov2cor(vcov(b0))
+
+# Extrair os vetores de estativativas bootstrap.
+str(b0)
+colnames(b0$t) <- names(b0$t0)
+best <- stack(as.data.frame(b0$t))
+names(best)
+
+st
+
+densityplot(~values | ind,
+ data = best,
+ scales = "free",
+ as.table = TRUE) +
+ layer(panel.abline(v = st$original[which.packet()] +
+ c(0, st$bootBias[which.packet()]),
+ col = c(1, 2),
+ lty = 2))
+
+# Intervalos de confiança.
+confint(b0, type = "norm")
+confint(b0, type = "basic")
+confint(b0, type = "perc")
+confint(b0, type = "bca")
+```
+
+# Monte Carlo
+
+## Processo pontual
+
+```{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,
+ axes = FALSE,
+ ann = FALSE)
+lines(x = c(0, 1, 1, 0, 0),
+ y = c(0, 0, 1, 1, 0),
+ lty = 2)
+
+# 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: localização uniforme no quadrado.
+M <- replicate(9999, {
+ # As localizações de n pontos.
+ loc <- cbind(x = runif(n),
+ y = runif(n))
+ # A menor distância entre eles sob H_0.
+ min(dist(loc))
+})
+
+# Concatena as estatísticas simuladas com a observada.
+M <- c(m, M)
+
+# Gráfico da distribuição acumulada empírica sob H_0.
+plot(ecdf(M))
+abline(h = c(0.025, 0.975), lty = 2, col = 2)
+abline(v = m, col = 2)
+
+# Gráfico da distribuição acumulada empírica sob H_0.
+plot(density(M))
+abline(v = m, col = 2)
+abline(v = quantile(M, c(0.025, 0.975)), lty = 2, col = 2)
+rug(M)
+
+# 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,
+ axes = FALSE,
+ ann = FALSE)
+lines(x = c(0, 1, 1, 0, 0),
+ y = c(0, 0, 1, 1, 0),
+ lty = 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(499, {
+ 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, col = "gray50")
+}
+lines(ecdf(D[, 1]), cex = 0, col = 2, lwd = 2)
+```
+
+## Amostragem por conjuntos ordenados
+
+```{r}
+#-----------------------------------------------------------------------
+# 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 final.
+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))
+```