Arquivo de ciência-de-dados

Prever o futuro do planeta: como processos gaussianos modelam as mudanças climáticas

19/12/202531/10/2025 Por antonino

Imagine que você é um cientista climático analisando décadas de dados de CO₂ do observatório de Mauna Loa no Havaí. Você vê padrões complexos: uma tendência de longo prazo de aumento, flutuações sazonais anuais e variações aleatórias. Prever como evoluirá essa curva vital para nosso planeta requer um método que capture todos esses padrões simultaneamente. O Processo Gaussiano é a ferramenta perfeita para esta tarefa, conseguindo modelar tendências, sazonalidades e incertezas em uma única framework coerente.

Como isso funciona na prática?

O exemplo de Mauna Loa demonstra o poder dos processos gaussianos para modelar séries temporais complexas com múltiplos componentes. O segredo está na combinação inteligente de kernels que capturam diferentes aspectos dos dados: um kernel para a tendência de longo prazo, outro para a sazonalidade anual e um terceiro para variações residuais. Diferentemente de métodos que exigem decomposição manual dos componentes, o GPR aprende automaticamente a importância relativa de cada padrão. Esta abordagem permite não apenas prever valores futuros, mas quantificar realisticamente a incerteza dessas previsões – informação crucial para políticas climáticas.

Mãos na massa: modelando CO₂ com kernels compostos

"""
Modelagem de dados de CO₂ de Mauna Loa com Processos Gaussianos
Demonstra como kernels compostos capturam tendências complexas
"""

import numpy as np
import matplotlib.pyplot as plt
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import RBF, ExpSineSquared, RationalQuadratic, WhiteKernel

# Carregando dados de CO₂ (simulados para exemplo)
# Na prática, você usaria dados reais de https://www.esrl.noaa.gov/gmd/ccgg/trends/
def carregar_dados_co2_simulados():
    """Simula dados similares aos de Mauna Loa"""
    np.random.seed(42)
    # Período de 1990-2020 em meses
    anos = np.linspace(1990, 2020, 30*12)
    
    # Componentes: tendência + sazonalidade + ruído
    tendencia = 350 + 2.0 * (anos - 1990)  # Aumento linear
    sazonalidade = 3.0 * np.sin(2 * np.pi * (anos - 1990))  # Sazonalidade anual
    ruido = np.random.normal(0, 0.5, len(anos))  # Ruído de medição
    
    co2 = tendencia + sazonalidade + ruido
    return anos.reshape(-1, 1), co2

X, y = carregar_dados_co2_simulados()

print(f"Dados de CO₂ carregados: {len(X)} meses de medições")
print(f"Variação de CO₂: {y.min():.1f} a {y.max():.1f} ppm")

# Kernel composto para capturar diferentes componentes temporais
kernel = (
    1.0 * RBF(length_scale=50.0) +                    # Tendência de longo prazo
    1.0 * ExpSineSquared(length_scale=1.0, periodicity=1.0) +  # Sazonalidade anual
    1.0 * RationalQuadratic(length_scale=1.0, alpha=1.0) +     # Variações de médio prazo
    1.0 * WhiteKernel(noise_level=0.1)               # Ruído de medição
)

# Criando e treinando o GPR
gpr = GaussianProcessRegressor(
    kernel=kernel,
    alpha=0.0,
    normalize_y=True,
    n_restarts_optimizer=9
)

# Dividindo em treino (até 2015) e teste (após 2015)
mask_treino = X.ravel() <= 2015
X_treino, y_treino = X[mask_treino], y[mask_treino]
X_teste, y_teste = X[~mask_treino], y[~mask_treino]

print(f"\nTreinando com {len(X_treino)} pontos (até 2015)")
print(f"Testando com {len(X_teste)} pontos (2016-2020)")

gpr.fit(X_treino, y_treino)

print(f"\nKernel otimizado: {gpr.kernel_}")

# Fazendo previsões para o futuro
X_futuro = np.linspace(2021, 2030, 120).reshape(-1, 1)
y_pred, sigma = gpr.predict(X_futuro, return_std=True)

# Visualizando resultados
plt.figure(figsize=(14, 10))

# Plot principal
plt.subplot(2, 1, 1)
plt.plot(X_treino, y_treino, 'b.', alpha=0.6, label='Dados de treino (1990-2015)')
plt.plot(X_teste, y_teste, 'r.', alpha=0.6, label='Dados de teste (2016-2020)')
plt.plot(X_futuro, y_pred, 'g-', linewidth=2, label='Previsão GPR (2021-2030)')
plt.fill_between(X_futuro.ravel(), 
                y_pred - 1.96*sigma, 
                y_pred + 1.96*sigma, 
                alpha=0.3, color='green', label='95% intervalo de confiança')

plt.xlabel('Ano')
plt.ylabel('CO₂ (ppm)')
plt.title('Previsão de CO₂ com Processo Gaussiano - Mauna Loa Simulado')
plt.legend()
plt.grid(True, alpha=0.3)

# Zoom na área de previsão
plt.subplot(2, 1, 2)
zoom_mask = (X.ravel() >= 2010) | (X_futuro.ravel() <= 2030)
X_zoom = np.concatenate([X[X.ravel() >= 2010], X_futuro])
y_zoom_actual = np.concatenate([y[X.ravel() >= 2010], [np.nan] * len(X_futuro)])
y_zoom_pred = np.concatenate([[np.nan] * sum(X.ravel() >= 2010), y_pred])

plt.plot(X_zoom, y_zoom_actual, 'ro-', alpha=0.7, label='Dados observados')
plt.plot(X_zoom, y_zoom_pred, 'g-', linewidth=2, label='Previsão GPR')
plt.fill_between(X_futuro.ravel(), 
                y_pred - 1.96*sigma, 
                y_pred + 1.96*sigma, 
                alpha=0.3, color='green', label='Incerteza')

plt.xlabel('Ano')
plt.ylabel('CO₂ (ppm)')
plt.title('Zoom: Previsões 2021-2030 e Incerteza')
plt.legend()
plt.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

# Análise das previsões
print(f"\nPrevisões para 2030:")
print(f"CO₂ esperado: {y_pred[-1]:.1f} ppm")
print(f"Intervalo de 95% confiança: {y_pred[-1] - 1.96*sigma[-1]:.1f} a {y_pred[-1] + 1.96*sigma[-1]:.1f} ppm")

100

101

102

103

104

105

106

107

108

109

"""

Modelagem de dados de CO₂ de Mauna Loa com Processos Gaussianos

Demonstra como kernels compostos capturam tendências complexas

"""

import numpy as np

import matplotlib.pyplot as plt

from sklearn.gaussian_process import GaussianProcessRegressor

from sklearn.gaussian_process.kernels import RBF, ExpSineSquared, RationalQuadratic, WhiteKernel

# Carregando dados de CO₂ (simulados para exemplo)

# Na prática, você usaria dados reais de https://www.esrl.noaa.gov/gmd/ccgg/trends/

def carregar_dados_co2_simulados():

"""Simula dados similares aos de Mauna Loa"""

np.random.seed(42)

# Período de 1990-2020 em meses

anos = np.linspace(1990, 2020, 30*12)

# Componentes: tendência + sazonalidade + ruído

tendencia = 350 + 2.0 * (anos - 1990) # Aumento linear

sazonalidade = 3.0 * np.sin(2 * np.pi * (anos - 1990)) # Sazonalidade anual

ruido = np.random.normal(0, 0.5, len(anos)) # Ruído de medição

co2 = tendencia + sazonalidade + ruido

return anos.reshape(-1, 1), co2

X, y = carregar_dados_co2_simulados()

print(f"Dados de CO₂ carregados: {len(X)} meses de medições")

print(f"Variação de CO₂: {y.min():.1f} a {y.max():.1f} ppm")

# Kernel composto para capturar diferentes componentes temporais

kernel = (

1.0 * RBF(length_scale=50.0) + # Tendência de longo prazo

1.0 * ExpSineSquared(length_scale=1.0, periodicity=1.0) + # Sazonalidade anual

1.0 * RationalQuadratic(length_scale=1.0, alpha=1.0) + # Variações de médio prazo

1.0 * WhiteKernel(noise_level=0.1) # Ruído de medição

)

# Criando e treinando o GPR

gpr = GaussianProcessRegressor(

kernel=kernel,

alpha=0.0,

normalize_y=True,

n_restarts_optimizer=9

)

# Dividindo em treino (até 2015) e teste (após 2015)

mask_treino = X.ravel() <= 2015

X_treino, y_treino = X[mask_treino], y[mask_treino]

X_teste, y_teste = X[~mask_treino], y[~mask_treino]

print(f"\nTreinando com {len(X_treino)} pontos (até 2015)")

print(f"Testando com {len(X_teste)} pontos (2016-2020)")

gpr.fit(X_treino, y_treino)

print(f"\nKernel otimizado: {gpr.kernel_}")

# Fazendo previsões para o futuro

X_futuro = np.linspace(2021, 2030, 120).reshape(-1, 1)

y_pred, sigma = gpr.predict(X_futuro, return_std=True)

# Visualizando resultados

plt.figure(figsize=(14, 10))

# Plot principal

plt.subplot(2, 1, 1)

plt.plot(X_treino, y_treino, 'b.', alpha=0.6, label='Dados de treino (1990-2015)')

plt.plot(X_teste, y_teste, 'r.', alpha=0.6, label='Dados de teste (2016-2020)')

plt.plot(X_futuro, y_pred, 'g-', linewidth=2, label='Previsão GPR (2021-2030)')

plt.fill_between(X_futuro.ravel(),

y_pred - 1.96*sigma,

y_pred + 1.96*sigma,

alpha=0.3, color='green', label='95% intervalo de confiança')

plt.xlabel('Ano')

plt.ylabel('CO₂ (ppm)')

plt.title('Previsão de CO₂ com Processo Gaussiano - Mauna Loa Simulado')

plt.legend()

plt.grid(True, alpha=0.3)

# Zoom na área de previsão

plt.subplot(2, 1, 2)

zoom_mask = (X.ravel() >= 2010) | (X_futuro.ravel() <= 2030)

X_zoom = np.concatenate([X[X.ravel() >= 2010], X_futuro])

y_zoom_actual = np.concatenate([y[X.ravel() >= 2010], [np.nan] * len(X_futuro)])

y_zoom_pred = np.concatenate([[np.nan] * sum(X.ravel() >= 2010), y_pred])

plt.plot(X_zoom, y_zoom_actual, 'ro-', alpha=0.7, label='Dados observados')

plt.plot(X_zoom, y_zoom_pred, 'g-', linewidth=2, label='Previsão GPR')

plt.fill_between(X_futuro.ravel(),

y_pred - 1.96*sigma,

y_pred + 1.96*sigma,

alpha=0.3, color='green', label='Incerteza')

plt.xlabel('Ano')

plt.ylabel('CO₂ (ppm)')

plt.title('Zoom: Previsões 2021-2030 e Incerteza')

plt.legend()

plt.grid(True, alpha=0.3)

plt.tight_layout()

plt.show()

# Análise das previsões

print(f"\nPrevisões para 2030:")

print(f"CO₂ esperado: {y_pred[-1]:.1f} ppm")

print(f"Intervalo de 95% confiança: {y_pred[-1] - 1.96*sigma[-1]:.1f} a {y_pred[-1] + 1.96*sigma[-1]:.1f} ppm")

Os detalhes que fazem diferença

O poder do GPR no exemplo de Mauna Loa está na combinação estratégica de kernels. O kernel ExpSineSquared é particularmente importante pois modela padrões periódicos como a sazonalidade anual do CO₂ (plantas absorvem CO₂ no verão, liberam no inverno). Contudo, a escolha dos hiperparâmetros iniciais é crucial – valores muito distantes do ótimo podem fazer a otimização falhar. Analogamente importante é entender que a incerteza cresce conforme previsões se afastam no futuro, refletindo honestamente nossa ignorância crescente. O kernel RationalQuadratic captura variações de médio prazo que não são nem totalmente periódicas nem totalmente aleatórias.

ExpSineSquared: Ideal para padrões periódicos como sazonalidade anual
RBF: Captura tendências suaves de longo prazo
RationalQuadratic: Flexível para padrões de médio prazo
WhiteKernel: Modela ruído de medição e variações não explicadas

Perguntas que os iniciantes fazem

Você deve estar se perguntando: “Por que usar tantos kernels diferentes?” Excelente questão! Cada kernel captura um tipo diferente de padrão nos dados. Usar apenas um kernel RBF, por exemplo, não conseguiria modelar adequadamente a sazonalidade anual. Uma confusão comum é pensar que mais kernels sempre significam melhor performance – na verdade, kernels desnecessários podem levar a overfitting. Outra dúvida frequente: “Como o GPR sabe qual parte é tendência e qual é sazonalidade?” A otimização dos hiperparâmetros automaticamente aprende a importância relativa de cada componente através da maximização da verossimilhança marginal.

Para onde ir agora?

Experimente aplicar GPR com kernels compostos em suas próprias séries temporais. Comece com dados climáticos, financeiros ou de negócios que exibam múltiplos padrões. Pratique a interpretação dos intervalos de confiança – eles não são apenas técnicos, mas contam uma história sobre o que sabemos e não sabemos. O momento “aha!” acontece quando você vê o modelo capturando automaticamente padrões complexos que você só conseguia ver manualmente antes.

Assuntos relacionados

Para dominar esta aplicação, estude:

Séries temporais: decomposição, estacionariedade e sazonalidade
Análise espectral: identificação de componentes periódicos
Kernels para séries temporais: propriedades e combinações
Mudanças climáticas: ciclos do carbono e dinâmica atmosférica
Incerteza em previsões: interpretação e comunicação de riscos

Referências que valem a pena

Modelos Lineares Generalizados: Busca por correspondência ortogonal

19/12/202518/10/2025 Por antonino

Anteriormente exploramos algoritmos como LARS e Lasso para problemas de regressão esparsa. Analogamente, a Busca por Correspondência Ortogonal (OMP) é outro algoritmo eficiente para aproximação esparsa que seleciona iterativamente as features mais correlacionadas com o resíduo atual.

Conceito Fundamental do OMP

Primordialmente, o OMP é um algoritmo guloso que resolve problemas de aproximação esparsa. Decerto, ele busca encontrar uma representação esparsa dos dados usando um número limitado de features (átomos) de um dicionário.

Conforme a documentação do scikit-learn, o OMP é particularmente útil quando sabemos antecipadamente o número de features que desejamos selecionar. Similarmente ao LARS, ele constrói a solução de forma incremental, mas com uma abordagem de projeção ortogonal.

Algoritmo OMP

O algoritmo opera através dos seguintes passos iterativos:

Inicializar o resíduo com o target original
Encontrar a feature mais correlacionada com o resíduo atual
Adicionar essa feature ao conjunto ativo
Resolver o problema de mínimos quadrados usando apenas as features ativas
Atualizar o resíduo subtraindo a contribuição das features selecionadas
Repetir até atingir o critério de parada

Formulação Matemática

O objetivo do OMP é resolver:

\(\min_{w} ||Xw – y||_2^2\)

Sujeito a:

\(||w||_0 \leq k\)

Onde:

X é a matriz de features
y é o vetor target
w são os coeficientes esparsos
k é o número máximo de features não-zero
||w||₀ é a norma L0 (número de elementos não-zero)

Implementações no Scikit-learn

Atualmente, o scikit-learn oferece duas implementações principais:

OrthogonalMatchingPursuit: Implementação padrão do OMP
OrthogonalMatchingPursuitCV: Versão com validação cruzada para seleção automática do parâmetro n_nonzero_coefs

Parâmetros Principais

Os principais parâmetros para ajuste no OMP são:

n_nonzero_coefs: Número máximo de coeficientes não-zero
tol: Tolerância para erro de aproximação
fit_intercept: Se deve calcular intercept
normalize: Se deve normalizar as features

Exemplo Prático: OMP em Ação

Ademais, vejamos um exemplo completo demonstrando o uso do Orthogonal Matching Pursuit:

'''
Aplicação do Orthogonal Matching Pursuit (OMP)
para Regressão com Seleção de Features
CÓDIGO CORRIGIDO - Problema com OrthogonalMatchingPursuitCV
'''

import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets import make_regression
from sklearn.linear_model import OrthogonalMatchingPursuit
from sklearn.linear_model import Lasso, LassoLars
from sklearn.model_selection import train_test_split, cross_val_score
from sklearn.preprocessing import StandardScaler
from sklearn.metrics import mean_squared_error, r2_score

# Gerar dataset com features esparsas
n_samples = 200
n_features = 50
n_informative = 10

print("Configuração do Dataset:")
print(f"Número de amostras: {n_samples}")
print(f"Número de features: {n_features}")
print(f"Features informativas: {n_informative}")

# Gerar dados com estrutura esparsa
X, y, true_coef = make_regression(n_samples=n_samples, n_features=n_features,
                                 n_informative=n_informative, noise=20,
                                 coef=True, random_state=42)

print(f"\nDimensões: X {X.shape}, y {y.shape}")

# Dividir em treino e teste
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3,
                                                    random_state=42)

# Normalizar os dados
scaler = StandardScaler()
X_train_scaled = scaler.fit_transform(X_train)
X_test_scaled = scaler.transform(X_test)

print(f"\n--- Comparação: OMP vs Outros Algoritmos Esparsos ---")

# Orthogonal Matching Pursuit com diferentes números de features
omp_5 = OrthogonalMatchingPursuit(n_nonzero_coefs=5)
omp_5.fit(X_train_scaled, y_train)
y_pred_omp_5 = omp_5.predict(X_test_scaled)
mse_omp_5 = mean_squared_error(y_test, y_pred_omp_5)
r2_omp_5 = r2_score(y_test, y_pred_omp_5)

omp_10 = OrthogonalMatchingPursuit(n_nonzero_coefs=10)
omp_10.fit(X_train_scaled, y_train)
y_pred_omp_10 = omp_10.predict(X_test_scaled)
mse_omp_10 = mean_squared_error(y_test, y_pred_omp_10)
r2_omp_10 = r2_score(y_test, y_pred_omp_10)

# Lasso para comparação
lasso = Lasso(alpha=0.1, max_iter=10000, random_state=42)
lasso.fit(X_train_scaled, y_train)
y_pred_lasso = lasso.predict(X_test_scaled)
mse_lasso = mean_squared_error(y_test, y_pred_lasso)
r2_lasso = r2_score(y_test, y_pred_lasso)

# LassoLars para comparação
lassolars = LassoLars(alpha=0.1, max_iter=1000)
lassolars.fit(X_train_scaled, y_train)
y_pred_lassolars = lassolars.predict(X_test_scaled)
mse_lassolars = mean_squared_error(y_test, y_pred_lassolars)
r2_lassolars = r2_score(y_test, y_pred_lassolars)

print(f"\nOMP (5 features):")
print(f"MSE: {mse_omp_5:.2f}, R²: {r2_omp_5:.4f}")
print(f"Features selecionadas: {np.sum(omp_5.coef_ != 0)}")

print(f"\nOMP (10 features):")
print(f"MSE: {mse_omp_10:.2f}, R²: {r2_omp_10:.4f}")
print(f"Features selecionadas: {np.sum(omp_10.coef_ != 0)}")

print(f"\nLasso:")
print(f"MSE: {mse_lasso:.2f}, R²: {r2_lasso:.4f}")
print(f"Features selecionadas: {np.sum(lasso.coef_ != 0)}")

print(f"\nLassoLars:")
print(f"MSE: {mse_lassolars:.2f}, R²: {r2_lassolars:.4f}")
print(f"Features selecionadas: {np.sum(lassolars.coef_ != 0)}")

# Implementação manual de validação cruzada para OMP
print(f"\n--- Validação Cruzada Manual para OMP ---")

n_features_range = range(1, 16)  # Limitado para evitar o erro
cv_scores = []

for n_feat in n_features_range:
    omp_cv = OrthogonalMatchingPursuit(n_nonzero_coefs=n_feat)
    scores = cross_val_score(omp_cv, X_train_scaled, y_train, 
                           cv=5, scoring='neg_mean_squared_error')
    cv_scores.append(-scores.mean())
    print(f"n_nonzero_coefs = {n_feat}: MSE CV = {-scores.mean():.2f}")

best_n_features = n_features_range[np.argmin(cv_scores)]
print(f"\nMelhor n_nonzero_coefs: {best_n_features}")

# Treinar OMP com melhor parâmetro
omp_best = OrthogonalMatchingPursuit(n_nonzero_coefs=best_n_features)
omp_best.fit(X_train_scaled, y_train)
y_pred_omp_best = omp_best.predict(X_test_scaled)
mse_omp_best = mean_squared_error(y_test, y_pred_omp_best)
r2_omp_best = r2_score(y_test, y_pred_omp_best)

print(f"\nOMP (Melhor - {best_n_features} features):")
print(f"MSE: {mse_omp_best:.2f}, R²: {r2_omp_best:.4f}")

# Visualização dos resultados
plt.figure(figsize=(15, 10))

# Gráfico 1: Comparação de coeficientes
plt.subplot(2, 3, 1)
features = range(len(true_coef))
n_plot = min(30, len(true_coef))
plt.plot(features[:n_plot], true_coef[:n_plot], 'o-', label='Verdadeiro', 
         linewidth=2, markersize=4, alpha=0.7)
plt.plot(features[:n_plot], omp_best.coef_[:n_plot], 's-', label='OMP Melhor', 
         alpha=0.8, markersize=3)
plt.plot(features[:n_plot], lasso.coef_[:n_plot], '^-', label='Lasso', 
         alpha=0.8, markersize=3)
plt.plot(features[:n_plot], lassolars.coef_[:n_plot], 'd-', label='LassoLars', 
         alpha=0.8, markersize=3)
plt.xlabel('Feature Index')
plt.ylabel('Valor do Coeficiente')
plt.title('Comparação dos Coeficientes')
plt.legend()
plt.grid(True, alpha=0.3)

# Gráfico 2: Validação Cruzada - MSE vs Número de Features
plt.subplot(2, 3, 2)
plt.plot(n_features_range, cv_scores, 'o-', linewidth=2, markersize=6)
plt.axvline(x=best_n_features, color='red', linestyle='--', 
            label=f'Melhor: {best_n_features}')
plt.xlabel('Número de Features (n_nonzero_coefs)')
plt.ylabel('MSE (Validação Cruzada)')
plt.title('Validação Cruzada - OMP')
plt.legend()
plt.grid(True, alpha=0.3)

# Gráfico 3: Features selecionadas
plt.subplot(2, 3, 3)
methods = ['OMP (5)', 'OMP (10)', 'OMP (Melhor)', 'Lasso', 'LassoLars']
n_selected = [np.sum(omp_5.coef_ != 0),
              np.sum(omp_10.coef_ != 0),
              np.sum(omp_best.coef_ != 0),
              np.sum(lasso.coef_ != 0),
              np.sum(lassolars.coef_ != 0)]

colors = ['lightcoral', 'lightcoral', 'red', 'lightgreen', 'skyblue']
plt.bar(methods, n_selected, color=colors, alpha=0.7)
plt.ylabel('Número de Features Selecionadas')
plt.title('Seleção de Features')
plt.xticks(rotation=45)
for i, v in enumerate(n_selected):
    plt.text(i, v + 0.5, str(v), ha='center', va='bottom')

# Gráfico 4: Comparação de performance
plt.subplot(2, 3, 4)
mses = [mse_omp_5, mse_omp_10, mse_omp_best, mse_lasso, mse_lassolars]
plt.bar(methods, mses, color=colors, alpha=0.7)
plt.ylabel('MSE')
plt.title('Comparação de Performance')
plt.xticks(rotation=45)
for i, v in enumerate(mses):
    plt.text(i, v + 0.5, f'{v:.1f}', ha='center', va='bottom')

# Gráfico 5: Ordem de seleção das features no OMP
plt.subplot(2, 3, 5)
# Para OMP, podemos analisar a ordem através da reconstrução
residual = y_train.copy()
selected_features = []
correlations_history = []

for step in range(omp_best.n_nonzero_coefs):
    correlations = np.abs(X_train_scaled.T @ residual)
    best_feature = np.argmax(correlations)
    selected_features.append(best_feature)
    correlations_history.append(correlations[best_feature])
    
    # Atualizar usando apenas features selecionadas
    X_active = X_train_scaled[:, selected_features]
    coef_active = np.linalg.lstsq(X_active, y_train, rcond=None)[0]
    residual = y_train - X_active @ coef_active

plt.plot(range(1, len(selected_features) + 1), selected_features, 'o-', linewidth=2)
plt.xlabel('Ordem de Seleção')
plt.ylabel('Índice da Feature')
plt.title('Ordem de Seleção das Features no OMP')
plt.grid(True, alpha=0.3)

# Gráfico 6: Correlação das features selecionadas
plt.subplot(2, 3, 6)
plt.plot(range(1, len(correlations_history) + 1), correlations_history, 's-', 
         linewidth=2, color='purple')
plt.xlabel('Ordem de Seleção')
plt.ylabel('Correlação com Resíduo')
plt.title('Correlação das Features Selecionadas')
plt.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

# Análise detalhada
print(f"\n--- Análise Detalhada ---")
true_informative = np.where(true_coef != 0)[0]
omp_selected = np.where(omp_best.coef_ != 0)[0]

print(f"Features informativas verdadeiras: {len(true_informative)}")
print(f"Features selecionadas pelo OMP: {len(omp_selected)}")
print(f"Features corretamente selecionadas: {len(np.intersect1d(true_informative, omp_selected))}")
print(f"Features informativas perdidas: {len(np.setdiff1d(true_informative, omp_selected))}")

# Métricas de seleção
def calculate_selection_metrics(true_idx, selected_idx):
    tp = len(np.intersect1d(true_idx, selected_idx))
    fp = len(np.setdiff1d(selected_idx, true_idx))
    fn = len(np.setdiff1d(true_idx, selected_idx))
    
    precision = tp / (tp + fp) if (tp + fp) > 0 else 0
    recall = tp / (tp + fn) if (tp + fn) > 0 else 0
    f1 = 2 * precision * recall / (precision + recall) if (precision + recall) > 0 else 0
    
    return precision, recall, f1

precision, recall, f1 = calculate_selection_metrics(true_informative, omp_selected)
print(f"\nMétricas de Seleção - OMP:")
print(f"Precisão: {precision:.3f}")
print(f"Recall: {recall:.3f}")
print(f"F1-Score: {f1:.3f}")

# Resumo final
print(f"\n--- Resumo Final ---")
best_model_idx = np.argmin([mse_omp_5, mse_omp_10, mse_omp_best, mse_lasso, mse_lassolars])
best_model_name = methods[best_model_idx]
print(f"Melhor modelo: {best_model_name} (MSE: {[mse_omp_5, mse_omp_10, mse_omp_best, mse_lasso, mse_lassolars][best_model_idx]:.2f})")
print(f"OMP com validação cruzada escolheu {best_n_features} features")
print(f"F1-Score do OMP: {f1:.3f}")

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

160

161

162

163

164

165

166

167

168

169

170

171

172

173

174

175

176

177

178

179

180

181

182

183

184

185

186

187

188

189

190

191

192

193

194

195

196

197

198

199

200

201

202

203

204

205

206

207

208

209

210

211

212

213

214

215

216

217

218

219

220

221

222

223

224

225

226

227

228

229

230

231

232

233

234

235

236

237

238

239

240

241

242

'''

Aplicação do Orthogonal Matching Pursuit (OMP)

para Regressão com Seleção de Features

CÓDIGO CORRIGIDO - Problema com OrthogonalMatchingPursuitCV

'''

import numpy as np

import matplotlib.pyplot as plt

from sklearn.datasets import make_regression

from sklearn.linear_model import OrthogonalMatchingPursuit

from sklearn.linear_model import Lasso, LassoLars

from sklearn.model_selection import train_test_split, cross_val_score

from sklearn.preprocessing import StandardScaler

from sklearn.metrics import mean_squared_error, r2_score

# Gerar dataset com features esparsas

n_samples = 200

n_features = 50

n_informative = 10

print("Configuração do Dataset:")

print(f"Número de amostras: {n_samples}")

print(f"Número de features: {n_features}")

print(f"Features informativas: {n_informative}")

# Gerar dados com estrutura esparsa

X, y, true_coef = make_regression(n_samples=n_samples, n_features=n_features,

n_informative=n_informative, noise=20,

coef=True, random_state=42)

print(f"\nDimensões: X {X.shape}, y {y.shape}")

# Dividir em treino e teste

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3,

random_state=42)

# Normalizar os dados

scaler = StandardScaler()

X_train_scaled = scaler.fit_transform(X_train)

X_test_scaled = scaler.transform(X_test)

print(f"\n--- Comparação: OMP vs Outros Algoritmos Esparsos ---")

# Orthogonal Matching Pursuit com diferentes números de features

omp_5 = OrthogonalMatchingPursuit(n_nonzero_coefs=5)

omp_5.fit(X_train_scaled, y_train)

y_pred_omp_5 = omp_5.predict(X_test_scaled)

mse_omp_5 = mean_squared_error(y_test, y_pred_omp_5)

r2_omp_5 = r2_score(y_test, y_pred_omp_5)

omp_10 = OrthogonalMatchingPursuit(n_nonzero_coefs=10)

omp_10.fit(X_train_scaled, y_train)

y_pred_omp_10 = omp_10.predict(X_test_scaled)

mse_omp_10 = mean_squared_error(y_test, y_pred_omp_10)

r2_omp_10 = r2_score(y_test, y_pred_omp_10)

# Lasso para comparação

lasso = Lasso(alpha=0.1, max_iter=10000, random_state=42)

lasso.fit(X_train_scaled, y_train)

y_pred_lasso = lasso.predict(X_test_scaled)

mse_lasso = mean_squared_error(y_test, y_pred_lasso)

r2_lasso = r2_score(y_test, y_pred_lasso)

# LassoLars para comparação

lassolars = LassoLars(alpha=0.1, max_iter=1000)

lassolars.fit(X_train_scaled, y_train)

y_pred_lassolars = lassolars.predict(X_test_scaled)

mse_lassolars = mean_squared_error(y_test, y_pred_lassolars)

r2_lassolars = r2_score(y_test, y_pred_lassolars)

print(f"\nOMP (5 features):")

print(f"MSE: {mse_omp_5:.2f}, R²: {r2_omp_5:.4f}")

print(f"Features selecionadas: {np.sum(omp_5.coef_ != 0)}")

print(f"\nOMP (10 features):")

print(f"MSE: {mse_omp_10:.2f}, R²: {r2_omp_10:.4f}")

print(f"Features selecionadas: {np.sum(omp_10.coef_ != 0)}")

print(f"\nLasso:")

print(f"MSE: {mse_lasso:.2f}, R²: {r2_lasso:.4f}")

print(f"Features selecionadas: {np.sum(lasso.coef_ != 0)}")

print(f"\nLassoLars:")

print(f"MSE: {mse_lassolars:.2f}, R²: {r2_lassolars:.4f}")

print(f"Features selecionadas: {np.sum(lassolars.coef_ != 0)}")

# Implementação manual de validação cruzada para OMP

print(f"\n--- Validação Cruzada Manual para OMP ---")

n_features_range = range(1, 16) # Limitado para evitar o erro

cv_scores = []

for n_feat in n_features_range:

omp_cv = OrthogonalMatchingPursuit(n_nonzero_coefs=n_feat)

scores = cross_val_score(omp_cv, X_train_scaled, y_train,

cv=5, scoring='neg_mean_squared_error')

cv_scores.append(-scores.mean())

print(f"n_nonzero_coefs = {n_feat}: MSE CV = {-scores.mean():.2f}")

best_n_features = n_features_range[np.argmin(cv_scores)]

print(f"\nMelhor n_nonzero_coefs: {best_n_features}")

# Treinar OMP com melhor parâmetro

omp_best = OrthogonalMatchingPursuit(n_nonzero_coefs=best_n_features)

omp_best.fit(X_train_scaled, y_train)

y_pred_omp_best = omp_best.predict(X_test_scaled)

mse_omp_best = mean_squared_error(y_test, y_pred_omp_best)

r2_omp_best = r2_score(y_test, y_pred_omp_best)

print(f"\nOMP (Melhor - {best_n_features} features):")

print(f"MSE: {mse_omp_best:.2f}, R²: {r2_omp_best:.4f}")

# Visualização dos resultados

plt.figure(figsize=(15, 10))

# Gráfico 1: Comparação de coeficientes

plt.subplot(2, 3, 1)

features = range(len(true_coef))

n_plot = min(30, len(true_coef))

plt.plot(features[:n_plot], true_coef[:n_plot], 'o-', label='Verdadeiro',

linewidth=2, markersize=4, alpha=0.7)

plt.plot(features[:n_plot], omp_best.coef_[:n_plot], 's-', label='OMP Melhor',

alpha=0.8, markersize=3)

plt.plot(features[:n_plot], lasso.coef_[:n_plot], '^-', label='Lasso',

alpha=0.8, markersize=3)

plt.plot(features[:n_plot], lassolars.coef_[:n_plot], 'd-', label='LassoLars',

alpha=0.8, markersize=3)

plt.xlabel('Feature Index')

plt.ylabel('Valor do Coeficiente')

plt.title('Comparação dos Coeficientes')

plt.legend()

plt.grid(True, alpha=0.3)

# Gráfico 2: Validação Cruzada - MSE vs Número de Features

plt.subplot(2, 3, 2)

plt.plot(n_features_range, cv_scores, 'o-', linewidth=2, markersize=6)

plt.axvline(x=best_n_features, color='red', linestyle='--',

label=f'Melhor: {best_n_features}')

plt.xlabel('Número de Features (n_nonzero_coefs)')

plt.ylabel('MSE (Validação Cruzada)')

plt.title('Validação Cruzada - OMP')

plt.legend()

plt.grid(True, alpha=0.3)

# Gráfico 3: Features selecionadas

plt.subplot(2, 3, 3)

methods = ['OMP (5)', 'OMP (10)', 'OMP (Melhor)', 'Lasso', 'LassoLars']

n_selected = [np.sum(omp_5.coef_ != 0),

np.sum(omp_10.coef_ != 0),

np.sum(omp_best.coef_ != 0),

np.sum(lasso.coef_ != 0),

np.sum(lassolars.coef_ != 0)]

colors = ['lightcoral', 'lightcoral', 'red', 'lightgreen', 'skyblue']

plt.bar(methods, n_selected, color=colors, alpha=0.7)

plt.ylabel('Número de Features Selecionadas')

plt.title('Seleção de Features')

plt.xticks(rotation=45)

for i, v in enumerate(n_selected):

plt.text(i, v + 0.5, str(v), ha='center', va='bottom')

# Gráfico 4: Comparação de performance

plt.subplot(2, 3, 4)

mses = [mse_omp_5, mse_omp_10, mse_omp_best, mse_lasso, mse_lassolars]

plt.bar(methods, mses, color=colors, alpha=0.7)

plt.ylabel('MSE')

plt.title('Comparação de Performance')

plt.xticks(rotation=45)

for i, v in enumerate(mses):

plt.text(i, v + 0.5, f'{v:.1f}', ha='center', va='bottom')

# Gráfico 5: Ordem de seleção das features no OMP

plt.subplot(2, 3, 5)

# Para OMP, podemos analisar a ordem através da reconstrução

residual = y_train.copy()

selected_features = []

correlations_history = []

for step in range(omp_best.n_nonzero_coefs):

correlations = np.abs(X_train_scaled.T @ residual)

best_feature = np.argmax(correlations)

selected_features.append(best_feature)

correlations_history.append(correlations[best_feature])

# Atualizar usando apenas features selecionadas

X_active = X_train_scaled[:, selected_features]

coef_active = np.linalg.lstsq(X_active, y_train, rcond=None)[0]

residual = y_train - X_active @ coef_active

plt.plot(range(1, len(selected_features) + 1), selected_features, 'o-', linewidth=2)

plt.xlabel('Ordem de Seleção')

plt.ylabel('Índice da Feature')

plt.title('Ordem de Seleção das Features no OMP')

plt.grid(True, alpha=0.3)

# Gráfico 6: Correlação das features selecionadas

plt.subplot(2, 3, 6)

plt.plot(range(1, len(correlations_history) + 1), correlations_history, 's-',

linewidth=2, color='purple')

plt.xlabel('Ordem de Seleção')

plt.ylabel('Correlação com Resíduo')

plt.title('Correlação das Features Selecionadas')

plt.grid(True, alpha=0.3)

plt.tight_layout()

plt.show()

# Análise detalhada

print(f"\n--- Análise Detalhada ---")

true_informative = np.where(true_coef != 0)[0]

omp_selected = np.where(omp_best.coef_ != 0)[0]

print(f"Features informativas verdadeiras: {len(true_informative)}")

print(f"Features selecionadas pelo OMP: {len(omp_selected)}")

print(f"Features corretamente selecionadas: {len(np.intersect1d(true_informative, omp_selected))}")

print(f"Features informativas perdidas: {len(np.setdiff1d(true_informative, omp_selected))}")

# Métricas de seleção

def calculate_selection_metrics(true_idx, selected_idx):

tp = len(np.intersect1d(true_idx, selected_idx))

fp = len(np.setdiff1d(selected_idx, true_idx))

fn = len(np.setdiff1d(true_idx, selected_idx))

precision = tp / (tp + fp) if (tp + fp) > 0 else 0

recall = tp / (tp + fn) if (tp + fn) > 0 else 0

f1 = 2 * precision * recall / (precision + recall) if (precision + recall) > 0 else 0

return precision, recall, f1

precision, recall, f1 = calculate_selection_metrics(true_informative, omp_selected)

print(f"\nMétricas de Seleção - OMP:")

print(f"Precisão: {precision:.3f}")

print(f"Recall: {recall:.3f}")

print(f"F1-Score: {f1:.3f}")

# Resumo final

print(f"\n--- Resumo Final ---")

best_model_idx = np.argmin([mse_omp_5, mse_omp_10, mse_omp_best, mse_lasso, mse_lassolars])

best_model_name = methods[best_model_idx]

print(f"Melhor modelo: {best_model_name} (MSE: {[mse_omp_5, mse_omp_10, mse_omp_best, mse_lasso, mse_lassolars][best_model_idx]:.2f})")

print(f"OMP com validação cruzada escolheu {best_n_features} features")

print(f"F1-Score do OMP: {f1:.3f}")

Vantagens do OMP

Embora existam vários algoritmos para regressão esparsa, o OMP oferece vantagens específicas:

Vantagens Principais

Controle direto: Especificação explícita do número de features
Eficiência computacional: Algoritmo guloso com complexidade controlada
Garantias teóricas: Boas propriedades de recuperação sob condições específicas
Interpretabilidade: Ordem de seleção fornece insights sobre importância

Casos de Uso Recomendados

O Orthogonal Matching Pursuit é particularmente eficaz em:

Compressed sensing: Recuperação de sinais esparsos
Seleção de features com orçamento fixo: Quando há limite no número de features
Problemas com dicionários grandes: Onde apenas poucos átomos são relevantes
Aplicações em tempo real: Onde eficiência computacional é crucial

Considerações Práticas

Algumas recomendações importantes para uso eficaz:

Use OrthogonalMatchingPursuitCV quando não souber o número ideal de features
Normalize os dados antes de aplicar OMP para melhor performance
Considere a correlação entre features, pois pode afetar a ordem de seleção
Para problemas muito grandes, verifique a escalabilidade do algoritmo

Enfim, o Orthogonal Matching Pursuit representa uma abordagem elegante e eficiente para problemas de aproximação esparsa, oferecendo controle direto sobre a esparsidade da solução e boas propriedades teóricas de recuperação.

Referência: https://scikit-learn.org/0.21/modules/linear_model.html#orthogonal-matching-pursuit-omp