"""
Explorando o kernel Exp-Seno-Quadrado para dados periódicos
Demonstra como capturar padrões cíclicos com variações suaves
"""
import numpy as np
import matplotlib.pyplot as plt
from sklearn.gaussian_process.kernels import ExpSineSquared, RBF
from sklearn.gaussian_process import GaussianProcessRegressor
# Criando dados com diferentes tipos de periodicidade
def criar_dados_periodicos():
"""Cria dados com vários tipos de padrões periódicos do mundo real"""
np.random.seed(42)
# Dados com periodicidade perfeita (para comparação)
t_perfeito = np.linspace(0, 20, 200).reshape(-1, 1)
y_perfeito = np.sin(2 * np.pi * t_perfeito.ravel() / 5) + 0.1 * np.random.normal(size=200)
# Dados com periodicidade variável (mais realista)
t_variavel = np.linspace(0, 20, 200).reshape(-1, 1)
periodo_variavel = 5 + 0.5 * np.sin(t_variavel.ravel() / 3) # Período que varia suavemente
y_variavel = np.sin(2 * np.pi * t_variavel.ravel() / periodo_variavel) + 0.1 * np.random.normal(size=200)
# Dados com múltiplas periodicidades (caso complexo)
t_multiplo = np.linspace(0, 30, 300).reshape(-1, 1)
y_multiplo = (np.sin(2 * np.pi * t_multiplo.ravel() / 3) + # Período curto
0.5 * np.sin(2 * np.pi * t_multiplo.ravel() / 8) + # Período médio
0.3 * np.sin(2 * np.pi * t_multiplo.ravel() / 15) + # Período longo
0.1 * np.random.normal(size=300))
# Dados do mundo real simulados: batimentos cardíacos
t_coracao = np.linspace(0, 10, 500).reshape(-1, 1)
# Simula batimentos com pequenas variações no ritmo
frequencia_base = 1.2 # Hz (72 batimentos por minuto)
variacao_ritmo = 0.1 * np.sin(2 * np.pi * t_coracao.ravel() / 4) # Variação respiratória
y_coracao = np.sin(2 * np.pi * (frequencia_base + variacao_ritmo) * t_coracao.ravel())
y_coracao += 0.05 * np.random.normal(size=500) # Ruído de medição
return (t_perfeito, y_perfeito), (t_variavel, y_variavel), (t_multiplo, y_multiplo), (t_coracao, y_coracao)
dados_perfeito, dados_variavel, dados_multiplo, dados_coracao = criar_dados_periodicos()
print("=== EXPLORANDO O KERNEL EXP-SENO-QUADRADO ===\n")
print("O kernel Exp-Seno-Quadrado é definido por:")
print("[latex]k(x_i, x_j) = \\exp\\left(-\\frac{2\\sin^2(\\pi d(x_i, x_j)/p)}{l^2}\\right)[/latex]")
print("Onde p é o período e l é o length_scale que controla a suavidade\n")
# Explorando diferentes parâmetros do kernel
parametros_testar = [
{'periodicity': 5.0, 'length_scale': 1.0},
{'periodicity': 5.0, 'length_scale': 0.5},
{'periodicity': 5.0, 'length_scale': 2.0},
{'periodicity': 3.0, 'length_scale': 1.0},
{'periodicity': 8.0, 'length_scale': 1.0}
]
fig, axes = plt.subplots(2, 3, figsize=(18, 12))
axes = axes.ravel()
t_dados, y_dados = dados_perfeito
for idx, params in enumerate(parametros_testar):
kernel_ess = ExpSineSquared(length_scale=params['length_scale'],
periodicity=params['periodicity'])
gpr = GaussianProcessRegressor(kernel=kernel_ess, alpha=0.0, random_state=42)
gpr.fit(t_dados, y_dados)
# Previsões
t_pred = np.linspace(0, 20, 300).reshape(-1, 1)
y_pred, sigma = gpr.predict(t_pred, return_std=True)
# Plot
axes[idx].plot(t_dados, y_dados, 'ro', alpha=0.4, markersize=2, label='Dados')
axes[idx].plot(t_pred, y_pred, 'b-', linewidth=2, label='Previsão ESS')
axes[idx].fill_between(t_pred.ravel(), y_pred - 1.96*sigma,
y_pred + 1.96*sigma, alpha=0.2, label='Incerteza 95%')
axes[idx].set_title(f'ESS: p={params["periodicity"]}, l={params["length_scale"]}\n'
f'LL: {gpr.log_marginal_likelihood():.2f}', fontsize=11)
axes[idx].grid(True, alpha=0.3)
axes[idx].legend(fontsize=9)
# Comparação com RBF no último subplot
kernel_rbf = RBF(length_scale=1.0)
gpr_rbf = GaussianProcessRegressor(kernel=kernel_rbf, alpha=0.0, random_state=42)
gpr_rbf.fit(t_dados, y_dados)
y_pred_rbf, sigma_rbf = gpr_rbf.predict(t_pred, return_std=True)
axes[5].plot(t_dados, y_dados, 'ro', alpha=0.4, markersize=2, label='Dados')
axes[5].plot(t_pred, y_pred_rbf, 'g-', linewidth=2, label='Previsão RBF')
axes[5].fill_between(t_pred.ravel(), y_pred_rbf - 1.96*sigma_rbf,
y_pred_rbf + 1.96*sigma_rbf, alpha=0.2, label='Incerteza 95%')
axes[5].set_title(f'RBF (comparação)\nLL: {gpr_rbf.log_marginal_likelihood():.2f}',
fontsize=11)
axes[5].grid(True, alpha=0.3)
axes[5].legend(fontsize=9)
plt.tight_layout()
plt.show()
# Testando em diferentes cenários periódicos
print("\n=== TESTANDO EM DIFERENTES CENÁRIOS PERIÓDICOS ===\n")
cenarios = {
'Periodicidade Perfeita': dados_perfeito,
'Periodicidade Variável': dados_variavel,
'Múltiplas Periodicidades': dados_multiplo,
'Batimentos Cardíacos': dados_coracao
}
resultados_comparacao = {}
fig, axes = plt.subplots(len(cenarios), 2, figsize=(15, 12))
for i, (nome_cenario, (t_data, y_data)) in enumerate(cenarios.items()):
# Kernel Exp-Seno-Quadrado com otimização
kernel_ess = ExpSineSquared(length_scale=1.0, periodicity=5.0)
gpr_ess = GaussianProcessRegressor(kernel=kernel_ess, random_state=42)
gpr_ess.fit(t_data, y_data)
# RBF para comparação
kernel_rbf = RBF(length_scale=1.0)
gpr_rbf = GaussianProcessRegressor(kernel=kernel_rbf, random_state=42)
gpr_rbf.fit(t_data, y_data)
# Previsões
t_pred = np.linspace(t_data.min(), t_data.max(), 300).reshape(-1, 1)
y_pred_ess, sigma_ess = gpr_ess.predict(t_pred, return_std=True)
y_pred_rbf, sigma_rbf = gpr_rbf.predict(t_pred, return_std=True)
# Plot Exp-Seno-Quadrado
axes[i, 0].plot(t_data, y_data, 'ro', alpha=0.5, markersize=2, label='Dados')
axes[i, 0].plot(t_pred, y_pred_ess, 'b-', linewidth=2, label='Previsão ESS')
axes[i, 0].fill_between(t_pred.ravel(), y_pred_ess - 1.96*sigma_ess,
y_pred_ess + 1.96*sigma_ess, alpha=0.2, label='Incerteza')
axes[i, 0].set_title(f'Exp-Seno-Quadrado - {nome_cenario}\n'
f'Kernel: {gpr_ess.kernel_}\n'
f'LL: {gpr_ess.log_marginal_likelihood():.2f}', fontsize=10)
axes[i, 0].grid(True, alpha=0.3)
axes[i, 0].legend(fontsize=8)
# Plot RBF
axes[i, 1].plot(t_data, y_data, 'ro', alpha=0.5, markersize=2, label='Dados')
axes[i, 1].plot(t_pred, y_pred_rbf, 'g-', linewidth=2, label='Previsão RBF')
axes[i, 1].fill_between(t_pred.ravel(), y_pred_rbf - 1.96*sigma_rbf,
y_pred_rbf + 1.96*sigma_rbf, alpha=0.2, label='Incerteza')
axes[i, 1].set_title(f'RBF - {nome_cenario}\n'
f'Kernel: {gpr_rbf.kernel_}\n'
f'LL: {gpr_rbf.log_marginal_likelihood():.2f}', fontsize=10)
axes[i, 1].grid(True, alpha=0.3)
axes[i, 1].legend(fontsize=8)
# Armazenar resultados
resultados_comparacao[nome_cenario] = {
'ESS_LL': gpr_ess.log_marginal_likelihood(),
'RBF_LL': gpr_rbf.log_marginal_likelihood(),
'ESS_kernel': str(gpr_ess.kernel_),
'RBF_kernel': str(gpr_rbf.kernel_)
}
plt.tight_layout()
plt.show()
# Análise dos resultados
print("=== ANÁLISE DOS RESULTADOS (Log-Likelihood) ===\n")
for cenario, resultados in resultados_comparacao.items():
ess_ll = resultados['ESS_LL']
rbf_ll = resultados['RBF_LL']
diferenca = ess_ll - rbf_ll
melhor = "Exp-Seno-Quadrado" if diferenca > 0 else "RBF"
print(f"{cenario}:")
print(f" Exp-Seno-Quadrado: {ess_ll:8.2f}")
print(f" RBF: {rbf_ll:8.2f}")
print(f" Diferença: {diferenca:8.2f} → {melhor} é melhor")
print(f" Kernel ESS otimizado: {resultados['ESS_kernel']}")
print()
# Visualização das funções de covariância
print("=== FUNÇÕES DE COVARIÂNCIA EXP-SENO-QUADRADO ===\n")
plt.figure(figsize=(15, 10))
# Plot das funções de covariância para diferentes períodos
plt.subplot(2, 2, 1)
distances = np.linspace(0, 10, 100)
periodos = [3.0, 5.0, 8.0, 12.0]
for periodo in periodos:
kernel = ExpSineSquared(length_scale=1.0, periodicity=periodo)
K = kernel(np.array([[0]]), distances.reshape(-1, 1))
plt.plot(distances, K, label=f'p = {periodo}', linewidth=2)
plt.title('Funções de Covariância - Diferentes Períodos')
plt.xlabel('Distância')
plt.ylabel('Covariância')
plt.legend()
plt.grid(True, alpha=0.3)
# Plot para diferentes length_scales
plt.subplot(2, 2, 2)
length_scales = [0.5, 1.0, 2.0, 5.0]
for ls in length_scales:
kernel = ExpSineSquared(length_scale=ls, periodicity=5.0)
K = kernel(np.array([[0]]), distances.reshape(-1, 1))
plt.plot(distances, K, label=f'l = {ls}', linewidth=2)
plt.title('Funções de Covariância - Diferentes Length Scales')
plt.xlabel('Distância')
plt.ylabel('Covariância')
plt.legend()
plt.grid(True, alpha=0.3)
# Comparação com RBF
plt.subplot(2, 2, 3)
kernel_ess = ExpSineSquared(length_scale=1.0, periodicity=5.0)
K_ess = kernel_ess(np.array([[0]]), distances.reshape(-1, 1))
kernel_rbf = RBF(length_scale=1.0)
K_rbf = kernel_rbf(np.array([[0]]), distances.reshape(-1, 1))
plt.plot(distances, K_ess, 'b-', linewidth=2, label='Exp-Seno-Quadrado (p=5)')
plt.plot(distances, K_rbf, 'g--', linewidth=2, label='RBF')
plt.title('Comparação: Exp-Seno-Quadrado vs RBF')
plt.xlabel('Distância')
plt.ylabel('Covariância')
plt.legend()
plt.grid(True, alpha=0.3)
# Aplicações do mundo real
plt.subplot(2, 2, 4)
aplicacoes = [
('Clima/Tempo', 'Sazonalidade anual/diária', 'p = 24 (horas), 365 (dias)'),
('Finanças', 'Ciclos econômicos/sazonais', 'p = 12 (meses), 4 (trimestres)'),
('Medicina', 'Batimentos cardíacos/respiração', 'p = 0.8-1.2 (segundos)'),
('Astronomia', 'Movimentos planetários/fases', 'p = 27.3 (lunar), 365 (solar)'),
('Engenharia', 'Vibrações/rotações mecânicas', 'p = 1/frequência')
]
y_pos = range(len(aplicacoes))
app_names = [app[0] for app in aplicacoes]
periodos_typical = [365, 12, 1, 27.3, 0.1] # Períodos típicos
plt.barh(y_pos, periodos_typical, color='lightsteelblue', alpha=0.7)
plt.yticks(y_pos, app_names)
plt.xlabel('Período Típico')
plt.title('Aplicações do Mundo Real')
for i, (nome, desc, periodo) in enumerate(aplicacoes):
plt.text(10, i, f'{desc}\n{periodo}', va='center', ha='left', fontsize=8)
plt.tight_layout()
plt.show()
# Detectando periodicidade automaticamente
print("=== DETECTANDO PERIODICIDADE AUTOMATICAMENTE ===\n")
def detectar_periodicidade(t, y):
"""Ferramenta simples para detectar periodicidade nos dados"""
from scipy import signal
# Usando FFT para encontrar frequências dominantes
y_detrended = y - np.mean(y)
fft = np.fft.fft(y_detrended)
freqs = np.fft.fftfreq(len(y), t[1]-t[0])
# Encontrando pico não-DC
magnitude = np.abs(fft)
magnitude[0] = 0 # Ignorar componente DC
idx_pico = np.argmax(magnitude[1:]) + 1
frequencia_dominante = np.abs(freqs[idx_pico])
periodo_estimado = 1.0 / frequencia_dominante if frequencia_dominante > 0 else 0
# Análise de autocorrelação
autocorr = signal.correlate(y_detrended, y_detrended, mode='full')
autocorr = autocorr[len(autocorr)//2:]
lag_pico = signal.find_peaks(autocorr, distance=10)[0]
if len(lag_pico) > 1:
periodo_autocorr = (lag_pico[1] - lag_pico[0]) * (t[1]-t[0])
else:
periodo_autocorr = 0
return periodo_estimado, periodo_autocorr
# Testando detecção em nossos dados
for nome, (t_data, y_data) in cenarios.items():
periodo_fft, periodo_autocorr = detectar_periodicidade(t_data.ravel(), y_data)
print(f"{nome}:")
print(f" Período estimado (FFT): {periodo_fft:.2f}")
print(f" Período estimado (Autocorr): {periodo_autocorr:.2f}")
print()
# Resumo prático
print("=== RESUMO PRÁTICO: QUANDO USAR EXP-SENO-QUADRADO ===\n")
vantagens_ess = [
"• Captura especificamente padrões periódicos e quase-periódicos",
"• Permite variações suaves na fase e amplitude dos ciclos",
"• Parâmetro de período é interpretável diretamente",
"• Excelente para dados sazonais, climáticos e biológicos",
"• Pode ser combinado com outros kernels para padrões complexos"
]
casos_uso_especificos = [
"Dados climáticos e meteorológicos (sazonalidade)",
"Sinais fisiológicos (batimentos cardíacos, respiração)",
"Dados econômicos e financeiros (ciclos sazonais)",
"Sinais astronômicos (movimentos planetários)",
"Vibrações mecânicas e sinais de engenharia"
]
print("VANTAGENS PRINCIPAIS:")
for vantagem in vantagens_ess:
print(vantagem)
print("\nCASOS DE USO RECOMENDADOS:")
for caso in casos_uso_especificos:
print(f" • {caso}")
print("\nCONFIGURAÇÃO INICIAL RECOMENDADA:")
print(" • Use detecção automática para estimar o período inicial")
print(" • Comece com length_scale = 1.0 e ajuste conforme necessário")
print(" • Combine com RBF para capturar tendências não-periódicas")
print(" • Use bounds realistas baseados no conhecimento do domínio")
