Prediction vs Confidence Interval

import matplotlib.pyplot as plt
import numpy as np

np.random.seed(42)
x = np.linspace(0, 10, 50)
y_true = 5 + 2*x
y = y_true + np.random.normal(0, 2, len(x))

# Fit
coef = np.polyfit(x, y, 1)
y_fit = np.polyval(coef, x)
residuals = y - y_fit
mse = np.mean(residuals**2)
se_fit = np.sqrt(mse * (1/len(x) + (x - np.mean(x))**2 / np.sum((x - np.mean(x))**2)))
se_pred = np.sqrt(mse * (1 + 1/len(x) + (x - np.mean(x))**2 / np.sum((x - np.mean(x))**2)))

fig, ax = plt.subplots(figsize=(10, 6), facecolor='#ffffff')
ax.set_facecolor('#ffffff')

# Prediction interval (wider)
ax.fill_between(x, y_fit - 1.96*se_pred, y_fit + 1.96*se_pred, alpha=0.15, color='#F5276C', label='95% Prediction')
# Confidence interval (narrower)
ax.fill_between(x, y_fit - 1.96*se_fit, y_fit + 1.96*se_fit, alpha=0.35, color='#27D3F5', label='95% Confidence')
ax.scatter(x, y, color='#374151', s=40, alpha=0.7, edgecolors='white', linewidth=0.5)
ax.plot(x, y_fit, color='#F5276C', linewidth=2.5, label='Fitted line')

ax.text(0.05, 0.95, f'y = {coef[1]:.2f} + {coef[0]:.2f}x\nR² = {1 - np.var(residuals)/np.var(y):.3f}', 
        transform=ax.transAxes, color='#1f2937', fontsize=10, va='top',
        bbox=dict(boxstyle='round', facecolor='#f3f4f6', edgecolor='#e5e7eb'))

ax.set_xlabel('x', color='#1f2937', fontsize=11, fontweight='500')
ax.set_ylabel('y', color='#1f2937', fontsize=11, fontweight='500')
ax.set_title('Confidence vs Prediction Intervals', color='#1f2937', fontsize=13, fontweight='600', pad=15)
ax.tick_params(colors='#374151', labelsize=9)
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
ax.spines['left'].set_color('#e5e7eb')
ax.spines['bottom'].set_color('#e5e7eb')
ax.yaxis.grid(True, color='#f3f4f6', linewidth=0.8)
ax.legend(facecolor='#ffffff', edgecolor='#e5e7eb', labelcolor='#1f2937', loc='lower right')

plt.tight_layout()
plt.show()