4 - Funções Wiggly¶

No capítulo anterior vimos como construir o nosso primeiro modelo de regressão linear, no qual estimamos a altura de uma pessoa usando a informação de seu próprio peso como um auxílio informativo. Temos, assim, como resultado do nosso modelo, a distribuição à posteriori dos parâmetros que foram estimados, \(\alpha\), \(\beta\) e o \(\sigma\).

Modelo Linear Simples

\[ altura_i = \alpha + \beta (peso_i - peso\_médio)\]

Das distribuições à posteriori dos parâmetros podemos pegar apenas algum valor "preciso" desses parâmetros, ou seja, geralmente nós tendemos pensar apenas com os valores esperados (valores médios) da distribuição do \(\alpha\) e do \(\beta\) e, com isso, traçar uma linha reta com esses valores.

Olhando para os resultados das inferências dos parâmetros, quando o peso está em seu valor médio, temos que \(\alpha\) significa o valor médio da altura!

Percebeu? Percebeu como isso é muito bonito! Agora temos uma interpretação decente para o \(\alpha\)!

Agora, o que significa o nosso \(\beta\)? Temos que para cada unidade de \(peso\) que aumentarmos a altura \(h\) também aumentará em \(\beta\) unidades. Assim, para nosso exemplo, cada kilo que aumentarmos, a altura também tende a aumentar, em média, \(0.90\)cm.

Mas sabemos que isso é insuficiente, porque queremos obter a incerteza a partir desse gráfico. A estatística bayesiana não lhe dará uma única estimativa pontual, mas sim, dará a SUA incerteza, que é comunicada a nós pela distribuição à posteriori. Essa incerteza é representada pelo número infinito de linhas, e cada uma dessas linhas são classificadas pela sua plausibilidade relativa (sua probabilidade) em comparação todas as outras linhas!

Warning

OBS: A partir de agora irei definir as variáveis em Inglês, mas os comentários continuarão em português para facilitar o entendimento seguir com o propósito do material.

Em ambientes profissionais é recomendado escrever tudo em inglês.

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import stan
import nest_asyncio

import logging
logging.getLogger("pystan").propagate=False

# Desbloqueio do asyncIO do jupyter
nest_asyncio.apply()

# Definindo o plano de fundo cinza claro para todos os gráficos feitos no matplotlib
plt.rcParams['axes.facecolor'] = 'lightgray'

# Lendo os dados

# Os dados podem serem obtidos em https://github.com/rmcelreath/rethinking/tree/master/data/Howell1.csv
df = pd.read_csv('./data/Howell1.csv', sep=';')  

weight = df.loc[df.age >= 30, 'weight'].values  # Filtro para pessoas com 30 anos ou mais 
height = df.loc[df.age >= 30, 'height'].values  # Filtro para pessoas com 30 anos ou mais

model_stan = """
    data {
        int N;
        vector[N] height;
        vector[N] weight;
    }
    
    parameters {
        real alpha;
        real beta;
        real<lower=0, upper=50> sigma;  // Priori implícita Uniforme(0, 50)
    }
    
    model {
        alpha ~ normal(178, 20);  // Priori para Alpha
        beta ~ lognormal(0, 1);  // Priori para Beta
        
        height ~ normal(alpha + beta * weight, sigma);
    }
"""

# =============================================================
#    Estimativa da altura explicada com a variável peso.
# =============================================================
# Reescrevendo o modelo anterior 

# Lembrando que estamos usando o (x_i - x_barra) e não apenas x_i
weight_adjust = weight - weight.mean() 

my_data = {
    'N': len(weight),
    'height': height,
    'weight': weight_adjust,
}

posteriori = stan.build(model_stan, data=my_data)
fit = posteriori.sample(num_chains=4, num_samples=1000)

# Recuperando os parâmetros
alpha = fit['alpha'].flatten()
beta = fit['beta'].flatten()
sigma = fit['sigma'].flatten()

# ==========================================
#  Plotando os histogramas dos Parâmentros
# ==========================================

fig, [ax1, ax2, ax3] = plt.subplots(1, 3, figsize=(17, 9))

# Parâmetro: alpha
ax1.hist(alpha, density=True, rwidth=0.9, bins=20)
ax1.grid(ls='--', color='white', linewidth=0.4)
ax1.set_title('Posteriori Alpha')

# Parâmetro: beta
ax2.hist(beta, density=True, rwidth=0.9, bins=20)
ax2.grid(ls='--', color='white', linewidth=0.4)
ax2.set_title('Posteriori Beta')

# Parâmetro: sigma
ax3.hist(sigma, density=True, rwidth=0.9, bins=20)
ax3.grid(ls='--', color='white', linewidth=0.4)
ax3.set_title('Posteriori Sigma')

plt.show()

A observação da nossa Incerteza¶

Agora nós iremos fazer a amostragem das posterioris. Uma das razões para usarmos esse procedimento de amostragem é que isso torna mais fácil de se pensar e, também, esse procedimento pode ser aplicado a todos os tipos de modelos possíveis bayesianos que quisermos ajustar.

Nossa distribuição à posteriori contém muitas linhas retas lá dentro, representadas pela amostragem dos valores da distribuição dos \(\alpha\) e dos \(\beta\) com o \(\sigma\) informando o desvio padrão.

Abaixo vamos ver as primeiras linhas do conjunto de dados dessas estimativas e as suas medidas de resumo:

print(pd.DataFrame({'alpha': alpha[:10], 'beta': beta[:10], 'sigma': sigma[:10]}))

        alpha      beta     sigma
154.765819  0.946624  4.904993
154.725060  0.894531  4.897424
154.969708  0.940925  5.181750
154.966983  0.978020  4.871757
154.991493  0.892440  5.335516
154.752034  0.866076  5.027855
154.205608  0.862118  5.100277
154.324713  0.838174  5.273384
154.761308  0.952450  5.090631
154.352351  0.883649  5.340434

Note

Para cada linha do conjunto de dados acima é uma linha reta no nosso modelo!

e o resumo dos dados:

pd.DataFrame({'alpha': alpha[:10], 'beta': beta[:10], 'sigma': sigma[:10]}).describe()

	alpha	beta	sigma
count	10.000000	10.000000	10.000000
mean	154.681508	0.905501	5.102402
std	0.287348	0.046124	0.178393
min	154.205608	0.838174	4.871757
25%	154.445528	0.870469	4.935708
50%	154.756671	0.893486	5.095454
75%	154.916692	0.945199	5.250476
max	154.991493	0.978020	5.340434

Nós temos muitas linhas! E as linhas mais plausíveis são as linhas que tem maior número de maneiras de reproduzir os dados que observamos. Essa plausabilidade é apresentada para nós através do maior acúmulo de linhas retas que se sobrepõem.

Quanto mais as linhas se sobreporem entre si em uma certa região, maior será a plausabilidade dessas retas descreverem os dados que observamos.

Note

O acúmulo dessas retas representa a nossa incerteza.

Foi primeira vez que conseguimos realmente observar, conscientemente, a nossa incerteza em um gráfico. Lindo demais!!!

Para termos uma melhor compreensão nesse momento, nós vamos supor que ao invés de nossa amostra ter 251 indivíduos, vamos supor que temos apenas alguns subconjuntos. Nosso objetivo aqui é mostrar que quanto mais informações temos (isto é, quanto maior a nossa amostra), menor será a nossa incerteza!

Para ficar mais claro, vamos simular a essa nossa estratégia e ver as diferenças dos acúmulos das linhas.

def generate_parcial_stan_models_results(N):
    """
    Ajustando as estimativas de um modelo linear 
    usando os primeiros N indivíduos da amostra.
    
    Parameters:
        N: Quantidade de indivíduos que iremos 
           utilizar na estimativa dos parâmetros.
    
    Return: 
       Estimativas do alpha, beta e o sigma, dado os dados. 
    """
    weight_adjust = weight - weight.mean() 

    my_data = {
        'N': N,
        'height': height[:N],
        'weight': weight_adjust[:N],
    }

    posteriori = stan.build(model_stan, data=my_data)
    fit = posteriori.sample(num_chains=4, num_samples=1000)

    # Recuperando os parâmetros
    alpha = fit['alpha'].flatten()
    beta = fit['beta'].flatten()
    sigma = fit['sigma'].flatten()
    
    return alpha, beta, sigma

# =====================================================
#   Rodando o modelo anterior utilizando diferentes N
# =====================================================

N_10 = 10     
alpha_10, beta_10, sigma_10 = generate_parcial_stan_models_results(N_10);


N_50 = 50     
alpha_50, beta_50, sigma_50 = generate_parcial_stan_models_results(N_50);


N_100 = 100     
alpha_100, beta_100, sigma_100 = generate_parcial_stan_models_results(N_100);


N = len(weight) # Com toda amostra disponível
alpha, beta, sigma = generate_parcial_stan_models_results(N);

Note

Uma das vantagens do gráfico gráfico de linhas (gráfico de espaguetes) é deixar claro que os limites formados pelas retas não tem significado algum.

# =============================
#     Gráfico de Espaguete
# =============================
# Vamos usar a variável weight_adjust igual para todos os 
# pesos sem perda de generalidade para essa análise.

# Plot dos dados altura x peso
fig, [[ax1, ax2], [ax3, ax4]] = plt.subplots(2, 2, figsize=(19, 13))

# Número de linhas retas que iremos plotar
qty_lines = 30

# ===============================================
#    Estimando as curvas usando 10 amostras
# ===============================================

ax1.scatter(weight[:N_10], height[:N_10])

for i in range(qty_lines):
    ax1.plot(weight, alpha_10[i] + beta_10[i] * weight_adjust, color='gray', linewidth=0.3)

ax1.grid(ls='--', color='white', linewidth=0.4)
ax1.set_title('Estimativa com N = ' + str(N_10) + ' pontos amostrados')
ax1.set_xlabel('Peso (Weight)')
ax1.set_ylabel('Altura (Height)')


# ===============================================
#    Estimando as curvas usando 50 amostras
# ===============================================

ax2.scatter(weight[:N_50], height[:N_50])

for i in range(qty_lines):
    ax2.plot(weight, alpha_50[i] + beta_50[i] * weight_adjust, color='gray', linewidth=0.3)

ax2.grid(ls='--', color='white', linewidth=0.4)
ax2.set_title('Estimativa com N = ' + str(N_50) + ' pontos amostrados')
ax2.set_xlabel('Peso (Weight)')
ax2.set_ylabel('Altura (Height)')


# ===============================================
#    Estimando as curvas usando 100 amostras
# ===============================================

ax3.scatter(weight[:N_100], height[:N_100])

for i in range(qty_lines):
    ax3.plot(weight, alpha_100[i] + beta_100[i] * weight_adjust, color='gray', linewidth=0.3)

ax3.grid(ls='--', color='white', linewidth=0.4)
ax3.set_title('Estimativa com N = ' + str(N_100) + ' pontos amostrados')
ax3.set_xlabel('Peso (Weight)')
ax3.set_ylabel('Altura (Height)')


# =================================================
#    Estimando as curvas usando todas as amostras
# =================================================

ax4.scatter(weight, height)

for i in range(qty_lines):
    ax4.plot(weight, alpha[i] + beta[i] * weight_adjust, color='gray', linewidth=0.3)

ax4.grid(ls='--', color='white', linewidth=0.4)
ax4.set_title('Estimativa com N = ' + str(N) + ' (todos) pontos amostrados')
ax4.set_xlabel('Peso (Weight)')
ax4.set_ylabel('Altura (Height)')

plt.show()

A divisão dos gráficos acima foi construída para conseguirmos percerber, de modo visual, que quando mais amostras tivermos coletado (isto é, quanto maior número de pontos azuis nós tivermos coletados como amostra), mais informações teremos e, portanto, teremos muito menos incerteza sobre o nosso objeto de estudo.

A nossa incerteza pode ser observada pela dispersão das curvas no gráfico, portanto quanto maior for a dispersão maior será a nossa incerteza sobre o que está acontecendo.

Construindo a Distribuição Preditiva de \(\mu_i\), dado um peso \(x_i\)¶

Modelo linear para a altura média \(\mu_i\) é:

\[ \mu_i = \alpha + \beta(x_i - \bar{x}) \]

Agora, a ideia básica é que, dado um valor específico do peso de um indivíduo, \(x_i\) (weight), podemos obter uma distribuição preditiva do \(\mu_i\). Essa distribuição preditiva nos informa quais as regiões de maior confiança que podemos esperar para a média da altura, de uma pessoa com o peso \(x_i\).

Para exemplificar, vamos supor que queremos estimar a altura (height) de uma pessoa com \(50 kg\), então:

\[ \mu_i = \alpha + \beta(50 - \bar{x}) \]

# ===========================================================
#    Construindo a distribuição preditiva de 𝜇 | x=50
# ===========================================================

mu_50_kg = alpha + beta * (50 - weight.mean())

plt.figure(figsize=(17, 9))
plt.hist(mu_50_kg, density=False, bins=80, rwidth=0.8)
plt.grid(ls='--', color='white', linewidth=0.3)

plt.title('Distribuição de $\mu_i$ para o peso de $50 kg$')
plt.ylabel('Plausabilidade')
plt.xlabel('$\mu | peso=50$')

plt.show()

Calculando a predição para todos os \(\mu_i\)¶

def HPDI(posterior_samples, credible_mass):
    
    # Calcula o maior intervalo de probabilidades a partir de uma amostra
    
    # Fonte: https://stackoverflow.com/questions/22284502/highest-posterior-density-region-and-central-credible-region
    # ** Refazer essa função para entender **
    
    sorted_points = sorted(posterior_samples)
    ciIdxInc = np.ceil(credible_mass * len(sorted_points)).astype('int')
    nCIs = len(sorted_points) - ciIdxInc
    ciWidth = [0]*nCIs
    
    for i in range(0, nCIs):
        ciWidth[i] = sorted_points[i + ciIdxInc] - sorted_points[i]
        HDImin = sorted_points[ciWidth.index(min(ciWidth))]
        HDImax = sorted_points[ciWidth.index(min(ciWidth))+ciIdxInc]

    return(HDImin, HDImax)

# =====================================
#    Calculando o HPDI dos dados
# =====================================

posterioris_dict_10 = {}  # Dicionário com todas as posterioris para cada um dos pesos (weight) variando de 25 á 70.
posterioris_dict_50 = {}  # Dicionário com todas as posterioris para cada um dos pesos (weight) variando de 25 á 70.
posterioris_dict_100 = {}  # Dicionário com todas as posterioris para cada um dos pesos (weight) variando de 25 á 70.
posterioris_dict = {}  # Dicionário com todas as posterioris para cada um dos pesos (weight) variando de 25 á 70.

for weight_i in range(25, 71):
    posterioris_dict_10[weight_i] = alpha_10 + beta_10 * (weight_i - weight.mean())
    posterioris_dict_50[weight_i] = alpha_50 + beta_50 * (weight_i - weight.mean())
    posterioris_dict_100[weight_i] = alpha_100 + beta_100 * (weight_i - weight.mean())
    posterioris_dict[weight_i] = alpha + beta * (weight_i - weight.mean())
    
# Gerando os dataframes
posterioris_10 = pd.DataFrame(posterioris_dict_10)
posterioris_50 = pd.DataFrame(posterioris_dict_50)
posterioris_100 = pd.DataFrame(posterioris_dict_100)
posterioris = pd.DataFrame(posterioris_dict)

# Calculando as médias das posterioris
posterioris_means_10 = posterioris_10.mean()
posterioris_means_50 = posterioris_50.mean()
posterioris_means_100 = posterioris_100.mean()
posterioris_means = posterioris.mean()

# Calculando os intervalos de HPDI das posterioris
posterioris_HPDIs_10 = []
posterioris_HPDIs_50 = []
posterioris_HPDIs_100 = []
posterioris_HPDIs = []

for weight_i in range(25, 71):
    posterioris_HPDIs_10.append(HPDI(posterioris_10[weight_i], 0.89))
    posterioris_HPDIs_50.append(HPDI(posterioris_50[weight_i], 0.89))
    posterioris_HPDIs_100.append(HPDI(posterioris_100[weight_i], 0.89))
    posterioris_HPDIs.append(HPDI(posterioris[weight_i], 0.89))
    
# Tranformando os dados em um array numpy    
posterioris_HPDIs_10 = np.array(posterioris_HPDIs_10)  
posterioris_HPDIs_50 = np.array(posterioris_HPDIs_50) 
posterioris_HPDIs_100 = np.array(posterioris_HPDIs_100)
posterioris_HPDIs = np.array(posterioris_HPDIs)

A Compatibilidade da Gravata Borboleta¶

Conforme os quatro gráficos plotados acima, iremos replotá-los porém agora mostrando o intervalo de compatibilidade (estilo gravata borboleta).

Ao olhar para esse gráfico é fácil cairmos na tentação de acharmos que os limites que escolhemos significa alguma coisa, ele não tem nenhum significado a não ser os limites de corte pelo HPDI que escolhemos.

# ==================================================
#       Plotando os gráfico de gravata borboletas
# ==================================================

fig, [[ax1, ax2], [ax3, ax4]] = plt.subplots(2, 2, figsize=(19,12))

# ===============================================
#    Estimando as curvas usando 10 amostras
# ===============================================
ax1.plot(range(25, 71), posterioris_HPDIs_10, color='darkgray', linewidth=0.5)
ax1.fill_between(range(25, 71), posterioris_HPDIs_10[:, 0], posterioris_HPDIs_10[:, 1], color='gray', alpha=0.4)
ax1.scatter(weight[:N_10], height[:N_10], alpha=0.5)

ax1.plot(range(25, 71), posterioris_means_10, color='black', linewidth=1)
ax1.grid(ls='--', color='white', linewidth=0.3)
ax1.set_title('Intervalos HPDI da posteriori da média altura com N = ' + str(N_10) + ' pontos de amostras')
ax1.set_ylabel('Altura (height)')
ax1.set_xlabel('Peso (weight)')

# ===============================================
#    Estimando as curvas usando 50 amostras
# ===============================================
ax2.plot(range(25, 71), posterioris_HPDIs_50, color='darkgray', linewidth=0.5)
ax2.fill_between(range(25, 71), posterioris_HPDIs_50[:, 0], posterioris_HPDIs_50[:, 1], color='gray', alpha=0.4)
ax2.scatter(weight[:N_50], height[:N_50], alpha=0.5)

ax2.plot(range(25, 71), posterioris_means_50, color='black', linewidth=1)
ax2.grid(ls='--', color='white', linewidth=0.3)
ax2.set_title('Intervalos HPDI da posteriori da média alturacom N = ' + str(N_50) + ' pontos de amostras')
ax2.set_ylabel('Altura (height)')
ax2.set_xlabel('Peso (weight)')

# ===============================================
#    Estimando as curvas usando 100 amostras
# ===============================================
ax3.plot(range(25, 71), posterioris_HPDIs_100, color='darkgray', linewidth=0.5)
ax3.fill_between(range(25, 71), posterioris_HPDIs_100[:, 0], posterioris_HPDIs_100[:, 1], color='gray', alpha=0.4)
ax3.scatter(weight[:N_100], height[:N_100], alpha=0.5)

ax3.plot(range(25, 71), posterioris_means_100, color='black', linewidth=1)
ax3.grid(ls='--', color='white', linewidth=0.3)
ax3.set_title('Intervalos HPDI da posteriori da média altura com N = ' + str(N_100) + ' pontos de amostras')
ax3.set_ylabel('Altura (height)')
ax3.set_xlabel('Peso (weight)')

# ==================================================
#    Estimando as curvas usando todas as amostras
# ==================================================
ax4.plot(range(25, 71), posterioris_HPDIs, color='darkgray', linewidth=0.5)
ax4.fill_between(range(25, 71), posterioris_HPDIs[:, 0], posterioris_HPDIs[:, 1], color='gray', alpha=0.4)
ax4.scatter(weight, height, alpha=0.2)

ax4.plot(range(25, 71), posterioris_means, color='black', linewidth=1)
ax4.grid(ls='--', color='white', linewidth=0.5)
ax4.set_title('Intervalos HPDI da posteriori da média altura com todos os pontos de amostras')
ax4.set_ylabel('Altura (height)')
ax4.set_xlabel('Peso (weight)')

plt.show()

Perceba que nos gráficos de borboletas nós usamos os cortes para os intervalos com \(89\%\). Isso é apenas para mostrar que não, isso mesmo, não existe nada especial nos cortes do HPDI.

Porém existem alguns cortes que são mais usados como o de \(50%\), \(90%\) ou \(95%\), corriqueiramente são usados em outras escolas estatísticas, são apenas cortes que queremos visualizar. Nada mais!

Até agora nós fizemos os gráficos apenas usando os parâmetros da média (\(\mu_i\)), mas podemos construir o gráfico de envelope para o sigma (\(\sigma\)).

Iremos construir ambos utilizando toda a nossa amostra.

# ============================================================================
#     Plotando o HPDI da distribuição preditiva da alutra e de sua média
# ============================================================================

# ===========================================
# Nota: Altere o valor para ver a diferença.
# ===========================================
HPDI_range = 0.93  # Define o tamanho do intervalo do HPDI. 

posterioris_height_HPDIs = []

for weight_i in range(25, 71):
    posteiori_height = np.random.normal(alpha + beta*(weight_i - weight.mean()), sigma)
    posterioris_height_HPDIs.append(HPDI(posteiori_height, HPDI_range))
    
posterioris_height_HPDIs = np.array(posterioris_height_HPDIs)

# Plontando o gráfico
plt.figure(figsize=(17, 9))

plt.plot(range(25, 71), posterioris_height_HPDIs, color='darkgray', linewidth=0.5)
plt.fill_between(range(25, 71), posterioris_height_HPDIs[:, 0], posterioris_height_HPDIs[:, 1], color='blue', alpha=0.4)
plt.fill_between(range(25, 71), posterioris_HPDIs[:, 0], posterioris_HPDIs[:, 1], color='blue', alpha=0.4)
plt.scatter(weight, height, color='red', alpha=0.4)

plt.plot(range(25, 71), posterioris_means, color='black', linewidth=1)
plt.grid(ls='--', color='white', linewidth=0.3)
plt.title('Intervalos HPDI da posteriori da altura - Usando  HPDI  com ' + str(int(HPDI_range*100)) + '% de massa de probabilidade')
plt.ylabel('Altura (height)')
plt.xlabel('Peso (weight)')

plt.show()

Como podemos ver no gráfico acima, temos duas regiões de confiança. A região mais interna - azul mais escura - nos indica que esse trecho contém \(89\%\) de probabilidade da média estar ali dentro, para cada um dos pesos.

Já a região mais externa - aquela faixa azul mais clara - nos indica que, com \(93\%\) probabilidade, será a região mais provável da altura de um pessoa, dado um determinado peso.

Como teste para uma melhor compreensão, faça alterações na variável HPDI_range para valores bem altos, por exemplo, \(0.99\) e veja o quanto a faixa mais externa cobrirá os dados (os pontos vermelhos). Para valores menores, isso também ocorre, por exemplo, \(0.7\).

Verifique!

Warning

A utilização de valores baixos, irá demandar cálculos mais complexos e por isso pode ser muito demorado.

Curvas a partir das Linhas Retas¶

O interessante de modelos lineares é que eles não são apenas lineares, podem ser curvas também. Isso é uma coisa meio enlouquecedora sobre o termo que é convencionalmente usado para desenhar linhas retas em gráficos.

Mas uma regressão linear é aditiva, vamos ter que essa equação para média será a soma de alguns parâmetros vezes algumas variáveis observadas. Assim teremos a soma de um monte de termos.

Temos uma equação aditiva e coisas aditivas são lineares na matemática.

Mas para nós, as palavras “adtiva” e “linear” são diferentes, então chamaremos esses modelos de Regressões Aditivas. Porque, podemos usar coisas que não parecem linhas retas. E, de agora até o final do capítulo, iremos fazer isso, desenhar curvas a partir de linhas retas.

Mas antes de continuar, gostaria de explicar o por que é interessante aprendermos a fazer esse tipo de modelagem?

Simples! Por que a Natureza não se limita a se manifestar por relações lineares entre duas variáveis! Os nossos intervalos podem ser aproximações úteis, mas geralmente são bobos.

No nosso caso, nós usamos as variáveis \(altura\) e \(peso\), porém apenas para indivíduos que tivessem \(30\) anos ou mais. Se fossêmos usar todos os indivíduos a relação do \(peso\) e a \(altura\) não mais seria uma linha reta, mas sim uma curva!

plt.figure(figsize=(17, 9))

# Indivíduos com menos de 18 anos de idade
height_0_18 = df.loc[df.age < 18, 'height'].values
weight_0_18 = df.loc[df.age < 18, 'weight'].values

# Indivíduos maiores (ou iguais) de 18 anos de idade - Adultos
height_18_ = df.loc[df.age >= 18, 'height'].values
weight_18_ = df.loc[df.age >= 18, 'weight'].values

# Plotando os gráficos
plt.scatter(weight_0_18, height_0_18, label='De 0 a 17 anos de idade', color='red', alpha=0.4)
plt.scatter(weight_18_, height_18_, label='(Adultos) 18+ anos de idade', color='blue', alpha=0.4)

plt.grid(ls='--', color='white', linewidth=0.3)
plt.legend()
plt.title('Altura x peso | Todos os indivíduos da amostra')
plt.xlabel('Peso')
plt.ylabel('Altura')

plt.show()

Vamos ver dois modos de construir os modelos de Regressão Aditivas:

Regressão Polinomial
- Usada de forma comum
- Geralmente apresenta um comportamento bem ruim

Esse é um tipo bem comum de regressão, mas também é um modelo bem ruim. Não há nada de errado usar esse tipo de regressão supondo que você entenda o que golem está fazendo e assim usar de forma responsável. Geralmente se utilizam esse tipo de modelo sem uma devida atenção. Iremos ver o motivo dele serem mal comportados em breve.

Splines
- São muito flexíveis
- Altamente Geocêntrica

Existem muitos tipos de splines, mas nós iremos ver apenas as splines de base, que são provavelmente as mais comuns. Software de desenhos como GIMP, Blender e etc. que possuem essa ferramenta com o nome de curvas de Bezier, que nada mais são do que splines de base.

As Splines são muito flexíveis, muito mais flexíveis do que as Regressões Polinomiais, porém não apresentam a patologia que os polinômios apresentam e, por isso, são considerados, na maioria dos casos, melhores que as regressões polinômiais.

Porém, tanto as splines quantos as regressões polinomiais são estratégias geocêntricas, não há nada de científico nelas. São apenas aproximações como o * modelo de epiciclos de Ptolomeu*.

Extrapolações com esses modelos podem trazer desastres nas predições, por isso é necessário checar e entender o que está acontecendo com o modelo.

	year	doy	temp	temp_upper	temp_lower
count	1215.000000	827.000000	1124.000000	1124.000000	1124.000000
mean	1408.000000	104.540508	6.141886	7.185151	5.098941
std	350.884596	6.407036	0.663648	0.992921	0.850350
min	801.000000	86.000000	4.670000	5.450000	0.750000
25%	1104.500000	100.000000	5.700000	6.480000	4.610000
50%	1408.000000	105.000000	6.100000	7.040000	5.145000
75%	1711.500000	109.000000	6.530000	7.720000	5.542500
max	2015.000000	124.000000	8.300000	12.100000	7.740000

Statistical Rethinking

4 - Funções Wiggly¶

A observação da nossa Incerteza¶

Construindo a Distribuição Preditiva de \(\mu_i\), dado um peso \(x_i\)¶

Calculando a predição para todos os \(\mu_i\)¶

A Compatibilidade da Gravata Borboleta¶

Curvas a partir das Linhas Retas¶

Regressão Polinomial¶

Modelo Parabólico da Altura¶

Padronizar os Preditores¶

Modelos Cúbicos¶

Modelo Regressão Cúbico para a Altura¶

Dores dos Polinômios¶

Splines¶

Agindo Localmente - B-splines¶

Festival das Flores de Cerejeira¶

Como as splines funcionam?¶

Knots - Como escolher os pontos?¶

Possibilidades das Splines¶