Bootstrap
Aula 12 — Módulo 2: Inferência Estatística
Dois distritos dominam a base
Manhattan e Brooklyn somam aproximadamente 85% dos anúncios.
Código do gráfico
A média é puxada pela cauda
Na base original, a mediana é US$ 106 e a média é US$ 153.
Para as análises seguintes, removemos 11 preços iguais a zero e o 1% superior. Permanecem 48.410 anúncios entre US$ 1 e US$ 799; essa limpeza muda explicitamente a população-alvo.
Código do gráfico
fig, axes = plt.subplots(1, 2, figsize=(11, 4.8))
sns.histplot(raw.loc[raw["price"] > 0, "price"], bins=70, color=BLUE, ax=axes[0])
axes[0].set_xlim(0, 1000)
axes[0].set_xlabel("Preço por noite (US$)")
axes[0].set_ylabel("Anúncios")
axes[0].set_title("Escala original")
sns.histplot(np.log10(raw.loc[raw["price"] > 0, "price"]), bins=45, color=ORANGE, ax=axes[1])
axes[1].set_xlabel("log10 do preço")
axes[1].set_ylabel("")
axes[1].set_title("Escala logarítmica")
fig.suptitle("Preços são fortemente assimétricos", fontweight="bold")
finish("distribuicao-precos.png")Localização importa
Mas comparar todos os anúncios mistura localização e tipo de acomodação.
Código do gráfico
order = ["Bronx", "Queens", "Staten Island", "Brooklyn", "Manhattan"]
plt.figure(figsize=(10, 5.4))
sns.boxplot(data=airbnb, x="neighbourhood_group", y="price", order=order,
showfliers=False, color=LIGHT, medianprops={"color": ORANGE, "linewidth": 3})
plt.xlabel("")
plt.ylabel("Preço por noite (US$)")
plt.title("A localização desloca toda a distribuição de preços")
finish("precos-distrito.png")O tipo de acomodação também importa
Manhattan tem uma proporção maior de acomodações inteiras.
Código do gráfico
composition = pd.crosstab(raw["neighbourhood_group"], raw["room_type"], normalize="index")
composition = composition.loc[["Manhattan", "Brooklyn"], ["Entire home/apt", "Private room", "Shared room"]]
composition.plot.bar(stacked=True, figsize=(8.8, 5), color=[BLUE, ORANGE, PURPLE])
plt.xticks(rotation=0)
plt.xlabel("")
plt.ylabel("Proporção")
plt.title("Comparar distritos exige controlar o tipo de acomodação")
plt.legend(title="Tipo", bbox_to_anchor=(1.02, 1), loc="upper left")
finish("composicao-tipo.png")A mediana acompanha o anúncio típico
A cauda extrema afeta pouco o ponto em que a curva acumulada cruza 50%.
Código do gráfico
plt.figure(figsize=(9.5, 5.2))
for values, label, color in [(brooklyn, "Brooklyn", ORANGE), (manhattan, "Manhattan", BLUE)]:
sns.ecdfplot(values, label=label, color=color, lw=3)
plt.axvline(np.median(values), color=color, ls="--", lw=2)
plt.xlim(0, 500)
plt.xlabel("Preço por noite (US$)")
plt.ylabel("Proporção acumulada")
plt.title("A mediana compara o centro sem obedecer à cauda")
plt.legend()
finish("ecdf-amostra.png")A densidade na mediana controla a precisão
As duas populações têm mediana zero. Com pouca massa perto do centro, pequenas mudanças na amostra fazem a mediana saltar entre valores distantes.
O bootstrap evita estimar explicitamente essas densidades, mas ainda exige uma amostra representativa e uma reamostragem compatível com a coleta.
Código do gráfico
def normal_pdf(x, mean, sd):
return np.exp(-0.5 * ((x - mean) / sd) ** 2) / (sd * np.sqrt(2 * np.pi))
x_grid = np.linspace(-4, 4, 800)
pdf_concentrada = normal_pdf(x_grid, 0, 1)
pdf_bimodal = (0.5 * normal_pdf(x_grid, -1.6, 0.8)
+ 0.5 * normal_pdf(x_grid, 1.6, 0.8))
median_rng = np.random.default_rng(1205)
n_rep_median, n_median = 5000, 51
sample_concentrada = median_rng.normal(0, 1, (n_rep_median, n_median))
componentes = median_rng.integers(0, 2, (n_rep_median, n_median))
medias_bimodais = np.where(componentes == 0, -1.6, 1.6)
sample_bimodal = median_rng.normal(medias_bimodais, 0.8)
medianas_concentrada = np.median(sample_concentrada, axis=1)
medianas_bimodal = np.median(sample_bimodal, axis=1)
fig, axes = plt.subplots(1, 2, figsize=(11.5, 4.8), sharex=True)
axes[0].plot(x_grid, pdf_concentrada, color=BLUE, lw=3,
label=r"alta densidade em $m$")
axes[0].plot(x_grid, pdf_bimodal, color=ORANGE, lw=3,
label=r"baixa densidade em $m$")
axes[0].axvline(0, color=INK, ls="--", lw=2, label="mediana comum: 0")
axes[0].set_title("Distribuições populacionais")
axes[0].set_xlabel("Valor")
axes[0].set_ylabel("Densidade")
def kde_on_grid(values, grid):
bandwidth = 1.06 * values.std(ddof=1) * len(values) ** (-1 / 5)
z = (grid[:, None] - values[None, :]) / bandwidth
return np.exp(-0.5 * z**2).mean(axis=1) / (bandwidth * np.sqrt(2*np.pi))
kde_concentrada = kde_on_grid(medianas_concentrada, x_grid)
kde_bimodal = kde_on_grid(medianas_bimodal, x_grid)
axes[1].fill_between(x_grid, kde_bimodal, color=ORANGE, alpha=0.30)
axes[1].plot(x_grid, kde_bimodal, color=ORANGE, lw=3,
label="baixa densidade em m")
axes[1].fill_between(x_grid, kde_concentrada, color=BLUE, alpha=0.38)
axes[1].plot(x_grid, kde_concentrada, color=BLUE, lw=3,
label="alta densidade em m")
axes[1].axvline(0, color=INK, ls="--", lw=2)
axes[1].set_title("Medianas de 5.000 amostras, n = 51")
axes[1].set_xlabel("Mediana amostral")
axes[1].set_ylabel("Densidade")
for ax in axes:
ax.set_xlim(-4, 4)
ax.legend()
fig.suptitle("Pouca densidade perto da mediana aumenta sua variabilidade",
fontweight="bold")
finish("densidade-mediana-precisao.png")A origem dos 63,2% no bootstrap
Repetições e ausências surgem porque o sorteio ocorre com reposição.
Para uma observação específica, a probabilidade de não aparecer em nenhum dos \(n\) sorteios é
\[ P(\text{ausente})=\left(1-\frac{1}{n}\right)^n. \]
Se \(D_n\) é o número de observações distintas, escrevemos \(D_n=\sum_{i=1}^n I_i\), onde \(I_i=1\) quando a observação \(i\) aparece. Logo,
\[ \begin{aligned} \mathbb{E}(D_n) &=\sum_{i=1}^n\mathbb{E}(I_i) \\ &=n\left[1-\left(1-\frac{1}{n}\right)^n\right]. \end{aligned} \]
Portanto,
\[ \frac{\mathbb{E}(D_n)}{n}\longrightarrow 1-e^{-1}\approx 0{,}632. \]
A incerteza ganha forma
A distribuição não precisa ser perfeitamente Normal para ser útil.
Código do gráfico
plt.figure(figsize=(9.3, 5.2))
sns.histplot(boot, bins=np.arange(18, 65, 1), color=BLUE)
plt.axvline(theta_hat, color=ORANGE, lw=3, label=f"diferença observada = US$ {theta_hat:.0f}")
plt.xlabel("Mediana Manhattan − mediana Brooklyn (US$)")
plt.ylabel("Réplicas")
plt.title("As réplicas revelam a incerteza da estimativa")
plt.legend()
finish("distribuicao-bootstrap.png")O intervalo percentil é direto
Um intervalo bootstrap percentil de nível \(1-\alpha\) é:
\[ \left[ q_{\alpha/2}(\widehat\theta^*), q_{1-\alpha/2}(\widehat\theta^*) \right]. \]
Ele usa os quantis da própria distribuição bootstrap.
\(q_\gamma(\widehat\theta^*)\) é o quantil de ordem \(\gamma\) das réplicas e \(\alpha\) é a probabilidade total deixada fora do intervalo; para 95%, \(\alpha=0{,}05\).
Os 95% centrais vão de US$ 30 a US$ 51 e permanecem acima de zero.
Código do gráfico
plt.figure(figsize=(9.3, 5.2))
sns.kdeplot(boot, fill=True, color=LIGHT, linewidth=0)
inside = (boot >= ci[0]) & (boot <= ci[1])
sns.kdeplot(boot[inside], fill=True, color=BLUE, linewidth=0)
plt.axvline(ci[0], color=ORANGE, ls="--", lw=3)
plt.axvline(ci[1], color=ORANGE, ls="--", lw=3)
plt.text(ci[0], .012, f"2,5%\nUS$ {ci[0]:.0f}", ha="right", color=ORANGE, fontweight="bold")
plt.text(ci[1], .012, f"97,5%\nUS$ {ci[1]:.0f}", ha="left", color=ORANGE, fontweight="bold")
plt.xlabel("Diferença de medianas (US$)")
plt.ylabel("Densidade")
plt.title("O intervalo percentil preserva os 95% centrais")
finish("intervalo-percentil.png")Quantas réplicas são suficientes?
Para intervalos, alguns milhares de réplicas costumam ser um ponto de partida prático.
Código do gráfico
checkpoints = np.arange(100, 3001, 100)
lower = np.array([np.quantile(boot[:b], .025) for b in checkpoints])
upper = np.array([np.quantile(boot[:b], .975) for b in checkpoints])
plt.figure(figsize=(9.4, 5.2))
plt.plot(checkpoints, lower, color=ORANGE, lw=2, label="limite inferior")
plt.plot(checkpoints, upper, color=BLUE, lw=2, label="limite superior")
plt.xlabel("Número de réplicas R")
plt.ylabel("Limite do IC (US$)")
plt.title("Mais réplicas reduzem o ruído de Monte Carlo")
plt.legend()
finish("estabilidade-b.png")Para médias, bootstrap e TCL ficam próximos
Com a amostra fixa, aumentar \(R\) reduz a oscilação de Monte Carlo. A proximidade com o TCL ocorre porque a diferença de médias tem distribuição aproximadamente Normal — não porque \(R\) obrigue os métodos a coincidir.
Código do gráfico
mean_hat = manhattan.mean() - brooklyn.mean()
se_tcl = np.sqrt(
manhattan.var(ddof=1) / len(manhattan)
+ brooklyn.var(ddof=1) / len(brooklyn)
)
ci_tcl = np.array([mean_hat - 1.96*se_tcl,
mean_hat + 1.96*se_tcl])
n_rep_max = 5000
mean_boot = (
rng.choice(manhattan, (n_rep_max, len(manhattan)), replace=True).mean(axis=1)
- rng.choice(brooklyn, (n_rep_max, len(brooklyn)), replace=True).mean(axis=1)
)
checkpoints = np.array([25, 50, 100, 200, 400, 800,
1200, 2000, 3000, 5000])
ci_boot_by_rep = np.array([
np.quantile(mean_boot[:n_rep], [.025, .975]) for n_rep in checkpoints
])
plt.plot(checkpoints, ci_boot_by_rep[:, 0], "o-", label="Bootstrap: inferior")
plt.plot(checkpoints, ci_boot_by_rep[:, 1], "o-", label="Bootstrap: superior")
plt.axhline(ci_tcl[0], ls="--", label="TCL: inferior")
plt.axhline(ci_tcl[1], ls="--", label="TCL: superior")
plt.xscale("log"); plt.legend()
finish("estabilidade-media-bootstrap-tcl.png")Mais dados ajudam mais que mais réplicas
Depois de estabilizar \(R\), ampliar a amostra é o que realmente estreita o intervalo.
Código do gráfico
sizes = np.array([50, 100, 200, 400, 600])
widths = []
for n in sizes:
widths_n = []
for repeat in range(6):
m = rng.choice(manhattan, n, replace=False)
b = rng.choice(brooklyn, n, replace=False)
values = bootstrap_difference(m, b, n_rep=250, seed=1000 + repeat + n)
q = np.quantile(values, [.025, .975])
widths_n.append(q[1] - q[0])
widths.append(widths_n)
plt.figure(figsize=(9.2, 5.2))
sns.boxplot(data=pd.DataFrame(widths, index=sizes).T, color=LIGHT,
medianprops={"color": ORANGE, "linewidth": 3})
plt.xlabel("Observações por distrito")
plt.ylabel("Largura do IC 95% (US$)")
plt.title("Mais dados estreitam o intervalo")
finish("tamanho-largura.png")Métodos podem discordar
Discordância relevante é um sinal para investigar assimetria, viés e tamanho amostral. Nesta amostra, os empates tornam o jackknife da mediana degenerado; usamos \(a=0\), e o BCa conserva apenas a correção de viés.
O cálculo detalhado do BCa e um exercício resolvido permanecem no material de estudo como leitura complementar.
Código do gráfico
from statistics import NormalDist
se_boot = boot.std(ddof=1)
# BCa: z0 vem das réplicas e a vem do jackknife
normal = NormalDist()
prop_below = np.clip(
np.mean(boot < theta_hat),
1 / (2 * len(boot)),
1 - 1 / (2 * len(boot)),
)
z0 = normal.inv_cdf(prop_below)
jack_m = np.array([
np.median(np.delete(manhattan, i)) - np.median(brooklyn)
for i in range(len(manhattan))
])
jack_b = np.array([
np.median(manhattan) - np.median(np.delete(brooklyn, i))
for i in range(len(brooklyn))
])
jack = np.concatenate([jack_m, jack_b])
jack_center = jack.mean()
jack_diff = jack_center - jack
jack_scale = np.sum(jack_diff**2)
a = 0.0 if jack_scale == 0 else (
np.sum(jack_diff**3) / (6 * jack_scale**1.5)
)
z_alpha = [normal.inv_cdf(.025), normal.inv_cdf(.975)]
alpha_bca = [
normal.cdf(z0 + (z0 + z) / (1 - a * (z0 + z)))
for z in z_alpha
]
ci_bca = np.quantile(boot, alpha_bca)
intervals = {
"Percentil": ci,
"Básico": np.array([2 * theta_hat - ci[1], 2 * theta_hat - ci[0]]),
"Normal": np.array([theta_hat - 1.96 * se_boot, theta_hat + 1.96 * se_boot]),
"BCa (a = 0)": ci_bca,
}
plt.figure(figsize=(9, 4.7))
for y, (name, limits) in enumerate(intervals.items()):
plt.plot(limits, [y, y], color=BLUE, lw=5)
plt.scatter(theta_hat, y, color=ORANGE, s=90, zorder=3)
plt.yticks(range(len(intervals)), intervals.keys())
plt.axvline(0, color=INK, ls="--")
plt.xlabel("Diferença de medianas (US$)")
plt.title("Métodos diferentes podem produzir limites diferentes")
finish("metodos-intervalo.png")O p-valor vem da distribuição bootstrap sob H0
Com \(R=10.000\) réplicas, nenhuma diferença nula atingiu \(|T_{obs}|=40\). Usando a correção de Monte Carlo,
\[ \widehat p_{boot} =\frac{1+\sum_{b=1}^{R}\mathbb 1\{|T_b^{*(0)}|\geq|T_{obs}|\}}{R+1} =\frac{1}{10.001}\approx0{,}0001. \]
\(T_{obs}\) é a diferença observada; \(T_b^{*(0)}\) é a diferença na réplica \(b\) gerada sob \(H_0\); e \(R\) é o total de réplicas. O resultado fornece forte evidência contra a igualdade das medianas, sob as hipóteses de amostragem e independência adotadas.
Código do gráfico e do p-valor
manhattan_h0 = manhattan - np.median(manhattan)
brooklyn_h0 = brooklyn - np.median(brooklyn)
test_rng = np.random.default_rng(1212)
n_rep_test = 10_000
idx_m = test_rng.integers(0, len(manhattan_h0), (n_rep_test, len(manhattan_h0)))
idx_b = test_rng.integers(0, len(brooklyn_h0), (n_rep_test, len(brooklyn_h0)))
boot_h0 = (np.median(manhattan_h0[idx_m], axis=1)
- np.median(brooklyn_h0[idx_b], axis=1))
extreme_h0 = np.abs(boot_h0) >= np.abs(theta_hat)
p_boot = (1 + extreme_h0.sum()) / (len(boot_h0) + 1)Duas distribuições, duas perguntas
O bootstrap usual fornece o intervalo do efeito. O bootstrap sob \(H_0\) fornece \(p\approx0{,}0044\), muito próximo de \(0{,}0041\) obtido pelo teste \(z\) na Aula 11. Aqui, ambos os resultados vêm da reamostragem direta das linhas 0/1; a pequena diferença vem do erro de Monte Carlo.
Validade significa cobertura em repetições
Um intervalo é válido quando sua cobertura real se aproxima do nível nominal:
\[ P_F\{\theta\in CI(X_1,\ldots,X_n)\}\approx 1-\alpha. \]
\(P_F\) representa amostras repetidas de \(F\) e \(1-\alpha\) é a cobertura nominal. Na prática, a população verdadeira não está disponível para essa verificação.
Código do gráfico
population_median = np.median(focus.loc[focus["neighbourhood_group"].eq("Manhattan"), "price"])
population_median -= np.median(focus.loc[focus["neighbourhood_group"].eq("Brooklyn"), "price"])
cover = []
for i in range(20):
m = rng.choice(focus.loc[focus["neighbourhood_group"].eq("Manhattan"), "price"], 250, replace=False)
b = rng.choice(focus.loc[focus["neighbourhood_group"].eq("Brooklyn"), "price"], 250, replace=False)
values = bootstrap_difference(m, b, n_rep=300, seed=3000 + i)
lo, hi = np.quantile(values, [.025, .975])
cover.append((lo, hi, lo <= population_median <= hi))
plt.figure(figsize=(9, 6.3))
for i, (lo, hi, covered) in enumerate(cover):
plt.plot([lo, hi], [i, i], color=BLUE if covered else ORANGE, lw=2)
plt.axvline(population_median, color=INK, lw=3)
plt.xlabel("Diferença de medianas (US$)")
plt.ylabel("Amostras repetidas")
plt.title("Validade significa cobertura no longo prazo")
finish("cobertura-bootstrap.png")Dependência e unidade de reamostragem
Anúncios do mesmo anfitrião podem compartilhar preço, localização e decisões de gestão. Reamostrar anúncios como independentes pode subestimar a incerteza. Quando a dependência ocorre por anfitrião, reamostramos anfitriões inteiros.
Código do gráfico
host_counts = sample.groupby(["neighbourhood_group", "host_id"]).size()
fig, axes = plt.subplots(1, 2, figsize=(10.5, 4.8))
sns.histplot(host_counts, discrete=True, color=BLUE, ax=axes[0])
axes[0].set_xlim(.5, 6.5)
axes[0].set_xlabel("Anúncios do mesmo anfitrião na amostra")
axes[0].set_ylabel("Anfitriões")
axes[0].set_title("Há agrupamento por anfitrião")
naive = bootstrap_difference(manhattan, brooklyn, n_rep=1000, seed=9)
cluster_estimates = []
for _ in range(100):
pieces = []
for borough in ["Manhattan", "Brooklyn"]:
group = sample[sample["neighbourhood_group"].eq(borough)]
hosts = group["host_id"].unique()
chosen = rng.choice(hosts, len(hosts), replace=True)
values = np.concatenate([group.loc[group["host_id"].eq(h), "price"].to_numpy() for h in chosen])
pieces.append(np.median(values))
cluster_estimates.append(pieces[0] - pieces[1])
sns.kdeplot(naive, color=BLUE, lw=3, label="por anúncio", ax=axes[1])
sns.kdeplot(cluster_estimates, color=ORANGE, lw=3, label="por anfitrião", ax=axes[1])
axes[1].set_xlabel("Diferença de medianas (US$)")
axes[1].set_ylabel("Densidade")
axes[1].set_title("A unidade de reamostragem importa")
axes[1].legend()
finish("dependencia-host.png")
print(f"linhas={len(raw)}; análise={len(airbnb)}; foco={len(focus)}")
print(f"medianas: Manhattan={np.median(manhattan):.0f}; Brooklyn={np.median(brooklyn):.0f}")
print(f"diferença={theta_hat:.0f}; IC95%=[{ci[0]:.0f}, {ci[1]:.0f}]")