From b76ca99bb81fa2fd6b97fbb8fec19f9bb352b46b Mon Sep 17 00:00:00 2001
From: Walmes Zeviani <walmes@ufpr.br>
Date: Wed, 4 Jul 2018 17:21:24 -0300
Subject: [PATCH] =?UTF-8?q?Conclui=20arquivo=20com=20exemplos=20de=20m?=
=?UTF-8?q?=C3=A9todos=20intensivos.?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit
---
scripts/metodos-intensivos.Rmd | 237 +++++++++++++++++++++++----------
1 file changed, 168 insertions(+), 69 deletions(-)
diff --git a/scripts/metodos-intensivos.Rmd b/scripts/metodos-intensivos.Rmd
index 3b29cd0..896c941 100644
--- a/scripts/metodos-intensivos.Rmd
+++ b/scripts/metodos-intensivos.Rmd
@@ -19,18 +19,44 @@ opts_chunk$set(cache = TRUE,
message = FALSE)
```
+Esse documento contém exemplos de aplicação de métodos
+computacionalmente intensivos em problemas de inferência estatística.
+Tais problemas são tipicamente testes de hipóteses, estimação intervalar
+e estudo de propriedades de estimadores. Os métodos exemplificados são:
+
+ * Testes de aleatorização
+ * Métodos Jackknife.
+ * Métodos bootstrap.
+ * Métodos Monte Carlo.
+
+As sessões a seguir estão dividas em exemplos.
+
# Testes de aleatorização
## Uma senhora toma chá
+Visite
+
+ * <https://en.wikipedia.org/wiki/Lady_tasting_tea>.
+ * Série de vídeos.
+ * <https://www.youtube.com/watch?v=vYVr50hjFbQ>.
+ * <https://www.youtube.com/watch?v=xh20btybjp4>.
+ * <https://www.youtube.com/watch?v=JMDycaD-YwY>.
+ * <https://www.youtube.com/watch?v=dQT6T3Zo-oA>.
+ * <https://www.youtube.com/watch?v=9sgQMHtLPNU>.
+ * <https://www.youtube.com/watch?v=lgs7d5saFFc>.
+
```{r}
#-----------------------------------------------------------------------
-# Uma senhora toma chá.
+# O curioso caso da senhora que toma chá.
# Quantidade de maneiras de gerar dois grupos de 4 xícaras usando 8.
n_poss <- choose(8, 4)
n_poss
+# E se fossem 12 xícaras, 6 de cada?
+choose(12, 6)
+
# De 4 xícaras selecionadas para um grupo, X representa o número de
# acertos.
n_acertar <- sapply(4:0,
@@ -40,15 +66,18 @@ n_acertar <- sapply(4:0,
})
n_acertar
-# Pr(Acerto total) = 1/70 < 5%.
-cumsum(n_acertar)/n_poss
+# Pr(Acerto total) = Pr(X = 4) = 1/70 < 5%.
+cbind(X = 4:0, "Pr(X = x)" = n_acertar/n_poss)
# Comprovando por simulação.
+# Uma particular configuração das xícaras.
v <- rep(0:1, each = 4)
+
+# Com qual frequência, ao acaso, a mesma configuração ocorre?
mean(replicate(n = 100000,
expr = all(sample(v) == v)))
-# Essa forma de fazer a conta é mais realista mas dá na mesma.
+# Essa forma de fazer é mais realista mas o resultado é o memso.
mean(replicate(n = 100000,
expr = all(sample(v) == sample(v))))
```
@@ -59,10 +88,11 @@ mean(replicate(n = 100000,
#-----------------------------------------------------------------------
# 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)
+# Comprimentos de mandíbulas de cães pré-históricos.
+m <- c(120, 107, 110, 116, 114, 111, 113, 117, 114, 112) # machos.
+f <- c(110, 111, 107, 108, 110, 105, 107, 106, 111, 111) # fêmeas.
+# Distribuição acumulada empírica.
plot(ecdf(m), xlim = range(c(m, f)), col = "cyan")
lines(ecdf(f), col = "magenta")
rug(m, col = "cyan")
@@ -72,11 +102,11 @@ rug(f, col = "magenta")
d <- mean(m) - mean(f)
d
-# Todos as combinações possíveis.
+# Todos as combinações possíveis entre 20 objetos tomados 10 a 10.
choose(n = 20, k = 10)
-#--------------------------------------------
-# Com todas as combinações possíveis (exaustão).
+#-----------------------------------------------------------------------
+# Solução com todas as permutações possíveis (exaustão).
# Para construir todas as combinações possíveis.
k <- combn(x = 1:20, m = 10)
@@ -85,33 +115,34 @@ dim(k)
# Vetor com os valores dos dois grupos concatenados.
mf <- c(m, f)
-# Vetor cheio de zeros.
+# Vetor de zeros.
g <- integer(20)
-# Calcula a diferença para todo arranjo possível.
+# Calcula a diferença para cada uma das permutações.
D <- apply(k,
MARGIN = 2,
FUN = function(i) {
- # Alguns elementos do vetor convertidos para 1.
+ # Alguns zeros são convertidos para 1 criando 2 grupos.
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.
+# Histograma da distribuição da estatística sob H_0.
hist(D, col = "gray50")
rug(D)
abline(v = d, col = 2)
+# Distribuição acumulada empírica da estatística sob H_0.
plot(ecdf(D), cex = 0)
rug(D)
abline(v = d, col = 2)
-# P-valor do teste.
+# P-valor do teste (bilateral).
2 * sum(D >= d)/length(D)
-#--------------------------------------------
-# Com simulação (não necessáriamente exaustivo).
+#-----------------------------------------------------------------------
+# Com reamostragem sem reposição (não necessáriamente exaustivo).
# Variáveis que indentifica os grupos.
g <- rep(1:2, each = 10)
@@ -136,7 +167,7 @@ plot(ecdf(D), cex = 0)
rug(D)
abline(v = d, col = 2)
-# P-valor do teste.
+# P-valor do teste (bilateral).
2 * sum(D >= d)/length(D)
```
@@ -149,12 +180,12 @@ abline(v = d, col = 2)
# 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.
+# Estatística de teste na amostra original: coeficiente de correlação de
+# Pearson.
r0 <- cor(x, y, method = "pearson")
r0
@@ -166,6 +197,11 @@ head(X)
# As 120 correlações obtidas para cada arranjo.
r <- apply(X, MARGIN = 1, FUN = cor, y = y, method = "spearman")
+# Gráfico da distribuição da estatística de teste.
+plot(ecdf(r), cex = 0)
+rug(r)
+abline(v = r0, col = 2)
+
# P-valor do teste.
2 * sum(r >= r0)/length(r)
```
@@ -181,25 +217,38 @@ r <- apply(X, MARGIN = 1, FUN = cor, y = y, method = "spearman")
#-----------------------------------------------------------------------
# Índice de Moran para medir dependência espacial.
-# Coordenadas dos eventos em uma malha regular.
+# Coordenadas dos eventos em uma malha regular 8 x 8.
x <- 1:8
y <- 1:8
-# Para criar a matriz de pesos.
-ind <- expand.grid(i = 1:length(x), j = 1:length(y))
+# Construção da matriz de pesos que determina a vizinhança entre
+# observações.
+ind <- expand.grid(i = 1:length(x),
+ j = 1:length(y))
+# Função que determina o peso entre duas localizações na malha.
f <- function(i, j) {
u <- min(3, sum(abs(ind[i, ] - ind[j, ])))
- c(0, 1, sqrt(1/2), 0)[u + 1]
+ w <- c(0, 1, sqrt(1/2), 0)[u + 1]
+ return(w)
}
-w <- matrix(0, nrow(ind), nrow(ind))
+# Cria os pesos, matriz (8^2) x (8^2) = 64 x 64.
+w <- matrix(0, nrow = nrow(ind), ncol = nrow(ind))
for (i in 1:nrow(ind)) {
for (j in 1:nrow(ind)) {
w[i, j] <- f(i, j)
}
}
+# Normaliza.
w <- w/sum(w)
-image(w)
+# Gráfico. Valores claros indicam maior peso entre observações.
+image(w, asp = 1, col = gray.colors(100))
+
+# Lógica do índica de Moran: correlação entre valores observados e
+# média dos vizinhos. Exemplo com valores simulados.
+xx <- rnorm(64)
+cor(cbind("Valores observados" = xx,
+ "Média dos vizinhos" = as.vector(xx %*% w)))
# Índice de Moran.
moran <- function(x, weights) {
@@ -211,7 +260,7 @@ moran <- function(x, weights) {
as.vector(z %*% weights %*% (z * sqrt(n/(n - 1))))
}
-# Teste de permutação com gráfico.
+# Teste de permutação com saída gráfica.
ppt <- function(z, w, N = 10000, stat, ...) {
# Índice de Moran por reamostragem.
sim <- replicate(N,
@@ -238,6 +287,9 @@ z <- matrix(rexp(length(x) * length(y),
length(x))
image(log(z), main = "Com dependência")
+cor(cbind("Valores observados" = as.vector(z),
+ "Média dos vizinhos" = as.vector(as.vector(z) %*% w)))
+
# Índice de Moran com dados originais.
stat <- moran(z, w)
stat
@@ -250,6 +302,9 @@ ppt(z, w, stat = stat, main = "I de Moran", xlab = "I")
z <- matrix(rnorm(length(x) * length(y), 0, 1/2), length(x))
image(z, main = "Sem dependência")
+cor(cbind("Valores observados" = as.vector(z),
+ "Média dos vizinhos" = as.vector(as.vector(z) %*% w)))
+
# Índice de Moran com dados originais.
stat <- moran(z, w)
stat
@@ -266,9 +321,12 @@ ppt(z, w, stat = stat, main = "I de Moran", xlab = "I")
#-----------------------------------------------------------------------
# Aplicação de Jackknife em modelos de regressão não linear.
-library(latticeExtra)
library(car)
library(alr3)
+library(latticeExtra)
+
+# Documentação do conjunto de dados.
+# help(segreg)
# Curva "palpite".
start <- list(th0 = 75, th1 = 0.5, th2 = 50)
@@ -277,7 +335,7 @@ xyplot(C ~ Temp, data = segreg) +
0 * (x < th2), lty = 2),
data = start)
-# Ajuste.
+# Ajuste do modelo não linear aos dados (modelo platô-linear).
n0 <- nls(C ~ th0 + th1 * (Temp - th2) * (Temp >= th2) +
0 * (Temp < th2),
data = segreg,
@@ -286,7 +344,7 @@ n0 <- nls(C ~ th0 + th1 * (Temp - th2) * (Temp >= th2) +
# Estimativas e medidas de ajuste.
summary(n0)
-# Observados e preditos.
+# Valores observados curva ajustada.
xyplot(C ~ Temp, data = segreg) +
layer(panel.curve(th0 + th1 * (x - th2) * (x >= th2) +
0 * (x < th2)),
@@ -295,11 +353,10 @@ xyplot(C ~ Temp, data = segreg) +
#-----------------------------------------------------------------------
# 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.
+# Matriz vazia para armazenar os pseudo-valores durante iterações.
thetaj <- matrix(0, nrow = n, ncol = length(theta))
colnames(thetaj) <- names(theta)
@@ -314,16 +371,19 @@ for (i in 1:n) {
thetaj[i, ] <- n * theta - (n - 1) * coef(n1)
}
-# Os primeiros pseudo valores.
+# Alguns pseudo-valores.
head(thetaj)
+tail(thetaj)
-# Estimativas pontuais.
+# Estimativas pontuais de Jackknife.
jk <- colMeans(thetaj)
-cbind(MLE = theta, Jack = jk, Vício = theta - jk)
+rbind(MLE = theta,
+ Jackknife = jk,
+ Vício = theta - jk)
# Erros padrões.
-cbind(MLE = summary(n0)$coefficients[, "Std. Error"],
- Jack = apply(thetaj, MARGIN = 2, sd)/sqrt(n))
+rbind(MLE = summary(n0)$coefficients[, "Std. Error"],
+ Jackknife = apply(thetaj, MARGIN = 2, sd)/sqrt(n))
# Pacote com função pronta para Jackknife para modelos não lineares.
library(nlstools)
@@ -340,6 +400,8 @@ summary(j0)
#-----------------------------------------------------------------------
# Inferência para a DL 50.
+# Dados no pacote labestData: https://github.com/pet-estatistica/labestData.
+
# library(labestData)
# data(PaulaTb3.12)
# str(PaulaTb3.12)
@@ -349,40 +411,50 @@ summary(j0)
# 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")
+ 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)
+# Visualização dos dados.
+xyplot(resp ~ vol,
+ data = PaulaTb3.12,
+ groups = resp)
+
+# Ajuste do modelo considerando resposta binária.
m0 <- glm(resp ~ vol,
data = PaulaTb3.12,
family = binomial)
summary(m0)
-# DL_50.
+# Estimativa da DL_50 a partir do modelo ajustado.
dl <- -coef(m0)[1]/coef(m0)[2]
-# str(m0$family)$link
-
+# Curva ajustada.
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úmero de observações e vetor de índices.
n <- nrow(PaulaTb3.12)
i <- 1:n
-# Estimativas por Jackknife.
+# Estimativas por Jackknife: pseudo-valores.
DL <- sapply(i,
FUN = function(i) {
m0 <- glm(resp ~ vol,
@@ -394,24 +466,30 @@ DL <- sapply(i,
n * dl - (n - 1) * DL
})
+# Estimativa pontual e erro-padrão por Jackknife.
m <- mean(DL)
e <- sd(DL)/sqrt(n)
+# Comparando com as estimativas de Wald (usa método Delta).
cbind(Est = rbind(MLE = dl,
- Jack = m),
+ Jackknife = m),
StdErr = rbind(MLE = car::deltaMethod(m0, g = "-(Intercept)/vol")["SE"],
Jack = e))
+# Distribuição acumulada empírica.
plot(ecdf(DL))
rug(DL)
+abline(v = c(m, dl), lty = c(1, 2), col = 2)
+# Densidade empírica.
plot(density(DL))
rug(DL)
+abline(v = c(m, dl), lty = c(1, 2), col = 2)
# NOTE: alguma pista sobre porque a distribuição fica bimodal?
```
-## Uma situação problemática
+## Condição problemática para uso do método Jackknife
```{r}
#-----------------------------------------------------------------------
@@ -419,6 +497,7 @@ rug(DL)
# Ordena o vetor de valores.
x <- sort(precip)
+x
# Calcula a mediana com os dados originais.
M <- median(x)
@@ -453,6 +532,9 @@ length(x)
library(alr3)
str(turk0)
+# Documentação do conjunto de dados.
+# help(turk0)
+
# Gráfico dos valores observados.
plot(Gain ~ A, data = turk0)
@@ -481,14 +563,16 @@ cbind(coef(n0), coef(n1))
#-----------------------------------------------------------------------
# Criando funções que fazem reamostragem dos dados e ajustam o modelo.
+# Número de observações e índices.
n <- nrow(turk0)
i <- 1:n
-myboot0 <- function(formula, start) {
+# Função que obtém estimativas bootstrap.
+myboot <- function(formula, start) {
n0 <- nls(formula = formula,
data = turk0[sample(i, n, replace = TRUE), ],
start = start)
- coef(n0)
+ return(coef(n0))
}
#-----------------------------------------------------------------------
@@ -496,18 +580,28 @@ myboot0 <- function(formula, start) {
# B = 999 reamostragens.
b0 <- replicate(999,
- myboot0(Gain ~ Int + Ass * A/(Mei + A), start))
+ myboot(Gain ~ Int + Ass * A/(Mei + A), start))
b1 <- replicate(999,
- myboot0(Gain ~ Int + Ass * (1 - 2^(-A/Mei)), start))
+ myboot(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)
+# Estimativas originais.
+mle_est <- rbind(n0 = coef(n0),
+ n0 = coef(n1))
+
+# Estimativas bootstrap.
+bot_est <- rbind(n0 = rowMeans(b0),
+ n1 = rowMeans(b1))
+
+mle_est
+bot_est
+
# Calcula o vício absoluto (bootstrap - original).
-rbind(n0 = rowMeans(b0) - coef(n0),
- n1 = rowMeans(b1) - coef(n1))
+bot_est - mle_est
# Variância do estimador.
rbind(n0 = apply(b0, MARGIN = 1, var),
@@ -531,8 +625,12 @@ rbind(n0 = apply(b0, 1, FUN = eqm),
bts <- cbind(rep(1:2, each = ncol(b0)),
as.data.frame(rbind(t(b0), t(b1))))
+
+# Matriz de diagramas de dispersão das estimativas nas amostras
+# bootstrap.
splom(bts[, -1], groups = bts[, 1])
+# Cada uma das curvas que foram ajustas às amostras bootstrap.
par(mfrow = c(1, 2))
plot(Gain ~ A, data = turk0, type = "n")
apply(b0,
@@ -567,8 +665,12 @@ p1 <- apply(b1, MARGIN = 2,
b[1] + b[2] * (1 - 2^(-0.2/b[3]))
})
-cbind(Est = rbind(mean(p0), mean(p1)),
- StdErr = rbind(sd(p0), sd(p1)))
+# Estimativa e erro-padrão para o valor predito em x = 0.4.
+est <- cbind(rbind(mean(p0), mean(p1)),
+ rbind(sd(p0), sd(p1)))
+dimnames(est) <- list(c("n0", "n1"),
+ c("Est", "StdErr"))
+est
```
## Usando o pacote `boot`
@@ -578,7 +680,6 @@ cbind(Est = rbind(mean(p0), mean(p1)),
# Usando o pacote boot.
library(boot)
-library(latticeExtra)
# dput(as.list(coef(n1)))
start <- list(Int = 622.958054146589,
@@ -661,8 +762,9 @@ confint(b0, type = "bca")
#-----------------------------------------------------------------------
# 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))
+# layout(1)
+# plot(x = NULL, y = 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))
@@ -835,9 +937,6 @@ 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)) +
--
GitLab