Comparing scatter_smooth with R’s smooth.spline

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import time
from scatter_smooth.datasets import load_bikeshare
from scatter_smooth import SplineSmoother

%load_ext rpy2.ipython

Bike = pd.read_csv(load_bikeshare())
Bike.head()

df = 7
# 'bikers' is 'cnt' in the original dataset, ISLP might have renamed it or we use 'cnt'
if 'bikers' not in Bike.columns:
    Bike['bikers'] = Bike['cnt']

hr_numeric = pd.to_numeric(Bike['hr'])
bikers = Bike['bikers']

# Sort by hour for cleaner plotting lines
sorted_idx = np.argsort(hr_numeric)
x_plot = hr_numeric.iloc[sorted_idx].unique()

# Fit model
spl_py = SplineSmoother(x=hr_numeric, df=df)
spl_py.smooth(bikers)

# Predict
y_py = spl_py.predict(x_plot)

%%R -i hr_numeric -i bikers -o y_r -i df
# Fit model in R
fit_r <- smooth.spline(hr_numeric, bikers, df=df)

# Predict at unique hours
# unique() in R returns unsorted, but we want to match x_plot order
x_vals <- sort(unique(hr_numeric))
pred_r <- predict(fit_r, x_vals)
y_r <- pred_r$y

fig, ax = plt.subplots(figsize=(10, 6))
ax.scatter(hr_numeric + np.random.normal(0, 0.1, len(hr_numeric)), bikers, 
            s=1, c='lightgray', alpha=0.5, label='Data')
ax.plot(x_plot, y_py, 'b-', lw=3, label='Python (scatter_smooth)', alpha=0.8)
ax.plot(x_plot, y_r, 'r--', lw=3, label='R (smooth.spline)', alpha=0.8)
ax.set_xlabel("Hour")
ax.set_ylabel("Number of Bikers")
ax.set_title("Comparison of Smoothing Splines (df={df})")
ax.legend()
plt.show()

# Numerical comparison
# Note: R might handle repeated x values slightly differently (using weights)
# scatter_smooth handles them naturally in the basis construction.
diff = np.mean(np.abs(y_py - y_r))
print(f"Mean Absolute Difference: {diff:.6f}")

# Python Timing
t_py = %timeit -o -n 10 -r 3 SplineSmoother(x=hr_numeric, df=10).smooth(bikers)

# R Timing

%%R -i hr_numeric -i bikers
if (!require("microbenchmark", quietly = TRUE)) {
    install.packages("microbenchmark", repos="http://cloud.r-project.org")
}
library(microbenchmark)
# R Timing
summary(microbenchmark(
  smooth.spline(hr_numeric, bikers, df=10),
  times=10), unit='milliseconds')[,-1]

%%R
fit_gcv <- smooth.spline(hr_numeric, bikers, cv=FALSE) # cv=FALSE implies GCV
cat("Selected df (R):", fit_gcv$df, "
")
cat("Selected lambda (R):", fit_gcv$lambda, "
")

# Initialize fitter with data
# Note: We use the internal C++ fitter for speed if available
fitter = SplineSmoother(x=hr_numeric, knots=np.unique(hr_numeric))

# Solve for GCV
best_lam = fitter.solve_gcv(bikers)
print(f"Selected lambda (Python): {best_lam}")

# Get corresponding df
# We can access the internal fitter to compute DF for verification
if fitter._cpp_fitter:
    best_lam_scaled = best_lam / fitter.x_scale_**3
    best_df = fitter._cpp_fitter.compute_df(best_lam_scaled)
    print(f"Selected df (Python): {best_df}")
else:
    best_df = "N/A (C++ extension not available)"

fig, ax = plt.subplots(figsize=(10, 6))
ax.scatter(hr_numeric + np.random.normal(0, 0.1, len(hr_numeric)), bikers, 
            s=1, c='lightgray')
# fitter is already fitted with best_lam by solve_gcv
ax.plot(x_plot, fitter.predict(x_plot), 'g-', lw=3, label=f'Optimal GCV (df={best_df:.2f})')
ax.set_xlabel("Hour")
ax.set_ylabel("Number of Bikers")
ax.legend()
plt.show()

# Fit model on log(bikers)
log_bikers = np.log(bikers + 1) # Add 1 to avoid log(0)
spl_log = SplineSmoother(x=hr_numeric, df=df)
spl_log.smooth(log_bikers)

fig, ax = plt.subplots(figsize=(10, 6))
ax.scatter(hr_numeric + np.random.normal(0, 0.1, len(hr_numeric)), log_bikers, 
            s=1, c='lightgray')
ax.plot(x_plot, spl_log.predict(x_plot), 'purple', lw=3, label='Log-Smoothing Spline (df={df})')
ax.set_xlabel("Hour")
ax.set_ylabel("Log(Number of Bikers + 1)")
ax.legend()
plt.show()

# Synthetic data
n_syn = 500
df = 20
rng = np.random.default_rng(42)
x_syn = np.sort(rng.uniform(0, 10, n_syn))
y_syn = np.sin(np.sqrt(x_syn)*2) + rng.normal(0, 1, n_syn)
x_plot = np.linspace(x_syn.min() - 1, x_syn.max() + 1, 200)

# Python Fitting
spline = SplineSmoother(x=x_syn, df=df)
spline.smooth(y_syn)
y_plot_py = spline.predict(x_plot)

# Python Timing (using all unique x as knots)
print(f"Python Timing (n={n_syn}, all knots):")
%timeit SplineSmoother(x=x_syn, df=df).smooth(y_syn)

%%R -i x_syn -i y_syn -i x_plot -o y_plot_R -i df
cat(sprintf("R Timing (n=%d):\n", length(x_syn)))

# Fit once to inspect knots
fit_r_syn <- smooth.spline(x_syn, y_syn, df=df, all.knots=TRUE)
n_knots_r <- fit_r_syn$fit$nk
y_plot_R = predict(fit_r_syn, x_plot)$y
cat(sprintf("Number of knots used by R: %d\n", n_knots_r))

summary(microbenchmark(
  smooth.spline(x_syn, y_syn, df=df, all.knots=TRUE),
  times=100), unit='milliseconds')[,-1]

fig, ax = plt.subplots(figsize=(10, 6))
ax.scatter(x_syn, y_syn, 
            s=15, c='lightgray', alpha=0.7, label='Data')
ax.plot(x_plot, y_plot_py, 'b-', lw=3, label='Python (scatter_smooth)', alpha=0.8)
ax.plot(x_plot, y_plot_R, 'r--', lw=3, label='R (smooth.spline)', alpha=0.8)
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_title(f"Comparison of Smoothing Splines (df={df}) - Reinsch Engine")
ax.legend()
plt.show()

# Python Timing with reduced knots
df = 20
n_syn = 10000
x_syn = np.sort(rng.uniform(0, 10, n_syn))
y_syn = np.sin(np.sqrt(x_syn)*2) + rng.normal(0, 1, n_syn)
n_knots_reduced = 500 # Similar to R's behavior

# Fit for plotting
spline = SplineSmoother(x=x_syn, df=df, n_knots=n_knots_reduced)
spline.smooth(y_syn)
y_plot_py = spline.predict(x_plot)

%%timeit 
SplineSmoother(x=x_syn, df=df, n_knots=n_knots_reduced).smooth(y_syn)

%%R -i x_syn -i y_syn -i n_knots_reduced -o y_plot_R -i df
cat(sprintf("R Timing (n=%d):\n", length(x_syn)))

# Fit once to inspect knots
fit_r_syn <- smooth.spline(x_syn, y_syn, df=df, nknots=n_knots_reduced)
y_plot_R = predict(fit_r_syn, x_plot)$y
n_knots_r <- fit_r_syn$fit$nk
cat(sprintf("Number of knots used by R: %d\n", n_knots_r))

summary(microbenchmark(
  smooth.spline(x_syn, y_syn, df=df, nknots=n_knots_reduced),
  times=100), unit='milliseconds')[,-1]

fig, ax = plt.subplots(figsize=(10, 6))
# Subsample data for plotting if N is large
plot_idx = np.random.choice(len(x_syn), size=min(len(x_syn), 1000), replace=False)
ax.scatter(x_syn[plot_idx], y_syn[plot_idx], 
            s=5, c='lightgray', alpha=0.5, label='Data (Subsampled)')
ax.plot(x_plot, y_plot_py, 'b-', lw=3, label='Python (scatter_smooth)', alpha=0.8)
ax.plot(x_plot, y_plot_R, 'r--', lw=3, label='R (smooth.spline)', alpha=0.8)
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_title(f"Comparison of Smoothing Splines (df={df}) - BSpline Engine")
ax.legend()
plt.show()