Skip to content

Instantly share code, notes, and snippets.

@ahwillia
Last active December 14, 2023 10:21
Show Gist options
  • Star 2 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save ahwillia/511097a0968bf05a2579db0eab353393 to your computer and use it in GitHub Desktop.
Save ahwillia/511097a0968bf05a2579db0eab353393 to your computer and use it in GitHub Desktop.
Generate M-spline functions in Python
"""
Python code to generate M-splines.
References
----------
Ramsay, J. O. (1988). Monotone regression splines in action.
Statistical science, 3(4), 425-441.
"""
import numpy as np
import matplotlib.pyplot as plt
def mspline_grid(order, num_basis_funcs, nt):
"""
Generate a set of M-spline basis functions with evenly
spaced knots.
Parameters
----------
order : int
Order parameter of the splines.
num_basis_funcs : int
Number of desired basis functions. Note that we
require num_basis_funcs >= order.
nt : int
Number of points to evaluate the basis functions.
Returns
-------
spine_basis : array
Matrix with shape (num_basis_funcs, nt), holding the
desired spline basis functions.
"""
# Determine number of interior knots.
num_interior_knots = num_basis_funcs - order
if num_interior_knots < 0:
raise ValueError(
"Spline `order` parameter cannot be larger "
"than `num_basis_funcs` parameter."
)
# Fine grid of numerically evaluated points.
x = np.linspace(0, 1 - 1e-6, nt)
# Set of spline knots. We need to add extra knots to
# the end to handle boundary conditions for higher-order
# spline bases. See Ramsay (1988) cited above.
#
# Note - this is poorly explained on most corners of the
# internet that I've found.
knots = np.concatenate((
np.zeros(order - 1),
np.linspace(0, 1, num_interior_knots + 2),
np.ones(order - 1),
))
# Evaluate and stack each basis function.
return np.row_stack(
[mspline(x, order, i, knots) for i in range(num_basis_funcs)]
)
def mspline(x, k, i, T):
"""
Compute M-spline basis function `i` at points `x` for a spline
basis of order-`k` with knots `T`.
Parameters
----------
x : array
Vector holding points to evaluate the spline.
"""
# Boundary conditions.
if (T[i + k] - T[i]) < 1e-6:
return np.zeros_like(x)
# Special base case of first-order spline basis.
elif k == 1:
v = np.zeros_like(x)
v[(x >= T[i]) & (x < T[i + 1])] = 1 / (T[i + 1] - T[i])
return v
# General case, defined recursively
else:
return k * (
(x - T[i]) * mspline(x, k - 1, i, T)
+ (T[i + k] - x) * mspline(x, k - 1, i + 1, T)
) / ((k-1) * (T[i + k] - T[i]))
# Test code
if __name__ == "__main__":
# Plot basis funcs, varying the order parameter.
fig, axes = plt.subplots(1, 5)
for k, ax in enumerate(axes):
order = k + 1
ax.plot(mspline_grid(order, 7, 1000).T)
ax.set_yticks([])
ax.set_title(f"order-{order}")
fig.tight_layout()
plt.show()
@ahwillia
Copy link
Author

image

@MichaelTwiton
Copy link

This is great. I was also a bit confused about the boundary conditions in the paper.
I think you have an extra factor of k in line 91 (you already accounted for it in line 87).
Thanks!

@ahwillia
Copy link
Author

Thanks -- good catch!

@MichaelTwiton
Copy link

You don't happen to have a reference on implementing I-splines from the same paper in Python by any chance, do you?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment