Gaussian Process Modeling for IOPS Web Server KPI¶
Objective: Build a robust Gaussian Process model for operational time series forecasting and anomaly detection using GPyTorch with sparse variational approximations.
Key Innovation: Student-t likelihood + composite periodic kernels trained on ALL data (including anomalies) for production-ready robustness.
0. Auto-Reload Configuration¶
Hot Reload: Enable automatic reloading of library code (src/) without kernel restart.
In [ ]:
Copied!
# Auto-reload: Picks up library changes without kernel restart
%load_ext autoreload
%autoreload 2
# Auto-reload: Picks up library changes without kernel restart
%load_ext autoreload
%autoreload 2
1. Environment Setup¶
In [ ]:
Copied!
# Environment Setup
# Local: Uses installed hellocloud
# Colab: Installs from GitHub
try:
import hellocloud
except ImportError:
!pip install -q git+https://github.com/nehalecky/hello-cloud.git
import hellocloud
# Environment Setup
# Local: Uses installed hellocloud
# Colab: Installs from GitHub
try:
import hellocloud
except ImportError:
!pip install -q git+https://github.com/nehalecky/hello-cloud.git
import hellocloud
In [ ]:
Copied!
# Core imports
import numpy as np
# Polars replaced with PySpark
# PyTorch and GPyTorch
import torch
import gpytorch
from gpytorch.likelihoods import StudentTLikelihood, GaussianLikelihood
# Cloud simulation library (our GP implementation)
from hellocloud.modeling.gaussian_process import (
CompositePeriodicKernel,
SparseGPModel,
initialize_inducing_points,
train_gp_model,
save_model,
load_model,
compute_metrics,
compute_anomaly_metrics,
compute_prediction_intervals,
)
# Visualization
import matplotlib.pyplot as plt
import seaborn as sns
# Metrics
from sklearn.metrics import (
precision_score, recall_score, f1_score, roc_auc_score,
roc_curve
)
from scipy.stats import norm, t as student_t
# Configuration
import warnings
warnings.filterwarnings('ignore')
sns.set_style("whitegrid")
sns.set_context("notebook", font_scale=1.1)
# Set random seeds
torch.manual_seed(42)
np.random.seed(42)
# Device selection
# NOTE: MPS (Apple Silicon) is not used - GPyTorch's variational inference requires
# operations (_linalg_eigh) not yet implemented in PyTorch's MPS backend.
# Using CPU is reliable and still performs well for this dataset size.
if torch.cuda.is_available():
device = torch.device('cuda')
print(f"✓ Using CUDA GPU: {torch.cuda.get_device_name(0)}")
else:
device = torch.device('cpu')
print(f"✓ Using CPU")
# Use float32 for all operations (standard for deep learning)
dtype = torch.float32
# Core imports
import numpy as np
# Polars replaced with PySpark
# PyTorch and GPyTorch
import torch
import gpytorch
from gpytorch.likelihoods import StudentTLikelihood, GaussianLikelihood
# Cloud simulation library (our GP implementation)
from hellocloud.modeling.gaussian_process import (
CompositePeriodicKernel,
SparseGPModel,
initialize_inducing_points,
train_gp_model,
save_model,
load_model,
compute_metrics,
compute_anomaly_metrics,
compute_prediction_intervals,
)
# Visualization
import matplotlib.pyplot as plt
import seaborn as sns
# Metrics
from sklearn.metrics import (
precision_score, recall_score, f1_score, roc_auc_score,
roc_curve
)
from scipy.stats import norm, t as student_t
# Configuration
import warnings
warnings.filterwarnings('ignore')
sns.set_style("whitegrid")
sns.set_context("notebook", font_scale=1.1)
# Set random seeds
torch.manual_seed(42)
np.random.seed(42)
# Device selection
# NOTE: MPS (Apple Silicon) is not used - GPyTorch's variational inference requires
# operations (_linalg_eigh) not yet implemented in PyTorch's MPS backend.
# Using CPU is reliable and still performs well for this dataset size.
if torch.cuda.is_available():
device = torch.device('cuda')
print(f"✓ Using CUDA GPU: {torch.cuda.get_device_name(0)}")
else:
device = torch.device('cpu')
print(f"✓ Using CPU")
# Use float32 for all operations (standard for deep learning)
dtype = torch.float32
2. EDA Recap: Key Findings¶
From our exploratory data analysis (03_iops_web_server_eda.ipynb), we discovered:
Dataset Characteristics¶
- Training: 146,255 samples with 285 labeled anomalies (0.19%)
- Testing: 149,130 samples with 991 labeled anomalies (0.66%)
- KPI: IOPS (I/O operations per second) from production web server
Two-Scale Periodic Pattern¶
SLOW Component (Sawtooth Envelope):
- Period: ~1250 timesteps (≈21 hours)
- Shape: Linear ramp-up → sharp drop/reset
- Operational interpretation: Daily accumulation + overnight reset
FAST Component (Sinusoidal Carrier):
- Period: ~250 timesteps (≈4 hours)
- 5 complete oscillations per sawtooth cycle (valley to valley)
- Operational interpretation: Regular micro-cycles within operational windows
Statistical Properties¶
- Normal periods: Mean ≈ 34.4, Std ≈ 3.9
- Anomalous periods: Mean ≈ 39.4, Std ≈ 15.8 (4× larger variance)
- KS test: Highly significant distributional differences (D=0.417, p<0.001)
Modeling Implications¶
- Robust approach needed: Student-t likelihood handles heavy tails
- Sparse GP required: Dataset size (n=146,255) requires inducing points
- Composite kernel: Capture sawtooth × sinusoidal interaction
3. Data Loading and Preprocessing¶
In [ ]:
Copied!
# Load IOPS data from HuggingFace (same as EDA notebook)
import pandas as pd
base_url = "https://huggingface.co/datasets/AutonLab/Timeseries-PILE/resolve/main"
kpi_id = "KPI-05f10d3a-239c-3bef-9bdc-a2feeb0037aa"
# Load train and test splits
train_url = f"{base_url}/anomaly_detection/TSB-UAD-Public/IOPS/{kpi_id}.train.out"
test_url = f"{base_url}/anomaly_detection/TSB-UAD-Public/IOPS/{kpi_id}.test.out"
print("Downloading IOPS data from HuggingFace...")
train_pd = pd.read_csv(train_url, header=None, names=['value', 'label'])
test_pd = pd.read_csv(test_url, header=None, names=['value', 'label'])
# Convert to Polars
# NOTE: Dataset has no timestamps - we create sequential indices (1-minute intervals)
train_df = spark.createDataFrame({
'timestamp': np.arange(len(train_pd)),
'value': train_pd['value'].values,
'is_anomaly': train_pd['label'].values.astype(bool)
})
test_df = spark.createDataFrame({
'timestamp': np.arange(len(test_pd)),
'value': test_pd['value'].values,
'is_anomaly': test_pd['label'].values.astype(bool)
})
print(f"✓ Data loaded successfully")
print(f" Training: {len(train_df):,} samples")
print(f" Test: {len(test_df):,} samples")
# Load IOPS data from HuggingFace (same as EDA notebook)
import pandas as pd
base_url = "https://huggingface.co/datasets/AutonLab/Timeseries-PILE/resolve/main"
kpi_id = "KPI-05f10d3a-239c-3bef-9bdc-a2feeb0037aa"
# Load train and test splits
train_url = f"{base_url}/anomaly_detection/TSB-UAD-Public/IOPS/{kpi_id}.train.out"
test_url = f"{base_url}/anomaly_detection/TSB-UAD-Public/IOPS/{kpi_id}.test.out"
print("Downloading IOPS data from HuggingFace...")
train_pd = pd.read_csv(train_url, header=None, names=['value', 'label'])
test_pd = pd.read_csv(test_url, header=None, names=['value', 'label'])
# Convert to Polars
# NOTE: Dataset has no timestamps - we create sequential indices (1-minute intervals)
train_df = spark.createDataFrame({
'timestamp': np.arange(len(train_pd)),
'value': train_pd['value'].values,
'is_anomaly': train_pd['label'].values.astype(bool)
})
test_df = spark.createDataFrame({
'timestamp': np.arange(len(test_pd)),
'value': test_pd['value'].values,
'is_anomaly': test_pd['label'].values.astype(bool)
})
print(f"✓ Data loaded successfully")
print(f" Training: {len(train_df):,} samples")
print(f" Test: {len(test_df):,} samples")
In [ ]:
Copied!
# Extract arrays for modeling
X_train = train_df['timestamp'].to_numpy().reshape(-1, 1).astype(np.float64)
y_train = train_df['value'].to_numpy().astype(np.float64)
anomaly_train = train_df['is_anomaly'].to_numpy()
n_train = len(X_train)
n_anomalies_train = anomaly_train.sum()
print(f"Training data:")
print(f" Total samples: {n_train:,}")
print(f" Anomalies: {n_anomalies_train} ({100*n_anomalies_train/n_train:.2f}%)")
print(f" Time range: {X_train.min():.0f} → {X_train.max():.0f}")
# Extract arrays for modeling
X_train = train_df['timestamp'].to_numpy().reshape(-1, 1).astype(np.float64)
y_train = train_df['value'].to_numpy().astype(np.float64)
anomaly_train = train_df['is_anomaly'].to_numpy()
n_train = len(X_train)
n_anomalies_train = anomaly_train.sum()
print(f"Training data:")
print(f" Total samples: {n_train:,}")
print(f" Anomalies: {n_anomalies_train} ({100*n_anomalies_train/n_train:.2f}%)")
print(f" Time range: {X_train.min():.0f} → {X_train.max():.0f}")
In [ ]:
Copied!
# Extract test arrays
X_test = test_df['timestamp'].to_numpy().reshape(-1, 1).astype(np.float64)
y_test = test_df['value'].to_numpy().astype(np.float64)
anomaly_test = test_df['is_anomaly'].to_numpy()
n_test = len(X_test)
n_anomalies_test = anomaly_test.sum()
print(f"Test data:")
print(f" Total samples: {n_test:,}")
print(f" Anomalies: {n_anomalies_test} ({100*n_anomalies_test/n_test:.2f}%)")
print(f" Time range: {X_test.min():.0f} → {X_test.max():.0f}")
# Extract test arrays
X_test = test_df['timestamp'].to_numpy().reshape(-1, 1).astype(np.float64)
y_test = test_df['value'].to_numpy().astype(np.float64)
anomaly_test = test_df['is_anomaly'].to_numpy()
n_test = len(X_test)
n_anomalies_test = anomaly_test.sum()
print(f"Test data:")
print(f" Total samples: {n_test:,}")
print(f" Anomalies: {n_anomalies_test} ({100*n_anomalies_test/n_test:.2f}%)")
print(f" Time range: {X_test.min():.0f} → {X_test.max():.0f}")
In [ ]:
Copied!
# Normalize timestamps for numerical stability
X_min = X_train.min()
X_max = X_train.max()
X_range = X_max - X_min
X_train_norm = (X_train - X_min) / X_range
X_test_norm = (X_test - X_min) / X_range
print(f"Normalized timestamps:")
print(f" Training: {X_train_norm.min():.6f} → {X_train_norm.max():.6f}")
print(f" Test: {X_test_norm.min():.6f} → {X_test_norm.max():.6f}")
# Normalize timestamps for numerical stability
X_min = X_train.min()
X_max = X_train.max()
X_range = X_max - X_min
X_train_norm = (X_train - X_min) / X_range
X_test_norm = (X_test - X_min) / X_range
print(f"Normalized timestamps:")
print(f" Training: {X_train_norm.min():.6f} → {X_train_norm.max():.6f}")
print(f" Test: {X_test_norm.min():.6f} → {X_test_norm.max():.6f}")
In [ ]:
Copied!
# Create clean training set (exclude anomalies) for baseline model
mask_clean = ~anomaly_train.astype(bool)
X_train_clean = X_train_norm[mask_clean]
y_train_clean = y_train[mask_clean]
print(f"Clean training set (traditional approach):")
print(f" Samples: {len(X_train_clean):,} (excluded {n_anomalies_train} anomalies)")
# Create clean training set (exclude anomalies) for baseline model
mask_clean = ~anomaly_train.astype(bool)
X_train_clean = X_train_norm[mask_clean]
y_train_clean = y_train[mask_clean]
print(f"Clean training set (traditional approach):")
print(f" Samples: {len(X_train_clean):,} (excluded {n_anomalies_train} anomalies)")
In [ ]:
Copied!
# Convert to PyTorch tensors (use dtype from device selection)
X_train_t = torch.tensor(X_train_norm, dtype=dtype, device=device)
y_train_t = torch.tensor(y_train, dtype=dtype, device=device)
X_train_clean_t = torch.tensor(X_train_clean, dtype=dtype, device=device)
y_train_clean_t = torch.tensor(y_train_clean, dtype=dtype, device=device)
X_test_t = torch.tensor(X_test_norm, dtype=dtype, device=device)
y_test_t = torch.tensor(y_test, dtype=dtype, device=device)
print(f"Tensor shapes:")
print(f" X_train: {X_train_t.shape}, y_train: {y_train_t.shape}")
print(f" X_train_clean: {X_train_clean_t.shape}, y_train_clean: {y_train_clean_t.shape}")
print(f" X_test: {X_test_t.shape}, y_test: {y_test_t.shape}")
# Convert to PyTorch tensors (use dtype from device selection)
X_train_t = torch.tensor(X_train_norm, dtype=dtype, device=device)
y_train_t = torch.tensor(y_train, dtype=dtype, device=device)
X_train_clean_t = torch.tensor(X_train_clean, dtype=dtype, device=device)
y_train_clean_t = torch.tensor(y_train_clean, dtype=dtype, device=device)
X_test_t = torch.tensor(X_test_norm, dtype=dtype, device=device)
y_test_t = torch.tensor(y_test, dtype=dtype, device=device)
print(f"Tensor shapes:")
print(f" X_train: {X_train_t.shape}, y_train: {y_train_t.shape}")
print(f" X_train_clean: {X_train_clean_t.shape}, y_train_clean: {y_train_clean_t.shape}")
print(f" X_test: {X_test_t.shape}, y_test: {y_test_t.shape}")
In [ ]:
Copied!
# Numerical stability configuration for Cholesky decomposition
cholesky_jitter = 1e-3 # Diagonal regularization
cholesky_max_tries = 10 # Retry attempts if decomposition fails
print(f"✓ Numerical stability: jitter={cholesky_jitter}, max_tries={cholesky_max_tries}")
# Numerical stability configuration for Cholesky decomposition
cholesky_jitter = 1e-3 # Diagonal regularization
cholesky_max_tries = 10 # Retry attempts if decomposition fails
print(f"✓ Numerical stability: jitter={cholesky_jitter}, max_tries={cholesky_max_tries}")
4. Exact GP Baseline (Simple Approach)¶
Why Exact GP?¶
The sparse variational GP approach (sections 5-9) uses inducing points for scalability:
- Dataset: 146k timesteps
- Inducing points: M=200
- Spacing: ~730 timesteps apart
- Fast period: 250 timesteps
Problem: With inducing points 3× the period spacing, the variational approximation cannot capture fine periodic structure. The model smooths over the patterns.
Solution: Use exact GP on a subset of data (following GPyTorch simple example):
- No inducing points needed
- Can capture true periodic structure
- Computational cost O(n³) limits to n≈5-10k points
4.1 Subsample Training Data¶
In [ ]:
Copied!
# Subsample to make exact GP tractable
# Take every Nth point to get ~5000 samples
subsample_factor = 30 # 146,255 / 30 ≈ 4,875 points
# Create subsampled indices (evenly spaced)
subsample_indices = np.arange(0, len(X_train), subsample_factor)
X_train_sub = X_train[subsample_indices]
y_train_sub = y_train[subsample_indices]
anomaly_train_sub = anomaly_train[subsample_indices]
print(f"Subsampled training data:")
print(f" Original: {len(X_train):,} samples")
print(f" Subsampled: {len(X_train_sub):,} samples (every {subsample_factor}th point)")
print(f" Reduction: {100*(1-len(X_train_sub)/len(X_train)):.1f}%")
print(f" Anomalies: {anomaly_train_sub.sum()} ({100*anomaly_train_sub.sum()/len(X_train_sub):.2f}%)")
# Subsample to make exact GP tractable
# Take every Nth point to get ~5000 samples
subsample_factor = 30 # 146,255 / 30 ≈ 4,875 points
# Create subsampled indices (evenly spaced)
subsample_indices = np.arange(0, len(X_train), subsample_factor)
X_train_sub = X_train[subsample_indices]
y_train_sub = y_train[subsample_indices]
anomaly_train_sub = anomaly_train[subsample_indices]
print(f"Subsampled training data:")
print(f" Original: {len(X_train):,} samples")
print(f" Subsampled: {len(X_train_sub):,} samples (every {subsample_factor}th point)")
print(f" Reduction: {100*(1-len(X_train_sub)/len(X_train)):.1f}%")
print(f" Anomalies: {anomaly_train_sub.sum()} ({100*anomaly_train_sub.sum()/len(X_train_sub):.2f}%)")
4.1.1 Subsampling Validation: Does It Preserve Signal Structure?¶
In [ ]:
Copied!
# CRITICAL: Verify subsampling doesn't destroy the signal
# Save all outputs to file for easy reference
output_file = '../subsampling_validation.txt'
with open(output_file, 'w') as f:
f.write("=" * 80 + "\n")
f.write("SUBSAMPLING VALIDATION REPORT\n")
f.write("=" * 80 + "\n\n")
# Statistical comparison
f.write("Statistical Comparison:\n")
f.write("-" * 80 + "\n")
f.write(f"{'Metric':<20} {'Full Data':<15} {'Subsampled':<15} {'Difference':<15}\n")
f.write("-" * 80 + "\n")
metrics = {
'Mean': (y_train.mean(), y_train_sub.mean()),
'Std': (y_train.std(), y_train_sub.std()),
'Variance': (y_train.var(), y_train_sub.var()),
'Min': (y_train.min(), y_train_sub.min()),
'Max': (y_train.max(), y_train_sub.max()),
'Q25': (np.percentile(y_train, 25), np.percentile(y_train_sub, 25)),
'Median': (np.median(y_train), np.median(y_train_sub)),
'Q75': (np.percentile(y_train, 75), np.percentile(y_train_sub, 75)),
}
for metric, (full, sub) in metrics.items():
diff = ((sub - full) / full) * 100 if full != 0 else 0
line = f"{metric:<20} {full:<15.3f} {sub:<15.3f} {diff:>+6.1f}%\n"
f.write(line)
print(line, end='')
f.write("=" * 80 + "\n\n")
print()
# Autocorrelation comparison (critical for temporal structure)
from scipy.stats import pearsonr
# Compute autocorrelation at key lags
lags_to_check = [1, 10, 50, 250, 1250] # Include our identified periods
f.write("Autocorrelation Comparison:\n")
f.write("-" * 80 + "\n")
f.write(f"{'Lag':<10} {'Full Data':<15} {'Subsampled':<15} {'Difference':<15}\n")
f.write("-" * 80 + "\n")
print("Autocorrelation Comparison:")
print("-" * 80)
for lag in lags_to_check:
# Full data autocorrelation
if lag < len(y_train):
acf_full = pearsonr(y_train[:-lag], y_train[lag:])[0]
else:
acf_full = np.nan
# Subsampled autocorrelation (adjust lag for subsampling)
lag_sub = lag // subsample_factor
if lag_sub > 0 and lag_sub < len(y_train_sub):
acf_sub = pearsonr(y_train_sub[:-lag_sub], y_train_sub[lag_sub:])[0]
else:
acf_sub = np.nan
if not np.isnan(acf_full) and not np.isnan(acf_sub):
diff = acf_sub - acf_full
line = f"{lag:<10} {acf_full:<15.3f} {acf_sub:<15.3f} {diff:>+6.3f}\n"
f.write(line)
print(line, end='')
else:
line = f"{lag:<10} {str(acf_full):<15} {str(acf_sub):<15} {'N/A':<15}\n"
f.write(line)
print(line, end='')
f.write("=" * 80 + "\n")
print("=" * 80)
print(f"\n✓ Validation report saved to: {output_file}")
print(f" Use: cat {output_file}")
# CRITICAL: Verify subsampling doesn't destroy the signal
# Save all outputs to file for easy reference
output_file = '../subsampling_validation.txt'
with open(output_file, 'w') as f:
f.write("=" * 80 + "\n")
f.write("SUBSAMPLING VALIDATION REPORT\n")
f.write("=" * 80 + "\n\n")
# Statistical comparison
f.write("Statistical Comparison:\n")
f.write("-" * 80 + "\n")
f.write(f"{'Metric':<20} {'Full Data':<15} {'Subsampled':<15} {'Difference':<15}\n")
f.write("-" * 80 + "\n")
metrics = {
'Mean': (y_train.mean(), y_train_sub.mean()),
'Std': (y_train.std(), y_train_sub.std()),
'Variance': (y_train.var(), y_train_sub.var()),
'Min': (y_train.min(), y_train_sub.min()),
'Max': (y_train.max(), y_train_sub.max()),
'Q25': (np.percentile(y_train, 25), np.percentile(y_train_sub, 25)),
'Median': (np.median(y_train), np.median(y_train_sub)),
'Q75': (np.percentile(y_train, 75), np.percentile(y_train_sub, 75)),
}
for metric, (full, sub) in metrics.items():
diff = ((sub - full) / full) * 100 if full != 0 else 0
line = f"{metric:<20} {full:<15.3f} {sub:<15.3f} {diff:>+6.1f}%\n"
f.write(line)
print(line, end='')
f.write("=" * 80 + "\n\n")
print()
# Autocorrelation comparison (critical for temporal structure)
from scipy.stats import pearsonr
# Compute autocorrelation at key lags
lags_to_check = [1, 10, 50, 250, 1250] # Include our identified periods
f.write("Autocorrelation Comparison:\n")
f.write("-" * 80 + "\n")
f.write(f"{'Lag':<10} {'Full Data':<15} {'Subsampled':<15} {'Difference':<15}\n")
f.write("-" * 80 + "\n")
print("Autocorrelation Comparison:")
print("-" * 80)
for lag in lags_to_check:
# Full data autocorrelation
if lag < len(y_train):
acf_full = pearsonr(y_train[:-lag], y_train[lag:])[0]
else:
acf_full = np.nan
# Subsampled autocorrelation (adjust lag for subsampling)
lag_sub = lag // subsample_factor
if lag_sub > 0 and lag_sub < len(y_train_sub):
acf_sub = pearsonr(y_train_sub[:-lag_sub], y_train_sub[lag_sub:])[0]
else:
acf_sub = np.nan
if not np.isnan(acf_full) and not np.isnan(acf_sub):
diff = acf_sub - acf_full
line = f"{lag:<10} {acf_full:<15.3f} {acf_sub:<15.3f} {diff:>+6.3f}\n"
f.write(line)
print(line, end='')
else:
line = f"{lag:<10} {str(acf_full):<15} {str(acf_sub):<15} {'N/A':<15}\n"
f.write(line)
print(line, end='')
f.write("=" * 80 + "\n")
print("=" * 80)
print(f"\n✓ Validation report saved to: {output_file}")
print(f" Use: cat {output_file}")
In [ ]:
Copied!
# Visual validation: overlay subsampled points on full data
fig, axes = plt.subplots(3, 1, figsize=(18, 12))
# Plot 1: First 5000 timesteps - Full data vs subsampled
n_viz = 5000
axes[0].plot(X_train[:n_viz].flatten(), y_train[:n_viz], 'k-', linewidth=0.5, alpha=0.5, label='Full data')
axes[0].scatter(X_train_sub[:n_viz//subsample_factor], y_train_sub[:n_viz//subsample_factor],
c='red', s=20, alpha=0.7, zorder=5, label=f'Subsampled (every {subsample_factor}th)')
axes[0].set_title(f'Subsampling Validation: First {n_viz} Timesteps', fontsize=14, fontweight='bold')
axes[0].set_xlabel('Timestamp')
axes[0].set_ylabel('IOPS Value')
axes[0].legend()
axes[0].grid(alpha=0.3)
# Plot 2: Zoom into 2 periods of fast oscillation (~500 timesteps)
zoom_start = 10000
zoom_end = zoom_start + 500
zoom_indices = (X_train.flatten() >= zoom_start) & (X_train.flatten() < zoom_end)
zoom_indices_sub = (X_train_sub.flatten() >= zoom_start) & (X_train_sub.flatten() < zoom_end)
axes[1].plot(X_train[zoom_indices].flatten(), y_train[zoom_indices], 'k-', linewidth=1.5, alpha=0.7, label='Full data')
axes[1].scatter(X_train_sub[zoom_indices_sub], y_train_sub[zoom_indices_sub],
c='red', s=50, alpha=0.8, zorder=5, label=f'Subsampled points')
axes[1].set_title(f'Zoom: Timesteps {zoom_start}-{zoom_end} (~2 Fast Periods)', fontsize=14, fontweight='bold')
axes[1].set_xlabel('Timestamp')
axes[1].set_ylabel('IOPS Value')
axes[1].legend()
axes[1].grid(alpha=0.3)
# Plot 3: Distribution comparison (histogram + KDE)
axes[2].hist(y_train, bins=50, alpha=0.5, density=True, color='black', label='Full data')
axes[2].hist(y_train_sub, bins=50, alpha=0.5, density=True, color='red', label='Subsampled')
axes[2].set_title('Value Distribution Comparison', fontsize=14, fontweight='bold')
axes[2].set_xlabel('IOPS Value')
axes[2].set_ylabel('Density')
axes[2].legend()
axes[2].grid(alpha=0.3)
plt.tight_layout()
# Save plot
visual_plot_file = '../subsampling_visual_validation.png'
plt.savefig(visual_plot_file, dpi=150, bbox_inches='tight')
print(f"✓ Visual validation plot saved to: {visual_plot_file}")
plt.show()
# Visual validation: overlay subsampled points on full data
fig, axes = plt.subplots(3, 1, figsize=(18, 12))
# Plot 1: First 5000 timesteps - Full data vs subsampled
n_viz = 5000
axes[0].plot(X_train[:n_viz].flatten(), y_train[:n_viz], 'k-', linewidth=0.5, alpha=0.5, label='Full data')
axes[0].scatter(X_train_sub[:n_viz//subsample_factor], y_train_sub[:n_viz//subsample_factor],
c='red', s=20, alpha=0.7, zorder=5, label=f'Subsampled (every {subsample_factor}th)')
axes[0].set_title(f'Subsampling Validation: First {n_viz} Timesteps', fontsize=14, fontweight='bold')
axes[0].set_xlabel('Timestamp')
axes[0].set_ylabel('IOPS Value')
axes[0].legend()
axes[0].grid(alpha=0.3)
# Plot 2: Zoom into 2 periods of fast oscillation (~500 timesteps)
zoom_start = 10000
zoom_end = zoom_start + 500
zoom_indices = (X_train.flatten() >= zoom_start) & (X_train.flatten() < zoom_end)
zoom_indices_sub = (X_train_sub.flatten() >= zoom_start) & (X_train_sub.flatten() < zoom_end)
axes[1].plot(X_train[zoom_indices].flatten(), y_train[zoom_indices], 'k-', linewidth=1.5, alpha=0.7, label='Full data')
axes[1].scatter(X_train_sub[zoom_indices_sub], y_train_sub[zoom_indices_sub],
c='red', s=50, alpha=0.8, zorder=5, label=f'Subsampled points')
axes[1].set_title(f'Zoom: Timesteps {zoom_start}-{zoom_end} (~2 Fast Periods)', fontsize=14, fontweight='bold')
axes[1].set_xlabel('Timestamp')
axes[1].set_ylabel('IOPS Value')
axes[1].legend()
axes[1].grid(alpha=0.3)
# Plot 3: Distribution comparison (histogram + KDE)
axes[2].hist(y_train, bins=50, alpha=0.5, density=True, color='black', label='Full data')
axes[2].hist(y_train_sub, bins=50, alpha=0.5, density=True, color='red', label='Subsampled')
axes[2].set_title('Value Distribution Comparison', fontsize=14, fontweight='bold')
axes[2].set_xlabel('IOPS Value')
axes[2].set_ylabel('Density')
axes[2].legend()
axes[2].grid(alpha=0.3)
plt.tight_layout()
# Save plot
visual_plot_file = '../subsampling_visual_validation.png'
plt.savefig(visual_plot_file, dpi=150, bbox_inches='tight')
print(f"✓ Visual validation plot saved to: {visual_plot_file}")
plt.show()
In [ ]:
Copied!
# Frequency domain analysis: Power spectral density comparison
from scipy import signal
# Compute PSD for full data (use Welch's method for long series)
freqs_full, psd_full = signal.welch(y_train, fs=1.0, nperseg=min(2048, len(y_train)//4))
# Compute PSD for subsampled data
freqs_sub, psd_sub = signal.welch(y_train_sub, fs=1.0/subsample_factor, nperseg=min(2048, len(y_train_sub)//4))
# Plot power spectral density
fig, ax = plt.subplots(figsize=(16, 6))
ax.semilogy(freqs_full, psd_full, 'k-', linewidth=2, alpha=0.7, label='Full data PSD')
ax.semilogy(freqs_sub, psd_sub, 'r--', linewidth=2, alpha=0.7, label='Subsampled PSD')
# Mark expected frequencies
expected_fast_freq = 1.0 / 250 # Fast period
expected_slow_freq = 1.0 / 1250 # Slow period
ax.axvline(expected_fast_freq, color='blue', linestyle=':', linewidth=2, alpha=0.5, label=f'Expected fast freq (1/250)')
ax.axvline(expected_slow_freq, color='green', linestyle=':', linewidth=2, alpha=0.5, label=f'Expected slow freq (1/1250)')
# Nyquist frequency for subsampled data
nyquist_sub = 1.0 / (2 * subsample_factor)
ax.axvline(nyquist_sub, color='orange', linestyle='--', linewidth=2, alpha=0.7, label=f'Subsampled Nyquist (1/{2*subsample_factor})')
ax.set_xlabel('Frequency (cycles/timestep)', fontsize=12)
ax.set_ylabel('Power Spectral Density', fontsize=12)
ax.set_title('Frequency Content: Full vs Subsampled Data', fontsize=14, fontweight='bold')
ax.legend()
ax.grid(alpha=0.3, which='both')
ax.set_xlim([freqs_full[1], 0.05]) # Focus on low frequencies
plt.tight_layout()
# Save PSD plot
psd_plot_file = '../subsampling_psd_analysis.png'
plt.savefig(psd_plot_file, dpi=150, bbox_inches='tight')
print(f"✓ PSD analysis plot saved to: {psd_plot_file}")
plt.show()
# Check if fast frequency is above Nyquist - SAVE TO FILE
aliasing_file = '../subsampling_aliasing_analysis.txt'
with open(aliasing_file, 'w') as f:
f.write("=" * 80 + "\n")
f.write("ALIASING ANALYSIS\n")
f.write("=" * 80 + "\n")
f.write(f"Fast period: 250 timesteps → frequency: {expected_fast_freq:.6f} cycles/timestep\n")
f.write(f"Slow period: 1250 timesteps → frequency: {expected_slow_freq:.6f} cycles/timestep\n")
f.write(f"Subsampling factor: {subsample_factor}\n")
f.write(f"Subsampled Nyquist frequency: {nyquist_sub:.6f} cycles/timestep\n")
f.write("\n")
if expected_fast_freq > nyquist_sub:
verdict = "⚠️ ALIASING DETECTED!"
f.write(verdict + "\n")
f.write(f" Fast frequency ({expected_fast_freq:.6f}) > Nyquist ({nyquist_sub:.6f})\n")
f.write(f" Subsampling CANNOT capture 250-timestep oscillations!\n")
f.write(f" Need subsampling factor ≤ {int(250/2)} to avoid aliasing\n")
else:
verdict = "✓ No aliasing - fast frequency below Nyquist"
f.write(verdict + "\n")
f.write("\n")
f.write("NOTE: No aliasing doesn't guarantee signal preservation!\n")
f.write("Check visual validation plots to confirm pattern capture.\n")
f.write("=" * 80 + "\n")
# Print to console as well
print()
with open(aliasing_file, 'r') as f:
print(f.read())
print(f"\n✓ Aliasing analysis saved to: {aliasing_file}")
print(f" Use: cat {aliasing_file}")
# Frequency domain analysis: Power spectral density comparison
from scipy import signal
# Compute PSD for full data (use Welch's method for long series)
freqs_full, psd_full = signal.welch(y_train, fs=1.0, nperseg=min(2048, len(y_train)//4))
# Compute PSD for subsampled data
freqs_sub, psd_sub = signal.welch(y_train_sub, fs=1.0/subsample_factor, nperseg=min(2048, len(y_train_sub)//4))
# Plot power spectral density
fig, ax = plt.subplots(figsize=(16, 6))
ax.semilogy(freqs_full, psd_full, 'k-', linewidth=2, alpha=0.7, label='Full data PSD')
ax.semilogy(freqs_sub, psd_sub, 'r--', linewidth=2, alpha=0.7, label='Subsampled PSD')
# Mark expected frequencies
expected_fast_freq = 1.0 / 250 # Fast period
expected_slow_freq = 1.0 / 1250 # Slow period
ax.axvline(expected_fast_freq, color='blue', linestyle=':', linewidth=2, alpha=0.5, label=f'Expected fast freq (1/250)')
ax.axvline(expected_slow_freq, color='green', linestyle=':', linewidth=2, alpha=0.5, label=f'Expected slow freq (1/1250)')
# Nyquist frequency for subsampled data
nyquist_sub = 1.0 / (2 * subsample_factor)
ax.axvline(nyquist_sub, color='orange', linestyle='--', linewidth=2, alpha=0.7, label=f'Subsampled Nyquist (1/{2*subsample_factor})')
ax.set_xlabel('Frequency (cycles/timestep)', fontsize=12)
ax.set_ylabel('Power Spectral Density', fontsize=12)
ax.set_title('Frequency Content: Full vs Subsampled Data', fontsize=14, fontweight='bold')
ax.legend()
ax.grid(alpha=0.3, which='both')
ax.set_xlim([freqs_full[1], 0.05]) # Focus on low frequencies
plt.tight_layout()
# Save PSD plot
psd_plot_file = '../subsampling_psd_analysis.png'
plt.savefig(psd_plot_file, dpi=150, bbox_inches='tight')
print(f"✓ PSD analysis plot saved to: {psd_plot_file}")
plt.show()
# Check if fast frequency is above Nyquist - SAVE TO FILE
aliasing_file = '../subsampling_aliasing_analysis.txt'
with open(aliasing_file, 'w') as f:
f.write("=" * 80 + "\n")
f.write("ALIASING ANALYSIS\n")
f.write("=" * 80 + "\n")
f.write(f"Fast period: 250 timesteps → frequency: {expected_fast_freq:.6f} cycles/timestep\n")
f.write(f"Slow period: 1250 timesteps → frequency: {expected_slow_freq:.6f} cycles/timestep\n")
f.write(f"Subsampling factor: {subsample_factor}\n")
f.write(f"Subsampled Nyquist frequency: {nyquist_sub:.6f} cycles/timestep\n")
f.write("\n")
if expected_fast_freq > nyquist_sub:
verdict = "⚠️ ALIASING DETECTED!"
f.write(verdict + "\n")
f.write(f" Fast frequency ({expected_fast_freq:.6f}) > Nyquist ({nyquist_sub:.6f})\n")
f.write(f" Subsampling CANNOT capture 250-timestep oscillations!\n")
f.write(f" Need subsampling factor ≤ {int(250/2)} to avoid aliasing\n")
else:
verdict = "✓ No aliasing - fast frequency below Nyquist"
f.write(verdict + "\n")
f.write("\n")
f.write("NOTE: No aliasing doesn't guarantee signal preservation!\n")
f.write("Check visual validation plots to confirm pattern capture.\n")
f.write("=" * 80 + "\n")
# Print to console as well
print()
with open(aliasing_file, 'r') as f:
print(f.read())
print(f"\n✓ Aliasing analysis saved to: {aliasing_file}")
print(f" Use: cat {aliasing_file}")
In [ ]:
Copied!
# Normalize subsampled data (same normalization as full dataset)
X_train_sub_norm = (X_train_sub - X_min) / X_range
# Convert to PyTorch tensors
X_train_sub_t = torch.tensor(X_train_sub_norm, dtype=dtype, device=device)
y_train_sub_t = torch.tensor(y_train_sub, dtype=dtype, device=device)
print(f"Subsampled tensor shapes:")
print(f" X: {X_train_sub_t.shape}")
print(f" y: {y_train_sub_t.shape}")
print(f" Normalized range: [{X_train_sub_norm.min():.6f}, {X_train_sub_norm.max():.6f}]")
# Normalize subsampled data (same normalization as full dataset)
X_train_sub_norm = (X_train_sub - X_min) / X_range
# Convert to PyTorch tensors
X_train_sub_t = torch.tensor(X_train_sub_norm, dtype=dtype, device=device)
y_train_sub_t = torch.tensor(y_train_sub, dtype=dtype, device=device)
print(f"Subsampled tensor shapes:")
print(f" X: {X_train_sub_t.shape}")
print(f" y: {y_train_sub_t.shape}")
print(f" Normalized range: [{X_train_sub_norm.min():.6f}, {X_train_sub_norm.max():.6f}]")
4.2 Exact GP Model Definition¶
In [ ]:
Copied!
from gpytorch.models import ExactGP
from gpytorch.means import ConstantMean
from gpytorch.kernels import ScaleKernel, RBFKernel, PeriodicKernel
from gpytorch.distributions import MultivariateNormal
from gpytorch.mlls import ExactMarginalLogLikelihood
from gpytorch.likelihoods import GaussianLikelihood
class ExactGPModel(ExactGP):
"""
Exact GP model (no inducing points) - follows GPyTorch simple example pattern.
This is computationally expensive (O(n³)) but can capture fine structure
without variational approximation.
"""
def __init__(self, train_x, train_y, likelihood, kernel_type='rbf'):
super().__init__(train_x, train_y, likelihood)
# Mean function - initialize to data mean
self.mean_module = ConstantMean()
# Kernel selection
if kernel_type == 'rbf':
# Simple RBF kernel (baseline)
self.covar_module = ScaleKernel(RBFKernel())
elif kernel_type == 'periodic':
# RBF + Periodic (multi-scale patterns)
slow_period = 1250 / X_range
fast_period = 250 / X_range
slow_periodic = ScaleKernel(PeriodicKernel())
slow_periodic.base_kernel.period_length = slow_period
slow_periodic.base_kernel.raw_period_length.requires_grad = False
fast_periodic = ScaleKernel(PeriodicKernel())
fast_periodic.base_kernel.period_length = fast_period
fast_periodic.base_kernel.raw_period_length.requires_grad = False
rbf_component = ScaleKernel(RBFKernel())
# Additive kernel
self.covar_module = slow_periodic + fast_periodic + rbf_component
else:
raise ValueError(f"Unknown kernel_type: {kernel_type}")
def forward(self, x):
mean_x = self.mean_module(x)
covar_x = self.covar_module(x)
return MultivariateNormal(mean_x, covar_x)
print("✓ Exact GP model class defined")
from gpytorch.models import ExactGP
from gpytorch.means import ConstantMean
from gpytorch.kernels import ScaleKernel, RBFKernel, PeriodicKernel
from gpytorch.distributions import MultivariateNormal
from gpytorch.mlls import ExactMarginalLogLikelihood
from gpytorch.likelihoods import GaussianLikelihood
class ExactGPModel(ExactGP):
"""
Exact GP model (no inducing points) - follows GPyTorch simple example pattern.
This is computationally expensive (O(n³)) but can capture fine structure
without variational approximation.
"""
def __init__(self, train_x, train_y, likelihood, kernel_type='rbf'):
super().__init__(train_x, train_y, likelihood)
# Mean function - initialize to data mean
self.mean_module = ConstantMean()
# Kernel selection
if kernel_type == 'rbf':
# Simple RBF kernel (baseline)
self.covar_module = ScaleKernel(RBFKernel())
elif kernel_type == 'periodic':
# RBF + Periodic (multi-scale patterns)
slow_period = 1250 / X_range
fast_period = 250 / X_range
slow_periodic = ScaleKernel(PeriodicKernel())
slow_periodic.base_kernel.period_length = slow_period
slow_periodic.base_kernel.raw_period_length.requires_grad = False
fast_periodic = ScaleKernel(PeriodicKernel())
fast_periodic.base_kernel.period_length = fast_period
fast_periodic.base_kernel.raw_period_length.requires_grad = False
rbf_component = ScaleKernel(RBFKernel())
# Additive kernel
self.covar_module = slow_periodic + fast_periodic + rbf_component
else:
raise ValueError(f"Unknown kernel_type: {kernel_type}")
def forward(self, x):
mean_x = self.mean_module(x)
covar_x = self.covar_module(x)
return MultivariateNormal(mean_x, covar_x)
print("✓ Exact GP model class defined")
4.3 Train Exact GP with RBF Kernel (Baseline)¶
In [ ]:
Copied!
# Initialize model and likelihood
likelihood_exact_rbf = GaussianLikelihood()
model_exact_rbf = ExactGPModel(X_train_sub_t, y_train_sub_t, likelihood_exact_rbf, kernel_type='rbf')
# Initialize mean to data mean
model_exact_rbf.mean_module.constant.data.fill_(y_train_sub.mean())
# Move to device
model_exact_rbf = model_exact_rbf.to(device)
likelihood_exact_rbf = likelihood_exact_rbf.to(device)
print("✓ Exact GP (RBF kernel) initialized")
print(f" Kernel: RBF")
print(f" Training samples: {len(X_train_sub_t):,}")
print(f" Mean initialized to: {model_exact_rbf.mean_module.constant.item():.3f}")
# Initialize model and likelihood
likelihood_exact_rbf = GaussianLikelihood()
model_exact_rbf = ExactGPModel(X_train_sub_t, y_train_sub_t, likelihood_exact_rbf, kernel_type='rbf')
# Initialize mean to data mean
model_exact_rbf.mean_module.constant.data.fill_(y_train_sub.mean())
# Move to device
model_exact_rbf = model_exact_rbf.to(device)
likelihood_exact_rbf = likelihood_exact_rbf.to(device)
print("✓ Exact GP (RBF kernel) initialized")
print(f" Kernel: RBF")
print(f" Training samples: {len(X_train_sub_t):,}")
print(f" Mean initialized to: {model_exact_rbf.mean_module.constant.item():.3f}")
In [ ]:
Copied!
# Training loop - standard exact GP training (like GPyTorch example)
model_exact_rbf.train()
likelihood_exact_rbf.train()
# Use Adam optimizer
optimizer = torch.optim.Adam(model_exact_rbf.parameters(), lr=0.1)
# Marginal log likelihood
mll = ExactMarginalLogLikelihood(likelihood_exact_rbf, model_exact_rbf)
# Training
n_epochs_exact = 50
losses_exact_rbf = []
print(f"Training exact GP (RBF) for {n_epochs_exact} epochs...")
for epoch in range(n_epochs_exact):
optimizer.zero_grad()
output = model_exact_rbf(X_train_sub_t)
loss = -mll(output, y_train_sub_t)
loss.backward()
optimizer.step()
losses_exact_rbf.append(loss.item())
if (epoch + 1) % 10 == 0:
print(f" Epoch {epoch+1}/{n_epochs_exact} - Loss: {loss.item():.3f}")
print(f"✓ Training complete - Final loss: {losses_exact_rbf[-1]:.3f}")
# Training loop - standard exact GP training (like GPyTorch example)
model_exact_rbf.train()
likelihood_exact_rbf.train()
# Use Adam optimizer
optimizer = torch.optim.Adam(model_exact_rbf.parameters(), lr=0.1)
# Marginal log likelihood
mll = ExactMarginalLogLikelihood(likelihood_exact_rbf, model_exact_rbf)
# Training
n_epochs_exact = 50
losses_exact_rbf = []
print(f"Training exact GP (RBF) for {n_epochs_exact} epochs...")
for epoch in range(n_epochs_exact):
optimizer.zero_grad()
output = model_exact_rbf(X_train_sub_t)
loss = -mll(output, y_train_sub_t)
loss.backward()
optimizer.step()
losses_exact_rbf.append(loss.item())
if (epoch + 1) % 10 == 0:
print(f" Epoch {epoch+1}/{n_epochs_exact} - Loss: {loss.item():.3f}")
print(f"✓ Training complete - Final loss: {losses_exact_rbf[-1]:.3f}")
4.4 Evaluate Exact GP (RBF Baseline)¶
In [ ]:
Copied!
# Make predictions on test set
model_exact_rbf.eval()
likelihood_exact_rbf.eval()
with torch.no_grad(), gpytorch.settings.fast_pred_var():
pred_dist_rbf = likelihood_exact_rbf(model_exact_rbf(X_test_t))
mean_exact_rbf = pred_dist_rbf.mean.cpu().numpy()
std_exact_rbf = pred_dist_rbf.stddev.cpu().numpy()
# Compute metrics
from sklearn.metrics import mean_squared_error, mean_absolute_error, r2_score
rmse_rbf = np.sqrt(mean_squared_error(y_test, mean_exact_rbf))
mae_rbf = mean_absolute_error(y_test, mean_exact_rbf)
r2_rbf = r2_score(y_test, mean_exact_rbf)
print("=" * 70)
print("EXACT GP (RBF KERNEL) - BASELINE RESULTS")
print("=" * 70)
print(f"RMSE: {rmse_rbf:.3f}")
print(f"MAE: {mae_rbf:.3f}")
print(f"R²: {r2_rbf:.3f}")
print()
print(f"Prediction range: [{mean_exact_rbf.min():.2f}, {mean_exact_rbf.max():.2f}]")
print(f"Data range: [{y_test.min():.2f}, {y_test.max():.2f}]")
print(f"Prediction variance: {mean_exact_rbf.var():.3f}")
print(f"Data variance: {y_test.var():.3f}")
print("=" * 70)
# Make predictions on test set
model_exact_rbf.eval()
likelihood_exact_rbf.eval()
with torch.no_grad(), gpytorch.settings.fast_pred_var():
pred_dist_rbf = likelihood_exact_rbf(model_exact_rbf(X_test_t))
mean_exact_rbf = pred_dist_rbf.mean.cpu().numpy()
std_exact_rbf = pred_dist_rbf.stddev.cpu().numpy()
# Compute metrics
from sklearn.metrics import mean_squared_error, mean_absolute_error, r2_score
rmse_rbf = np.sqrt(mean_squared_error(y_test, mean_exact_rbf))
mae_rbf = mean_absolute_error(y_test, mean_exact_rbf)
r2_rbf = r2_score(y_test, mean_exact_rbf)
print("=" * 70)
print("EXACT GP (RBF KERNEL) - BASELINE RESULTS")
print("=" * 70)
print(f"RMSE: {rmse_rbf:.3f}")
print(f"MAE: {mae_rbf:.3f}")
print(f"R²: {r2_rbf:.3f}")
print()
print(f"Prediction range: [{mean_exact_rbf.min():.2f}, {mean_exact_rbf.max():.2f}]")
print(f"Data range: [{y_test.min():.2f}, {y_test.max():.2f}]")
print(f"Prediction variance: {mean_exact_rbf.var():.3f}")
print(f"Data variance: {y_test.var():.3f}")
print("=" * 70)
In [ ]:
Copied!
# Visualize first 1000 predictions
n_viz = 1000
fig, axes = plt.subplots(2, 1, figsize=(16, 10))
# Plot 1: Predictions vs actual
axes[0].plot(y_test[:n_viz], 'k-', linewidth=1, label='Actual', alpha=0.7)
axes[0].plot(mean_exact_rbf[:n_viz], 'b--', linewidth=1.5, label='Exact GP (RBF) predictions')
axes[0].fill_between(
np.arange(n_viz),
mean_exact_rbf[:n_viz] - 2*std_exact_rbf[:n_viz],
mean_exact_rbf[:n_viz] + 2*std_exact_rbf[:n_viz],
alpha=0.2,
color='blue',
label='95% CI'
)
axes[0].set_title('Exact GP (RBF): First 1000 Test Points', fontsize=14, fontweight='bold')
axes[0].set_xlabel('Time step')
axes[0].set_ylabel('IOPS Value')
axes[0].legend()
axes[0].grid(alpha=0.3)
# Plot 2: Residuals
residuals_rbf = y_test[:n_viz] - mean_exact_rbf[:n_viz]
axes[1].scatter(np.arange(n_viz), residuals_rbf, alpha=0.3, s=10, color='red')
axes[1].axhline(y=0, color='k', linestyle='--', linewidth=2)
axes[1].set_title('Residuals', fontsize=14, fontweight='bold')
axes[1].set_xlabel('Time step')
axes[1].set_ylabel('Residual (Actual - Predicted)')
axes[1].grid(alpha=0.3)
plt.tight_layout()
plt.show()
# Visualize first 1000 predictions
n_viz = 1000
fig, axes = plt.subplots(2, 1, figsize=(16, 10))
# Plot 1: Predictions vs actual
axes[0].plot(y_test[:n_viz], 'k-', linewidth=1, label='Actual', alpha=0.7)
axes[0].plot(mean_exact_rbf[:n_viz], 'b--', linewidth=1.5, label='Exact GP (RBF) predictions')
axes[0].fill_between(
np.arange(n_viz),
mean_exact_rbf[:n_viz] - 2*std_exact_rbf[:n_viz],
mean_exact_rbf[:n_viz] + 2*std_exact_rbf[:n_viz],
alpha=0.2,
color='blue',
label='95% CI'
)
axes[0].set_title('Exact GP (RBF): First 1000 Test Points', fontsize=14, fontweight='bold')
axes[0].set_xlabel('Time step')
axes[0].set_ylabel('IOPS Value')
axes[0].legend()
axes[0].grid(alpha=0.3)
# Plot 2: Residuals
residuals_rbf = y_test[:n_viz] - mean_exact_rbf[:n_viz]
axes[1].scatter(np.arange(n_viz), residuals_rbf, alpha=0.3, s=10, color='red')
axes[1].axhline(y=0, color='k', linestyle='--', linewidth=2)
axes[1].set_title('Residuals', fontsize=14, fontweight='bold')
axes[1].set_xlabel('Time step')
axes[1].set_ylabel('Residual (Actual - Predicted)')
axes[1].grid(alpha=0.3)
plt.tight_layout()
plt.show()
5. Sparse Variational GP Approach (Original - For Comparison)¶
Approach Comparison¶
| Aspect | Robust (Recommended) | Traditional (Baseline) |
|---|---|---|
| Training Data | ALL 146,255 samples | 145,970 samples (exclude anomalies) |
| Likelihood | Student-t (ν=4) | Gaussian |
| Philosophy | Outliers are real data | Outliers corrupt training |
| Robustness | Heavy tails handle extremes | Assumes normality |
| Production | Trained on operational reality | Trained on sanitized data |
Why Student-t Likelihood?¶
The Student-t distribution has heavy tails controlled by degrees of freedom (ν):
$$ p(y | \mu, \sigma, \nu) = \frac{\Gamma(\frac{\nu+1}{2})}{\Gamma(\frac{\nu}{2})\sqrt{\pi\nu}\sigma} \left(1 + \frac{1}{\nu}\left(\frac{y-\mu}{\sigma}\right)^2\right)^{-\frac{\nu+1}{2}} $$
- ν=4: Good balance between robustness and efficiency
- Lower ν → Heavier tails → More robust to outliers
- As ν→∞ → Converges to Gaussian
Key advantage: Outliers contribute less to parameter learning, allowing the model to learn patterns despite occasional anomalies.
5. Using the Gaussian Process Library¶
NOTE: The implementation details (kernel design, model architecture, training loops) have been extracted to the cloud_sim.ml_models.gaussian_process library for production use.
For design rationale and implementation details, see: docs/modeling/gaussian-process-design.md
5.1 Library Architecture Overview¶
The GP module provides:
CompositePeriodicKernel: Multi-scale periodic kernel- SLOW component: Sawtooth envelope (1250 steps ≈ 21 hours)
- FAST component: Sinusoidal carrier (250 steps ≈ 4 hours)
- ADDITIVE structure for numerical stability
SparseGPModel: Variational sparse GP with O(nm²) complexity- Uses inducing points for scalability
- Learns inducing locations during training
- Cholesky variational distribution
Training utilities:
train_gp_model,save_model,load_model- Mini-batch training with ELBO objective
- Maximum numerical stability settings
- Progress tracking and model persistence
Evaluation metrics:
compute_metrics,compute_anomaly_metrics- Point accuracy (RMSE, MAE, R²)
- Uncertainty calibration (coverage, sharpness)
- Anomaly detection (precision, recall, F1, AUC-ROC)
5.2 Initialize Models Using Library¶
In [ ]:
Copied!
# Initialize inducing points (evenly spaced across training data)
M = 200 # Number of inducing points
inducing_points = initialize_inducing_points(
X_train=X_train_t,
num_inducing=M,
method="evenly_spaced"
)
print(f"✓ Inducing points initialized:")
print(f" Count: {M}")
print(f" Shape: {inducing_points.shape}")
print(f" Range: {inducing_points.min():.6f} → {inducing_points.max():.6f}")
# Initialize inducing points (evenly spaced across training data)
M = 200 # Number of inducing points
inducing_points = initialize_inducing_points(
X_train=X_train_t,
num_inducing=M,
method="evenly_spaced"
)
print(f"✓ Inducing points initialized:")
print(f" Count: {M}")
print(f" Shape: {inducing_points.shape}")
print(f" Range: {inducing_points.min():.6f} → {inducing_points.max():.6f}")
In [ ]:
Copied!
# Compute training data statistics for proper initialization
y_train_mean = y_train.mean()
y_train_std = y_train.std()
y_train_var = y_train.var()
print(f"Training data statistics:")
print(f" Mean: {y_train_mean:.3f}")
print(f" Std: {y_train_std:.3f}")
print(f" Variance: {y_train_var:.3f}")
print()
print("⚠️ CRITICAL: Proper initialization required for variational GP!")
print(" Default initialization (mean=0, outputscale=1) causes scale mismatch")
print(" with variational inference's KL penalty, leading to underfitting.")
# Compute training data statistics for proper initialization
y_train_mean = y_train.mean()
y_train_std = y_train.std()
y_train_var = y_train.var()
print(f"Training data statistics:")
print(f" Mean: {y_train_mean:.3f}")
print(f" Std: {y_train_std:.3f}")
print(f" Variance: {y_train_var:.3f}")
print()
print("⚠️ CRITICAL: Proper initialization required for variational GP!")
print(" Default initialization (mean=0, outputscale=1) causes scale mismatch")
print(" with variational inference's KL penalty, leading to underfitting.")
In [ ]:
Copied!
# Create robust model with Student-t likelihood
model_robust = SparseGPModel(
inducing_points=inducing_points,
slow_period=1250 / X_range, # Sawtooth envelope period
fast_period=250 / X_range, # Sinusoidal carrier period
rbf_lengthscale=0.1
).to(device)
likelihood_robust = StudentTLikelihood(
deg_free_prior=gpytorch.priors.NormalPrior(4.0, 1.0)
).to(device)
# Initialize mean function to training data mean
model_robust.mean_module.constant.data.fill_(y_train_mean)
# Initialize kernel outputscales based on data variance
# Distribute variance across the three kernel components
variance_per_component = y_train_var / 3.0
initial_outputscale = variance_per_component
model_robust.covar_module.slow_periodic.outputscale = initial_outputscale
model_robust.covar_module.fast_periodic.outputscale = initial_outputscale
model_robust.covar_module.rbf.outputscale = initial_outputscale
print("✓ Robust model initialized (Student-t likelihood)")
print(f" Training on ALL {len(X_train_t):,} samples (including anomalies)")
print(f" Mean initialized to: {model_robust.mean_module.constant.item():.3f}")
print(f" Outputscales initialized to: {initial_outputscale:.3f} each")
# Create robust model with Student-t likelihood
model_robust = SparseGPModel(
inducing_points=inducing_points,
slow_period=1250 / X_range, # Sawtooth envelope period
fast_period=250 / X_range, # Sinusoidal carrier period
rbf_lengthscale=0.1
).to(device)
likelihood_robust = StudentTLikelihood(
deg_free_prior=gpytorch.priors.NormalPrior(4.0, 1.0)
).to(device)
# Initialize mean function to training data mean
model_robust.mean_module.constant.data.fill_(y_train_mean)
# Initialize kernel outputscales based on data variance
# Distribute variance across the three kernel components
variance_per_component = y_train_var / 3.0
initial_outputscale = variance_per_component
model_robust.covar_module.slow_periodic.outputscale = initial_outputscale
model_robust.covar_module.fast_periodic.outputscale = initial_outputscale
model_robust.covar_module.rbf.outputscale = initial_outputscale
print("✓ Robust model initialized (Student-t likelihood)")
print(f" Training on ALL {len(X_train_t):,} samples (including anomalies)")
print(f" Mean initialized to: {model_robust.mean_module.constant.item():.3f}")
print(f" Outputscales initialized to: {initial_outputscale:.3f} each")
In [ ]:
Copied!
# Compute clean training data statistics for proper initialization
y_train_clean_mean = y_train_clean.mean()
y_train_clean_var = y_train_clean.var()
print(f"Clean training data statistics:")
print(f" Mean: {y_train_clean_mean:.3f}")
print(f" Variance: {y_train_clean_var:.3f}")
# Compute clean training data statistics for proper initialization
y_train_clean_mean = y_train_clean.mean()
y_train_clean_var = y_train_clean.var()
print(f"Clean training data statistics:")
print(f" Mean: {y_train_clean_mean:.3f}")
print(f" Variance: {y_train_clean_var:.3f}")
In [ ]:
Copied!
# Create traditional model with Gaussian likelihood (baseline)
inducing_points_clean = initialize_inducing_points(
X_train=X_train_clean_t,
num_inducing=M,
method="evenly_spaced"
)
model_traditional = SparseGPModel(
inducing_points=inducing_points_clean,
slow_period=1250 / X_range,
fast_period=250 / X_range,
rbf_lengthscale=0.1
).to(device)
likelihood_traditional = GaussianLikelihood().to(device)
# Initialize mean function to clean training data mean
model_traditional.mean_module.constant.data.fill_(y_train_clean_mean)
# Initialize kernel outputscales based on clean data variance
variance_per_component_clean = y_train_clean_var / 3.0
model_traditional.covar_module.slow_periodic.outputscale = variance_per_component_clean
model_traditional.covar_module.fast_periodic.outputscale = variance_per_component_clean
model_traditional.covar_module.rbf.outputscale = variance_per_component_clean
print("✓ Traditional model initialized (Gaussian likelihood)")
print(f" Training on {len(X_train_clean_t):,} samples (excluded {n_anomalies_train} anomalies)")
print(f" Mean initialized to: {model_traditional.mean_module.constant.item():.3f}")
print(f" Outputscales initialized to: {variance_per_component_clean:.3f} each")
# Create traditional model with Gaussian likelihood (baseline)
inducing_points_clean = initialize_inducing_points(
X_train=X_train_clean_t,
num_inducing=M,
method="evenly_spaced"
)
model_traditional = SparseGPModel(
inducing_points=inducing_points_clean,
slow_period=1250 / X_range,
fast_period=250 / X_range,
rbf_lengthscale=0.1
).to(device)
likelihood_traditional = GaussianLikelihood().to(device)
# Initialize mean function to clean training data mean
model_traditional.mean_module.constant.data.fill_(y_train_clean_mean)
# Initialize kernel outputscales based on clean data variance
variance_per_component_clean = y_train_clean_var / 3.0
model_traditional.covar_module.slow_periodic.outputscale = variance_per_component_clean
model_traditional.covar_module.fast_periodic.outputscale = variance_per_component_clean
model_traditional.covar_module.rbf.outputscale = variance_per_component_clean
print("✓ Traditional model initialized (Gaussian likelihood)")
print(f" Training on {len(X_train_clean_t):,} samples (excluded {n_anomalies_train} anomalies)")
print(f" Mean initialized to: {model_traditional.mean_module.constant.item():.3f}")
print(f" Outputscales initialized to: {variance_per_component_clean:.3f} each")
In [ ]:
Copied!
import os
# Check if saved models exist
robust_model_path = '../models/gp_robust_model.pth'
traditional_model_path = '../models/gp_traditional_model.pth'
models_exist = os.path.exists(robust_model_path) and os.path.exists(traditional_model_path)
if models_exist:
print("✓ Found saved models - loading from disk...")
print(f" Robust: {robust_model_path}")
print(f" Traditional: {traditional_model_path}")
print("\nTo retrain from scratch, delete the models/ directory")
else:
print("✗ No saved models found - will train from scratch")
print("\nModels will be saved to ../models/ after training")
import os
# Check if saved models exist
robust_model_path = '../models/gp_robust_model.pth'
traditional_model_path = '../models/gp_traditional_model.pth'
models_exist = os.path.exists(robust_model_path) and os.path.exists(traditional_model_path)
if models_exist:
print("✓ Found saved models - loading from disk...")
print(f" Robust: {robust_model_path}")
print(f" Traditional: {traditional_model_path}")
print("\nTo retrain from scratch, delete the models/ directory")
else:
print("✗ No saved models found - will train from scratch")
print("\nModels will be saved to ../models/ after training")
6.2 Train or Load Robust Model (Student-t Likelihood)¶
In [ ]:
Copied!
# Load model if it exists, otherwise train from scratch
if models_exist:
model_robust, likelihood_robust, checkpoint_robust = load_model(
load_path=robust_model_path,
likelihood_class=StudentTLikelihood,
device=device,
# Kernel parameters for backward compatibility with old checkpoints
slow_period=1250 / X_range,
fast_period=250 / X_range,
rbf_lengthscale=0.1
)
losses_robust = checkpoint_robust['losses']
else:
# Train robust model using library (with device-specific jitter)
losses_robust = train_gp_model(
model=model_robust,
likelihood=likelihood_robust,
X_train=X_train_t,
y_train=y_train_t,
n_epochs=100,
batch_size=2048,
learning_rate=0.01,
cholesky_jitter=cholesky_jitter,
cholesky_max_tries=cholesky_max_tries,
verbose=True
)
# Save trained model
save_model(
model=model_robust,
likelihood=likelihood_robust,
save_path=robust_model_path,
losses=losses_robust,
metadata={'dataset': 'IOPS', 'approach': 'robust'}
)
# Load model if it exists, otherwise train from scratch
if models_exist:
model_robust, likelihood_robust, checkpoint_robust = load_model(
load_path=robust_model_path,
likelihood_class=StudentTLikelihood,
device=device,
# Kernel parameters for backward compatibility with old checkpoints
slow_period=1250 / X_range,
fast_period=250 / X_range,
rbf_lengthscale=0.1
)
losses_robust = checkpoint_robust['losses']
else:
# Train robust model using library (with device-specific jitter)
losses_robust = train_gp_model(
model=model_robust,
likelihood=likelihood_robust,
X_train=X_train_t,
y_train=y_train_t,
n_epochs=100,
batch_size=2048,
learning_rate=0.01,
cholesky_jitter=cholesky_jitter,
cholesky_max_tries=cholesky_max_tries,
verbose=True
)
# Save trained model
save_model(
model=model_robust,
likelihood=likelihood_robust,
save_path=robust_model_path,
losses=losses_robust,
metadata={'dataset': 'IOPS', 'approach': 'robust'}
)
6.3 Train or Load Traditional Model (Gaussian Likelihood)¶
In [ ]:
Copied!
# Load model if it exists, otherwise train from scratch
if models_exist:
model_traditional, likelihood_traditional, checkpoint_trad = load_model(
load_path=traditional_model_path,
likelihood_class=GaussianLikelihood,
device=device,
# Kernel parameters for backward compatibility with old checkpoints
slow_period=1250 / X_range,
fast_period=250 / X_range,
rbf_lengthscale=0.1
)
losses_traditional = checkpoint_trad['losses']
else:
# Train traditional model using library (with device-specific jitter)
losses_traditional = train_gp_model(
model=model_traditional,
likelihood=likelihood_traditional,
X_train=X_train_clean_t,
y_train=y_train_clean_t,
n_epochs=100,
batch_size=2048,
learning_rate=0.01,
cholesky_jitter=cholesky_jitter,
cholesky_max_tries=cholesky_max_tries,
verbose=True
)
# Save trained model
save_model(
model=model_traditional,
likelihood=likelihood_traditional,
save_path=traditional_model_path,
losses=losses_traditional,
metadata={'dataset': 'IOPS', 'approach': 'traditional'}
)
# Load model if it exists, otherwise train from scratch
if models_exist:
model_traditional, likelihood_traditional, checkpoint_trad = load_model(
load_path=traditional_model_path,
likelihood_class=GaussianLikelihood,
device=device,
# Kernel parameters for backward compatibility with old checkpoints
slow_period=1250 / X_range,
fast_period=250 / X_range,
rbf_lengthscale=0.1
)
losses_traditional = checkpoint_trad['losses']
else:
# Train traditional model using library (with device-specific jitter)
losses_traditional = train_gp_model(
model=model_traditional,
likelihood=likelihood_traditional,
X_train=X_train_clean_t,
y_train=y_train_clean_t,
n_epochs=100,
batch_size=2048,
learning_rate=0.01,
cholesky_jitter=cholesky_jitter,
cholesky_max_tries=cholesky_max_tries,
verbose=True
)
# Save trained model
save_model(
model=model_traditional,
likelihood=likelihood_traditional,
save_path=traditional_model_path,
losses=losses_traditional,
metadata={'dataset': 'IOPS', 'approach': 'traditional'}
)
6.4 Training Loss Comparison¶
In [ ]:
Copied!
# Plot training losses
fig, ax = plt.subplots(figsize=(12, 5))
ax.plot(losses_robust, linewidth=2, label='Robust (Student-t, all data)', color='steelblue')
ax.plot(losses_traditional, linewidth=2, label='Traditional (Gaussian, clean data)', color='coral')
ax.set_xlabel('Epoch', fontsize=12)
ax.set_ylabel('Negative ELBO Loss', fontsize=12)
ax.set_title('Training Loss Comparison', fontsize=14, fontweight='bold')
ax.legend(fontsize=11)
ax.grid(alpha=0.3)
plt.tight_layout()
plt.show()
# Plot training losses
fig, ax = plt.subplots(figsize=(12, 5))
ax.plot(losses_robust, linewidth=2, label='Robust (Student-t, all data)', color='steelblue')
ax.plot(losses_traditional, linewidth=2, label='Traditional (Gaussian, clean data)', color='coral')
ax.set_xlabel('Epoch', fontsize=12)
ax.set_ylabel('Negative ELBO Loss', fontsize=12)
ax.set_title('Training Loss Comparison', fontsize=14, fontweight='bold')
ax.legend(fontsize=11)
ax.grid(alpha=0.3)
plt.tight_layout()
plt.show()
In [ ]:
Copied!
# Switch to evaluation mode
model_robust.eval()
likelihood_robust.eval()
print("Generating predictions with robust model...")
print(f"Test samples: {len(X_test_t):,}")
# Batched prediction to prevent memory exhaustion
batch_size_pred = 4096 # Process in chunks
n_batches = int(np.ceil(len(X_test_t) / batch_size_pred))
mean_robust_list = []
std_robust_list = []
# Apply same numerical stability settings as training (device-specific)
with torch.no_grad(), \
gpytorch.settings.fast_pred_var(), \
gpytorch.settings.cholesky_jitter(cholesky_jitter), \
gpytorch.settings.cholesky_max_tries(cholesky_max_tries):
for batch_idx in range(n_batches):
start_idx = batch_idx * batch_size_pred
end_idx = min(start_idx + batch_size_pred, len(X_test_t))
X_batch = X_test_t[start_idx:end_idx]
# Get predictive distribution for batch
pred_dist_batch = likelihood_robust(model_robust(X_batch))
# Extract mean and variance (move to CPU, then to numpy)
mean_batch = pred_dist_batch.mean.cpu().numpy()
std_batch = pred_dist_batch.stddev.cpu().numpy()
# Flatten to 1D (handle any batch/feature dimensions)
mean_batch = mean_batch.flatten()
std_batch = std_batch.flatten()
# Sanity check: should match batch size
expected_size = end_idx - start_idx
if mean_batch.size != expected_size:
print(f" WARNING: Batch {batch_idx} size mismatch: got {mean_batch.size}, expected {expected_size}")
# Truncate or pad to expected size
mean_batch = mean_batch[:expected_size]
std_batch = std_batch[:expected_size]
mean_robust_list.append(mean_batch)
std_robust_list.append(std_batch)
if (batch_idx + 1) % 10 == 0 or (batch_idx + 1) == n_batches:
print(f" Processed {end_idx:,}/{len(X_test_t):,} samples ({100*end_idx/len(X_test_t):.1f}%)")
# Concatenate all batches
mean_robust = np.concatenate(mean_robust_list)
std_robust = np.concatenate(std_robust_list)
print(f"✓ Predictions generated")
print(f" Mean range: {mean_robust.min():.2f} → {mean_robust.max():.2f}")
print(f" Std range: {std_robust.min():.2f} → {std_robust.max():.2f}")
# Switch to evaluation mode
model_robust.eval()
likelihood_robust.eval()
print("Generating predictions with robust model...")
print(f"Test samples: {len(X_test_t):,}")
# Batched prediction to prevent memory exhaustion
batch_size_pred = 4096 # Process in chunks
n_batches = int(np.ceil(len(X_test_t) / batch_size_pred))
mean_robust_list = []
std_robust_list = []
# Apply same numerical stability settings as training (device-specific)
with torch.no_grad(), \
gpytorch.settings.fast_pred_var(), \
gpytorch.settings.cholesky_jitter(cholesky_jitter), \
gpytorch.settings.cholesky_max_tries(cholesky_max_tries):
for batch_idx in range(n_batches):
start_idx = batch_idx * batch_size_pred
end_idx = min(start_idx + batch_size_pred, len(X_test_t))
X_batch = X_test_t[start_idx:end_idx]
# Get predictive distribution for batch
pred_dist_batch = likelihood_robust(model_robust(X_batch))
# Extract mean and variance (move to CPU, then to numpy)
mean_batch = pred_dist_batch.mean.cpu().numpy()
std_batch = pred_dist_batch.stddev.cpu().numpy()
# Flatten to 1D (handle any batch/feature dimensions)
mean_batch = mean_batch.flatten()
std_batch = std_batch.flatten()
# Sanity check: should match batch size
expected_size = end_idx - start_idx
if mean_batch.size != expected_size:
print(f" WARNING: Batch {batch_idx} size mismatch: got {mean_batch.size}, expected {expected_size}")
# Truncate or pad to expected size
mean_batch = mean_batch[:expected_size]
std_batch = std_batch[:expected_size]
mean_robust_list.append(mean_batch)
std_robust_list.append(std_batch)
if (batch_idx + 1) % 10 == 0 or (batch_idx + 1) == n_batches:
print(f" Processed {end_idx:,}/{len(X_test_t):,} samples ({100*end_idx/len(X_test_t):.1f}%)")
# Concatenate all batches
mean_robust = np.concatenate(mean_robust_list)
std_robust = np.concatenate(std_robust_list)
print(f"✓ Predictions generated")
print(f" Mean range: {mean_robust.min():.2f} → {mean_robust.max():.2f}")
print(f" Std range: {std_robust.min():.2f} → {std_robust.max():.2f}")
In [ ]:
Copied!
# Compute prediction intervals using library function
nu_final = likelihood_robust.deg_free.item()
intervals_robust = compute_prediction_intervals(
mean=mean_robust,
std=std_robust,
confidence_levels=[0.95, 0.99],
distribution="student_t",
nu=nu_final
)
lower_95_robust, upper_95_robust = intervals_robust[0.95]
lower_99_robust, upper_99_robust = intervals_robust[0.99]
# Display quantiles for reference
q_95 = student_t.ppf(0.975, df=nu_final)
q_99 = student_t.ppf(0.995, df=nu_final)
print(f"✓ Prediction intervals computed (Student-t, ν={nu_final:.2f}):")
print(f" 95% quantile: ±{q_95:.3f}")
print(f" 99% quantile: ±{q_99:.3f}")
# Compute prediction intervals using library function
nu_final = likelihood_robust.deg_free.item()
intervals_robust = compute_prediction_intervals(
mean=mean_robust,
std=std_robust,
confidence_levels=[0.95, 0.99],
distribution="student_t",
nu=nu_final
)
lower_95_robust, upper_95_robust = intervals_robust[0.95]
lower_99_robust, upper_99_robust = intervals_robust[0.99]
# Display quantiles for reference
q_95 = student_t.ppf(0.975, df=nu_final)
q_99 = student_t.ppf(0.995, df=nu_final)
print(f"✓ Prediction intervals computed (Student-t, ν={nu_final:.2f}):")
print(f" 95% quantile: ±{q_95:.3f}")
print(f" 99% quantile: ±{q_99:.3f}")
8.2 Traditional Model Predictions¶
In [ ]:
Copied!
# Switch to evaluation mode
model_traditional.eval()
likelihood_traditional.eval()
print("Generating predictions with traditional model...")
print(f"Test samples: {len(X_test_t):,}")
# Batched prediction to prevent memory exhaustion
mean_traditional_list = []
std_traditional_list = []
# Apply same numerical stability settings as training (device-specific)
with torch.no_grad(), \
gpytorch.settings.fast_pred_var(), \
gpytorch.settings.cholesky_jitter(cholesky_jitter), \
gpytorch.settings.cholesky_max_tries(cholesky_max_tries):
for batch_idx in range(n_batches):
start_idx = batch_idx * batch_size_pred
end_idx = min(start_idx + batch_size_pred, len(X_test_t))
X_batch = X_test_t[start_idx:end_idx]
# Get predictive distribution for batch
pred_dist_batch = likelihood_traditional(model_traditional(X_batch))
# Extract mean and variance (move to CPU, then to numpy)
mean_batch = pred_dist_batch.mean.cpu().numpy()
std_batch = pred_dist_batch.stddev.cpu().numpy()
# Flatten to 1D (handle any batch/feature dimensions)
mean_batch = mean_batch.flatten()
std_batch = std_batch.flatten()
# Sanity check: should match batch size
expected_size = end_idx - start_idx
if mean_batch.size != expected_size:
print(f" WARNING: Batch {batch_idx} size mismatch: got {mean_batch.size}, expected {expected_size}")
# Truncate to expected size
mean_batch = mean_batch[:expected_size]
std_batch = std_batch[:expected_size]
mean_traditional_list.append(mean_batch)
std_traditional_list.append(std_batch)
if (batch_idx + 1) % 10 == 0 or (batch_idx + 1) == n_batches:
print(f" Processed {end_idx:,}/{len(X_test_t):,} samples ({100*end_idx/len(X_test_t):.1f}%)")
# Concatenate all batches
mean_traditional = np.concatenate(mean_traditional_list)
std_traditional = np.concatenate(std_traditional_list)
print(f"✓ Predictions generated")
print(f" Mean range: {mean_traditional.min():.2f} → {mean_traditional.max():.2f}")
print(f" Std range: {std_traditional.min():.2f} → {std_traditional.max():.2f}")
# Switch to evaluation mode
model_traditional.eval()
likelihood_traditional.eval()
print("Generating predictions with traditional model...")
print(f"Test samples: {len(X_test_t):,}")
# Batched prediction to prevent memory exhaustion
mean_traditional_list = []
std_traditional_list = []
# Apply same numerical stability settings as training (device-specific)
with torch.no_grad(), \
gpytorch.settings.fast_pred_var(), \
gpytorch.settings.cholesky_jitter(cholesky_jitter), \
gpytorch.settings.cholesky_max_tries(cholesky_max_tries):
for batch_idx in range(n_batches):
start_idx = batch_idx * batch_size_pred
end_idx = min(start_idx + batch_size_pred, len(X_test_t))
X_batch = X_test_t[start_idx:end_idx]
# Get predictive distribution for batch
pred_dist_batch = likelihood_traditional(model_traditional(X_batch))
# Extract mean and variance (move to CPU, then to numpy)
mean_batch = pred_dist_batch.mean.cpu().numpy()
std_batch = pred_dist_batch.stddev.cpu().numpy()
# Flatten to 1D (handle any batch/feature dimensions)
mean_batch = mean_batch.flatten()
std_batch = std_batch.flatten()
# Sanity check: should match batch size
expected_size = end_idx - start_idx
if mean_batch.size != expected_size:
print(f" WARNING: Batch {batch_idx} size mismatch: got {mean_batch.size}, expected {expected_size}")
# Truncate to expected size
mean_batch = mean_batch[:expected_size]
std_batch = std_batch[:expected_size]
mean_traditional_list.append(mean_batch)
std_traditional_list.append(std_batch)
if (batch_idx + 1) % 10 == 0 or (batch_idx + 1) == n_batches:
print(f" Processed {end_idx:,}/{len(X_test_t):,} samples ({100*end_idx/len(X_test_t):.1f}%)")
# Concatenate all batches
mean_traditional = np.concatenate(mean_traditional_list)
std_traditional = np.concatenate(std_traditional_list)
print(f"✓ Predictions generated")
print(f" Mean range: {mean_traditional.min():.2f} → {mean_traditional.max():.2f}")
print(f" Std range: {std_traditional.min():.2f} → {std_traditional.max():.2f}")
In [ ]:
Copied!
# Compute prediction intervals using library function
intervals_traditional = compute_prediction_intervals(
mean=mean_traditional,
std=std_traditional,
confidence_levels=[0.95, 0.99],
distribution="gaussian"
)
lower_95_traditional, upper_95_traditional = intervals_traditional[0.95]
lower_99_traditional, upper_99_traditional = intervals_traditional[0.99]
# Display quantiles for reference
q_95_gauss = norm.ppf(0.975)
q_99_gauss = norm.ppf(0.995)
print(f"✓ Prediction intervals computed (Gaussian):")
print(f" 95% quantile: ±{q_95_gauss:.3f}")
print(f" 99% quantile: ±{q_99_gauss:.3f}")
# Compute prediction intervals using library function
intervals_traditional = compute_prediction_intervals(
mean=mean_traditional,
std=std_traditional,
confidence_levels=[0.95, 0.99],
distribution="gaussian"
)
lower_95_traditional, upper_95_traditional = intervals_traditional[0.95]
lower_99_traditional, upper_99_traditional = intervals_traditional[0.99]
# Display quantiles for reference
q_95_gauss = norm.ppf(0.975)
q_99_gauss = norm.ppf(0.995)
print(f"✓ Prediction intervals computed (Gaussian):")
print(f" 95% quantile: ±{q_95_gauss:.3f}")
print(f" 99% quantile: ±{q_99_gauss:.3f}")
In [ ]:
Copied!
# Compute metrics using library function
metrics_robust = compute_metrics(
y_true=y_test,
y_pred=mean_robust,
lower_95=lower_95_robust,
upper_95=upper_95_robust,
lower_99=lower_99_robust,
upper_99=upper_99_robust,
model_name='Robust (Student-t)'
)
metrics_traditional = compute_metrics(
y_true=y_test,
y_pred=mean_traditional,
lower_95=lower_95_traditional,
upper_95=upper_95_traditional,
lower_99=lower_99_traditional,
upper_99=upper_99_traditional,
model_name='Traditional (Gaussian)'
)
# Create comparison DataFrame
metrics_df = spark.createDataFrame([metrics_robust, metrics_traditional])
metrics_df
# Compute metrics using library function
metrics_robust = compute_metrics(
y_true=y_test,
y_pred=mean_robust,
lower_95=lower_95_robust,
upper_95=upper_95_robust,
lower_99=lower_99_robust,
upper_99=upper_99_robust,
model_name='Robust (Student-t)'
)
metrics_traditional = compute_metrics(
y_true=y_test,
y_pred=mean_traditional,
lower_95=lower_95_traditional,
upper_95=upper_95_traditional,
lower_99=lower_99_traditional,
upper_99=upper_99_traditional,
model_name='Traditional (Gaussian)'
)
# Create comparison DataFrame
metrics_df = spark.createDataFrame([metrics_robust, metrics_traditional])
metrics_df
In [ ]:
Copied!
# Print detailed comparison
print("=" * 70)
print("MODEL EVALUATION METRICS")
print("=" * 70)
for metric_name in ['RMSE', 'MAE', 'R²', 'Coverage 95%', 'Coverage 99%', 'Sharpness 95%', 'Sharpness 99%']:
robust_val = metrics_robust[metric_name]
trad_val = metrics_traditional[metric_name]
# Format based on metric type
if 'Coverage' in metric_name:
print(f"{metric_name:20s} | Robust: {robust_val:6.1%} | Traditional: {trad_val:6.1%}")
else:
print(f"{metric_name:20s} | Robust: {robust_val:6.3f} | Traditional: {trad_val:6.3f}")
print("=" * 70)
# Print detailed comparison
print("=" * 70)
print("MODEL EVALUATION METRICS")
print("=" * 70)
for metric_name in ['RMSE', 'MAE', 'R²', 'Coverage 95%', 'Coverage 99%', 'Sharpness 95%', 'Sharpness 99%']:
robust_val = metrics_robust[metric_name]
trad_val = metrics_traditional[metric_name]
# Format based on metric type
if 'Coverage' in metric_name:
print(f"{metric_name:20s} | Robust: {robust_val:6.1%} | Traditional: {trad_val:6.1%}")
else:
print(f"{metric_name:20s} | Robust: {robust_val:6.3f} | Traditional: {trad_val:6.3f}")
print("=" * 70)
9.2 Detailed Diagnostic Analysis¶
In [ ]:
Copied!
# Run comprehensive diagnostics to understand model behavior
import sys
sys.path.insert(0, '..')
from diagnose_gp_results import diagnose_gp_predictions
# Generate diagnostic report
diagnose_gp_predictions(
y_test=y_test,
mean_robust=mean_robust,
mean_traditional=mean_traditional,
model_robust=model_robust,
model_traditional=model_traditional,
X_test_t=X_test_t,
save_path="../gp_diagnostics.txt"
)
# Display the report
with open("../gp_diagnostics.txt", "r") as f:
print(f.read())
# Run comprehensive diagnostics to understand model behavior
import sys
sys.path.insert(0, '..')
from diagnose_gp_results import diagnose_gp_predictions
# Generate diagnostic report
diagnose_gp_predictions(
y_test=y_test,
mean_robust=mean_robust,
mean_traditional=mean_traditional,
model_robust=model_robust,
model_traditional=model_traditional,
X_test_t=X_test_t,
save_path="../gp_diagnostics.txt"
)
# Display the report
with open("../gp_diagnostics.txt", "r") as f:
print(f.read())
9.3 Calibration Analysis¶
In [ ]:
Copied!
# Compute standardized residuals
residuals_robust = (y_test - mean_robust) / std_robust
residuals_traditional = (y_test - mean_traditional) / std_traditional
# Expected vs observed quantiles (Q-Q plot)
quantiles = np.linspace(0.01, 0.99, 99)
# For robust model, compare against Student-t
expected_quantiles_robust = student_t.ppf(quantiles, df=nu_final)
observed_quantiles_robust = np.quantile(residuals_robust, quantiles)
# For traditional model, compare against Gaussian
expected_quantiles_trad = norm.ppf(quantiles)
observed_quantiles_trad = np.quantile(residuals_traditional, quantiles)
# Plot calibration
fig, axes = plt.subplots(1, 2, figsize=(16, 6))
# Robust model calibration
axes[0].scatter(expected_quantiles_robust, observed_quantiles_robust, alpha=0.6, s=30, color='steelblue')
axes[0].plot([-4, 4], [-4, 4], 'r--', linewidth=2, label='Perfect Calibration')
axes[0].set_xlabel(f'Expected Quantiles (Student-t, ν={nu_final:.2f})', fontsize=12)
axes[0].set_ylabel('Observed Quantiles (Standardized Residuals)', fontsize=12)
axes[0].set_title('Robust Model Calibration', fontsize=14, fontweight='bold')
axes[0].legend()
axes[0].grid(alpha=0.3)
# Traditional model calibration
axes[1].scatter(expected_quantiles_trad, observed_quantiles_trad, alpha=0.6, s=30, color='coral')
axes[1].plot([-4, 4], [-4, 4], 'r--', linewidth=2, label='Perfect Calibration')
axes[1].set_xlabel('Expected Quantiles (Gaussian)', fontsize=12)
axes[1].set_ylabel('Observed Quantiles (Standardized Residuals)', fontsize=12)
axes[1].set_title('Traditional Model Calibration', fontsize=14, fontweight='bold')
axes[1].legend()
axes[1].grid(alpha=0.3)
plt.tight_layout()
plt.show()
# Compute standardized residuals
residuals_robust = (y_test - mean_robust) / std_robust
residuals_traditional = (y_test - mean_traditional) / std_traditional
# Expected vs observed quantiles (Q-Q plot)
quantiles = np.linspace(0.01, 0.99, 99)
# For robust model, compare against Student-t
expected_quantiles_robust = student_t.ppf(quantiles, df=nu_final)
observed_quantiles_robust = np.quantile(residuals_robust, quantiles)
# For traditional model, compare against Gaussian
expected_quantiles_trad = norm.ppf(quantiles)
observed_quantiles_trad = np.quantile(residuals_traditional, quantiles)
# Plot calibration
fig, axes = plt.subplots(1, 2, figsize=(16, 6))
# Robust model calibration
axes[0].scatter(expected_quantiles_robust, observed_quantiles_robust, alpha=0.6, s=30, color='steelblue')
axes[0].plot([-4, 4], [-4, 4], 'r--', linewidth=2, label='Perfect Calibration')
axes[0].set_xlabel(f'Expected Quantiles (Student-t, ν={nu_final:.2f})', fontsize=12)
axes[0].set_ylabel('Observed Quantiles (Standardized Residuals)', fontsize=12)
axes[0].set_title('Robust Model Calibration', fontsize=14, fontweight='bold')
axes[0].legend()
axes[0].grid(alpha=0.3)
# Traditional model calibration
axes[1].scatter(expected_quantiles_trad, observed_quantiles_trad, alpha=0.6, s=30, color='coral')
axes[1].plot([-4, 4], [-4, 4], 'r--', linewidth=2, label='Perfect Calibration')
axes[1].set_xlabel('Expected Quantiles (Gaussian)', fontsize=12)
axes[1].set_ylabel('Observed Quantiles (Standardized Residuals)', fontsize=12)
axes[1].set_title('Traditional Model Calibration', fontsize=14, fontweight='bold')
axes[1].legend()
axes[1].grid(alpha=0.3)
plt.tight_layout()
plt.show()
In [ ]:
Copied!
# Flag anomalies using prediction intervals
anomalies_95_robust = (y_test < lower_95_robust) | (y_test > upper_95_robust)
anomalies_99_robust = (y_test < lower_99_robust) | (y_test > upper_99_robust)
anomalies_95_traditional = (y_test < lower_95_traditional) | (y_test > upper_95_traditional)
anomalies_99_traditional = (y_test < lower_99_traditional) | (y_test > upper_99_traditional)
print(f"Anomalies detected:")
print(f" Robust (95%): {anomalies_95_robust.sum():,}")
print(f" Robust (99%): {anomalies_99_robust.sum():,}")
print(f" Traditional (95%): {anomalies_95_traditional.sum():,}")
print(f" Traditional (99%): {anomalies_99_traditional.sum():,}")
print(f" Ground truth: {anomaly_test.sum():,}")
# Flag anomalies using prediction intervals
anomalies_95_robust = (y_test < lower_95_robust) | (y_test > upper_95_robust)
anomalies_99_robust = (y_test < lower_99_robust) | (y_test > upper_99_robust)
anomalies_95_traditional = (y_test < lower_95_traditional) | (y_test > upper_95_traditional)
anomalies_99_traditional = (y_test < lower_99_traditional) | (y_test > upper_99_traditional)
print(f"Anomalies detected:")
print(f" Robust (95%): {anomalies_95_robust.sum():,}")
print(f" Robust (99%): {anomalies_99_robust.sum():,}")
print(f" Traditional (95%): {anomalies_95_traditional.sum():,}")
print(f" Traditional (99%): {anomalies_99_traditional.sum():,}")
print(f" Ground truth: {anomaly_test.sum():,}")
10.2 Anomaly Detection Metrics¶
In [ ]:
Copied!
# Compute anomaly detection metrics using library function
anomaly_metrics = [
compute_anomaly_metrics(
y_true_anomaly=anomaly_test,
y_pred_anomaly=anomalies_95_robust,
model_name='Robust (Student-t)',
threshold_name='95% Interval'
),
compute_anomaly_metrics(
y_true_anomaly=anomaly_test,
y_pred_anomaly=anomalies_99_robust,
model_name='Robust (Student-t)',
threshold_name='99% Interval'
),
compute_anomaly_metrics(
y_true_anomaly=anomaly_test,
y_pred_anomaly=anomalies_95_traditional,
model_name='Traditional (Gaussian)',
threshold_name='95% Interval'
),
compute_anomaly_metrics(
y_true_anomaly=anomaly_test,
y_pred_anomaly=anomalies_99_traditional,
model_name='Traditional (Gaussian)',
threshold_name='99% Interval'
),
]
anomaly_metrics_df = spark.createDataFrame(anomaly_metrics)
anomaly_metrics_df
# Compute anomaly detection metrics using library function
anomaly_metrics = [
compute_anomaly_metrics(
y_true_anomaly=anomaly_test,
y_pred_anomaly=anomalies_95_robust,
model_name='Robust (Student-t)',
threshold_name='95% Interval'
),
compute_anomaly_metrics(
y_true_anomaly=anomaly_test,
y_pred_anomaly=anomalies_99_robust,
model_name='Robust (Student-t)',
threshold_name='99% Interval'
),
compute_anomaly_metrics(
y_true_anomaly=anomaly_test,
y_pred_anomaly=anomalies_95_traditional,
model_name='Traditional (Gaussian)',
threshold_name='95% Interval'
),
compute_anomaly_metrics(
y_true_anomaly=anomaly_test,
y_pred_anomaly=anomalies_99_traditional,
model_name='Traditional (Gaussian)',
threshold_name='99% Interval'
),
]
anomaly_metrics_df = spark.createDataFrame(anomaly_metrics)
anomaly_metrics_df
10.3 ROC Curve Analysis¶
In [ ]:
Copied!
# Compute anomaly scores (standardized residuals)
anomaly_scores_robust = np.abs((y_test - mean_robust) / std_robust)
anomaly_scores_traditional = np.abs((y_test - mean_traditional) / std_traditional)
# Compute ROC curves
fpr_robust, tpr_robust, thresholds_robust = roc_curve(anomaly_test, anomaly_scores_robust)
fpr_traditional, tpr_traditional, thresholds_traditional = roc_curve(anomaly_test, anomaly_scores_traditional)
# Compute AUC-ROC
auc_robust = roc_auc_score(anomaly_test, anomaly_scores_robust)
auc_traditional = roc_auc_score(anomaly_test, anomaly_scores_traditional)
print(f"AUC-ROC Scores:")
print(f" Robust: {auc_robust:.4f}")
print(f" Traditional: {auc_traditional:.4f}")
# Compute anomaly scores (standardized residuals)
anomaly_scores_robust = np.abs((y_test - mean_robust) / std_robust)
anomaly_scores_traditional = np.abs((y_test - mean_traditional) / std_traditional)
# Compute ROC curves
fpr_robust, tpr_robust, thresholds_robust = roc_curve(anomaly_test, anomaly_scores_robust)
fpr_traditional, tpr_traditional, thresholds_traditional = roc_curve(anomaly_test, anomaly_scores_traditional)
# Compute AUC-ROC
auc_robust = roc_auc_score(anomaly_test, anomaly_scores_robust)
auc_traditional = roc_auc_score(anomaly_test, anomaly_scores_traditional)
print(f"AUC-ROC Scores:")
print(f" Robust: {auc_robust:.4f}")
print(f" Traditional: {auc_traditional:.4f}")
In [ ]:
Copied!
# Plot ROC curves
fig, ax = plt.subplots(figsize=(10, 8))
ax.plot(fpr_robust, tpr_robust, linewidth=2.5,
label=f'Robust (Student-t) - AUC = {auc_robust:.3f}', color='steelblue')
ax.plot(fpr_traditional, tpr_traditional, linewidth=2.5,
label=f'Traditional (Gaussian) - AUC = {auc_traditional:.3f}', color='coral')
ax.plot([0, 1], [0, 1], 'k--', linewidth=1.5, label='Random Classifier', alpha=0.5)
ax.set_xlabel('False Positive Rate', fontsize=13)
ax.set_ylabel('True Positive Rate', fontsize=13)
ax.set_title('ROC Curve: Anomaly Detection Performance', fontsize=15, fontweight='bold')
ax.legend(fontsize=12, loc='lower right')
ax.grid(alpha=0.3)
plt.tight_layout()
plt.show()
# Plot ROC curves
fig, ax = plt.subplots(figsize=(10, 8))
ax.plot(fpr_robust, tpr_robust, linewidth=2.5,
label=f'Robust (Student-t) - AUC = {auc_robust:.3f}', color='steelblue')
ax.plot(fpr_traditional, tpr_traditional, linewidth=2.5,
label=f'Traditional (Gaussian) - AUC = {auc_traditional:.3f}', color='coral')
ax.plot([0, 1], [0, 1], 'k--', linewidth=1.5, label='Random Classifier', alpha=0.5)
ax.set_xlabel('False Positive Rate', fontsize=13)
ax.set_ylabel('True Positive Rate', fontsize=13)
ax.set_title('ROC Curve: Anomaly Detection Performance', fontsize=15, fontweight='bold')
ax.legend(fontsize=12, loc='lower right')
ax.grid(alpha=0.3)
plt.tight_layout()
plt.show()
In [ ]:
Copied!
# Plot subset for clarity (first 5000 test points)
plot_start = 0
plot_end = 5000
X_plot = X_test[plot_start:plot_end].ravel()
y_plot = y_test[plot_start:plot_end]
mean_plot = mean_robust[plot_start:plot_end]
lower_95_plot = lower_95_robust[plot_start:plot_end]
upper_95_plot = upper_95_robust[plot_start:plot_end]
anomaly_plot = anomaly_test[plot_start:plot_end]
fig, ax = plt.subplots(figsize=(18, 6))
# Observations
ax.plot(X_plot, y_plot, 'k.', alpha=0.3, markersize=2, label='Observed', zorder=1)
# GP mean
ax.plot(X_plot, mean_plot, 'b-', linewidth=2, label='GP Mean (Robust)', zorder=3)
# 95% prediction interval
ax.fill_between(
X_plot,
lower_95_plot,
upper_95_plot,
alpha=0.25,
color='steelblue',
label='95% Prediction Interval',
zorder=2
)
# Highlight labeled anomalies
anomaly_mask = anomaly_plot.astype(bool)
ax.scatter(
X_plot[anomaly_mask],
y_plot[anomaly_mask],
color='red',
s=60,
marker='x',
linewidths=2.5,
label=f'Labeled Anomalies ({anomaly_mask.sum()})',
zorder=4
)
ax.set_xlabel('Timestamp', fontsize=12)
ax.set_ylabel('KPI Value (IOPS)', fontsize=12)
ax.set_title('Robust GP: Predictions with Uncertainty Quantification', fontsize=14, fontweight='bold')
ax.legend(loc='upper right', fontsize=11)
ax.grid(alpha=0.3)
plt.tight_layout()
plt.show()
# Plot subset for clarity (first 5000 test points)
plot_start = 0
plot_end = 5000
X_plot = X_test[plot_start:plot_end].ravel()
y_plot = y_test[plot_start:plot_end]
mean_plot = mean_robust[plot_start:plot_end]
lower_95_plot = lower_95_robust[plot_start:plot_end]
upper_95_plot = upper_95_robust[plot_start:plot_end]
anomaly_plot = anomaly_test[plot_start:plot_end]
fig, ax = plt.subplots(figsize=(18, 6))
# Observations
ax.plot(X_plot, y_plot, 'k.', alpha=0.3, markersize=2, label='Observed', zorder=1)
# GP mean
ax.plot(X_plot, mean_plot, 'b-', linewidth=2, label='GP Mean (Robust)', zorder=3)
# 95% prediction interval
ax.fill_between(
X_plot,
lower_95_plot,
upper_95_plot,
alpha=0.25,
color='steelblue',
label='95% Prediction Interval',
zorder=2
)
# Highlight labeled anomalies
anomaly_mask = anomaly_plot.astype(bool)
ax.scatter(
X_plot[anomaly_mask],
y_plot[anomaly_mask],
color='red',
s=60,
marker='x',
linewidths=2.5,
label=f'Labeled Anomalies ({anomaly_mask.sum()})',
zorder=4
)
ax.set_xlabel('Timestamp', fontsize=12)
ax.set_ylabel('KPI Value (IOPS)', fontsize=12)
ax.set_title('Robust GP: Predictions with Uncertainty Quantification', fontsize=14, fontweight='bold')
ax.legend(loc='upper right', fontsize=11)
ax.grid(alpha=0.3)
plt.tight_layout()
plt.show()
11.2 Pattern Reconstruction: Sawtooth × Sinusoidal¶
In [ ]:
Copied!
# Zoom into one complete sawtooth cycle (~1250 normalized steps)
cycle_length_norm = 1250 / X_range
start_idx = 1000
end_idx = start_idx + int(cycle_length_norm * len(X_test))
# Ensure we don't exceed bounds
end_idx = min(end_idx, len(X_test))
X_cycle = X_test[start_idx:end_idx].ravel()
y_cycle = y_test[start_idx:end_idx]
mean_cycle = mean_robust[start_idx:end_idx]
lower_cycle = lower_95_robust[start_idx:end_idx]
upper_cycle = upper_95_robust[start_idx:end_idx]
fig, ax = plt.subplots(figsize=(18, 6))
# Observations
ax.plot(X_cycle, y_cycle, 'k.', alpha=0.4, markersize=3, label='Observed')
# GP mean (should capture sawtooth × sine)
ax.plot(X_cycle, mean_cycle, 'b-', linewidth=2.5, label='GP Mean (captures pattern)')
# Prediction interval
ax.fill_between(X_cycle, lower_cycle, upper_cycle, alpha=0.25, color='steelblue', label='95% Interval')
# Annotate the 5 sinusoidal oscillations
if len(X_cycle) > 0:
cycle_width = (X_cycle[-1] - X_cycle[0]) / 5
for i in range(5):
segment_start = X_cycle[0] + i * cycle_width
segment_end = segment_start + cycle_width
ax.axvspan(segment_start, segment_end, alpha=0.08,
color=['orange', 'green', 'blue', 'purple', 'red'][i])
ax.set_xlabel('Timestamp', fontsize=12)
ax.set_ylabel('KPI Value (IOPS)', fontsize=12)
ax.set_title('Pattern Reconstruction: Sawtooth Envelope × 5 Sinusoidal Oscillations',
fontsize=14, fontweight='bold')
ax.legend(loc='upper left', fontsize=11)
ax.grid(alpha=0.3)
plt.tight_layout()
plt.show()
# Zoom into one complete sawtooth cycle (~1250 normalized steps)
cycle_length_norm = 1250 / X_range
start_idx = 1000
end_idx = start_idx + int(cycle_length_norm * len(X_test))
# Ensure we don't exceed bounds
end_idx = min(end_idx, len(X_test))
X_cycle = X_test[start_idx:end_idx].ravel()
y_cycle = y_test[start_idx:end_idx]
mean_cycle = mean_robust[start_idx:end_idx]
lower_cycle = lower_95_robust[start_idx:end_idx]
upper_cycle = upper_95_robust[start_idx:end_idx]
fig, ax = plt.subplots(figsize=(18, 6))
# Observations
ax.plot(X_cycle, y_cycle, 'k.', alpha=0.4, markersize=3, label='Observed')
# GP mean (should capture sawtooth × sine)
ax.plot(X_cycle, mean_cycle, 'b-', linewidth=2.5, label='GP Mean (captures pattern)')
# Prediction interval
ax.fill_between(X_cycle, lower_cycle, upper_cycle, alpha=0.25, color='steelblue', label='95% Interval')
# Annotate the 5 sinusoidal oscillations
if len(X_cycle) > 0:
cycle_width = (X_cycle[-1] - X_cycle[0]) / 5
for i in range(5):
segment_start = X_cycle[0] + i * cycle_width
segment_end = segment_start + cycle_width
ax.axvspan(segment_start, segment_end, alpha=0.08,
color=['orange', 'green', 'blue', 'purple', 'red'][i])
ax.set_xlabel('Timestamp', fontsize=12)
ax.set_ylabel('KPI Value (IOPS)', fontsize=12)
ax.set_title('Pattern Reconstruction: Sawtooth Envelope × 5 Sinusoidal Oscillations',
fontsize=14, fontweight='bold')
ax.legend(loc='upper left', fontsize=11)
ax.grid(alpha=0.3)
plt.tight_layout()
plt.show()
11.3 Residual Analysis¶
In [ ]:
Copied!
# Residual plots for both models
fig, axes = plt.subplots(2, 2, figsize=(16, 10))
# Robust model residuals
axes[0, 0].scatter(mean_robust, residuals_robust, alpha=0.3, s=10, color='steelblue')
axes[0, 0].axhline(y=0, color='red', linestyle='--', linewidth=2)
axes[0, 0].set_xlabel('Predicted Value', fontsize=11)
axes[0, 0].set_ylabel('Standardized Residuals', fontsize=11)
axes[0, 0].set_title('Robust: Residuals vs Predicted', fontsize=12, fontweight='bold')
axes[0, 0].grid(alpha=0.3)
axes[0, 1].hist(residuals_robust, bins=100, alpha=0.7, color='steelblue', edgecolor='black')
axes[0, 1].axvline(x=0, color='red', linestyle='--', linewidth=2)
axes[0, 1].set_xlabel('Standardized Residuals', fontsize=11)
axes[0, 1].set_ylabel('Frequency', fontsize=11)
axes[0, 1].set_title('Robust: Residual Distribution', fontsize=12, fontweight='bold')
axes[0, 1].grid(alpha=0.3)
# Traditional model residuals
axes[1, 0].scatter(mean_traditional, residuals_traditional, alpha=0.3, s=10, color='coral')
axes[1, 0].axhline(y=0, color='red', linestyle='--', linewidth=2)
axes[1, 0].set_xlabel('Predicted Value', fontsize=11)
axes[1, 0].set_ylabel('Standardized Residuals', fontsize=11)
axes[1, 0].set_title('Traditional: Residuals vs Predicted', fontsize=12, fontweight='bold')
axes[1, 0].grid(alpha=0.3)
axes[1, 1].hist(residuals_traditional, bins=100, alpha=0.7, color='coral', edgecolor='black')
axes[1, 1].axvline(x=0, color='red', linestyle='--', linewidth=2)
axes[1, 1].set_xlabel('Standardized Residuals', fontsize=11)
axes[1, 1].set_ylabel('Frequency', fontsize=11)
axes[1, 1].set_title('Traditional: Residual Distribution', fontsize=12, fontweight='bold')
axes[1, 1].grid(alpha=0.3)
plt.tight_layout()
plt.show()
# Residual plots for both models
fig, axes = plt.subplots(2, 2, figsize=(16, 10))
# Robust model residuals
axes[0, 0].scatter(mean_robust, residuals_robust, alpha=0.3, s=10, color='steelblue')
axes[0, 0].axhline(y=0, color='red', linestyle='--', linewidth=2)
axes[0, 0].set_xlabel('Predicted Value', fontsize=11)
axes[0, 0].set_ylabel('Standardized Residuals', fontsize=11)
axes[0, 0].set_title('Robust: Residuals vs Predicted', fontsize=12, fontweight='bold')
axes[0, 0].grid(alpha=0.3)
axes[0, 1].hist(residuals_robust, bins=100, alpha=0.7, color='steelblue', edgecolor='black')
axes[0, 1].axvline(x=0, color='red', linestyle='--', linewidth=2)
axes[0, 1].set_xlabel('Standardized Residuals', fontsize=11)
axes[0, 1].set_ylabel('Frequency', fontsize=11)
axes[0, 1].set_title('Robust: Residual Distribution', fontsize=12, fontweight='bold')
axes[0, 1].grid(alpha=0.3)
# Traditional model residuals
axes[1, 0].scatter(mean_traditional, residuals_traditional, alpha=0.3, s=10, color='coral')
axes[1, 0].axhline(y=0, color='red', linestyle='--', linewidth=2)
axes[1, 0].set_xlabel('Predicted Value', fontsize=11)
axes[1, 0].set_ylabel('Standardized Residuals', fontsize=11)
axes[1, 0].set_title('Traditional: Residuals vs Predicted', fontsize=12, fontweight='bold')
axes[1, 0].grid(alpha=0.3)
axes[1, 1].hist(residuals_traditional, bins=100, alpha=0.7, color='coral', edgecolor='black')
axes[1, 1].axvline(x=0, color='red', linestyle='--', linewidth=2)
axes[1, 1].set_xlabel('Standardized Residuals', fontsize=11)
axes[1, 1].set_ylabel('Frequency', fontsize=11)
axes[1, 1].set_title('Traditional: Residual Distribution', fontsize=12, fontweight='bold')
axes[1, 1].grid(alpha=0.3)
plt.tight_layout()
plt.show()
12. Summary and Recommendations¶
In [ ]:
Copied!
# Create comprehensive comparison table
comparison_data = {
'Metric': [
'Training Samples',
'Likelihood',
'Degrees of Freedom (ν)',
'---',
'RMSE',
'MAE',
'R²',
'---',
'95% Coverage',
'99% Coverage',
'95% Sharpness',
'99% Sharpness',
'---',
'Precision (95% interval)',
'Recall (95% interval)',
'F1-Score (95% interval)',
'AUC-ROC',
],
'Robust (Student-t)': [
f"{len(X_train_t):,} (all data)",
'Student-t',
f'{nu_final:.2f}',
'---',
f'{metrics_robust["RMSE"]:.3f}',
f'{metrics_robust["MAE"]:.3f}',
f'{metrics_robust["R²"]:.3f}',
'---',
f'{metrics_robust["Coverage 95%"]:.1%}',
f'{metrics_robust["Coverage 99%"]:.1%}',
f'{metrics_robust["Sharpness 95%"]:.3f}',
f'{metrics_robust["Sharpness 99%"]:.3f}',
'---',
f'{anomaly_metrics[0]["Precision"]:.3f}',
f'{anomaly_metrics[0]["Recall"]:.3f}',
f'{anomaly_metrics[0]["F1-Score"]:.3f}',
f'{auc_robust:.3f}',
],
'Traditional (Gaussian)': [
f"{len(X_train_clean_t):,} (exclude anomalies)",
'Gaussian',
'∞ (normal)',
'---',
f'{metrics_traditional["RMSE"]:.3f}',
f'{metrics_traditional["MAE"]:.3f}',
f'{metrics_traditional["R²"]:.3f}',
'---',
f'{metrics_traditional["Coverage 95%"]:.1%}',
f'{metrics_traditional["Coverage 99%"]:.1%}',
f'{metrics_traditional["Sharpness 95%"]:.3f}',
f'{metrics_traditional["Sharpness 99%"]:.3f}',
'---',
f'{anomaly_metrics[2]["Precision"]:.3f}',
f'{anomaly_metrics[2]["Recall"]:.3f}',
f'{anomaly_metrics[2]["F1-Score"]:.3f}',
f'{auc_traditional:.3f}',
]
}
comparison_df = spark.createDataFrame(comparison_data)
print("=" * 80)
print("COMPREHENSIVE MODEL COMPARISON")
print("=" * 80)
print(comparison_df.toPandas().to_string(index=False))
print("=" * 80)
# Create comprehensive comparison table
comparison_data = {
'Metric': [
'Training Samples',
'Likelihood',
'Degrees of Freedom (ν)',
'---',
'RMSE',
'MAE',
'R²',
'---',
'95% Coverage',
'99% Coverage',
'95% Sharpness',
'99% Sharpness',
'---',
'Precision (95% interval)',
'Recall (95% interval)',
'F1-Score (95% interval)',
'AUC-ROC',
],
'Robust (Student-t)': [
f"{len(X_train_t):,} (all data)",
'Student-t',
f'{nu_final:.2f}',
'---',
f'{metrics_robust["RMSE"]:.3f}',
f'{metrics_robust["MAE"]:.3f}',
f'{metrics_robust["R²"]:.3f}',
'---',
f'{metrics_robust["Coverage 95%"]:.1%}',
f'{metrics_robust["Coverage 99%"]:.1%}',
f'{metrics_robust["Sharpness 95%"]:.3f}',
f'{metrics_robust["Sharpness 99%"]:.3f}',
'---',
f'{anomaly_metrics[0]["Precision"]:.3f}',
f'{anomaly_metrics[0]["Recall"]:.3f}',
f'{anomaly_metrics[0]["F1-Score"]:.3f}',
f'{auc_robust:.3f}',
],
'Traditional (Gaussian)': [
f"{len(X_train_clean_t):,} (exclude anomalies)",
'Gaussian',
'∞ (normal)',
'---',
f'{metrics_traditional["RMSE"]:.3f}',
f'{metrics_traditional["MAE"]:.3f}',
f'{metrics_traditional["R²"]:.3f}',
'---',
f'{metrics_traditional["Coverage 95%"]:.1%}',
f'{metrics_traditional["Coverage 99%"]:.1%}',
f'{metrics_traditional["Sharpness 95%"]:.3f}',
f'{metrics_traditional["Sharpness 99%"]:.3f}',
'---',
f'{anomaly_metrics[2]["Precision"]:.3f}',
f'{anomaly_metrics[2]["Recall"]:.3f}',
f'{anomaly_metrics[2]["F1-Score"]:.3f}',
f'{auc_traditional:.3f}',
]
}
comparison_df = spark.createDataFrame(comparison_data)
print("=" * 80)
print("COMPREHENSIVE MODEL COMPARISON")
print("=" * 80)
print(comparison_df.toPandas().to_string(index=False))
print("=" * 80)
Key Findings¶
1. Pattern Reconstruction:
- ✅ Composite periodic kernel successfully captures sawtooth × sinusoidal pattern
- ✅ Visual inspection confirms 5 oscillations per sawtooth cycle
- ✅ GP mean tracks complex two-scale structure
2. Calibration Performance:
- ✅ Robust model: Better calibrated prediction intervals (closer to nominal 95%/99%)
- ⚠ Traditional model: May under/overestimate uncertainty due to Gaussian assumption
3. Anomaly Detection:
- ✅ Robust model: Higher recall (catches more true anomalies)
- ✅ AUC-ROC: Both models show strong discriminative ability
- ⚠ Trade-off: Precision vs Recall depends on interval threshold (95% vs 99%)
4. Production Readiness:
- ✅ Robust approach: Trained on real operational data (all samples)
- ✅ Student-t likelihood: Naturally handles outliers without pre-processing
- ✅ Sparse GP: Scalable to large datasets (O(nm²) complexity)
- ✅ GPyTorch: Production-proven library (Uber, Meta, Amazon)
Recommendations¶
For Production Deployment:
- Use robust approach (Student-t likelihood, train on all data)
- Adaptive thresholding: Tune 95%/99% intervals based on operational cost of false positives/negatives
- Online learning: Incrementally update model with new data
- Ensemble methods: Combine GP with other anomaly detectors for robustness
For Further Improvement: 5. Add covariates: Time-of-day, day-of-week features 6. Multi-output GP: Model multiple related KPIs jointly 7. Spectral mixture kernel: Automate period discovery 8. Deeper inducing points: Increase M for higher accuracy (trade-off: computational cost)
Next Steps:
- Deploy model via FastAPI endpoint
- Integrate with monitoring dashboard
- A/B test against existing anomaly detection system
- Collect feedback from operations team
References¶
- GPyTorch: https://gpytorch.ai/
- Sparse GPs (SVGP): Hensman et al. (2013), "Gaussian Processes for Big Data"
- Student-t Processes: Shah et al. (2014), "Student-t Processes as Alternatives to Gaussian Processes"
- Composite Kernels: Rasmussen & Williams (2006), "Gaussian Processes for Machine Learning"
In [ ]:
Copied!
In [ ]:
Copied!