如何用Scipy拟合单调递增的光滑样条函数

需要为数据集拟合单调递增的光滑样条函数，期望达到类似Matlab中单调递增样条的拟合效果。

已尝试的方法及问题 #

1. scipy.interpolate.interp1d三次样条拟合 #

代码实现：

from scipy.interpolate import interp1d
import matplotlib.pyplot as plt

x = [0., 0.75, 1.8, 2.25, 3.75, 4.5, 6.45, 6.75, 7.5, 8.325, 10.875, 11.25, 12.525, 12.75, 15., 20.85, 21.]
y = [2.83811035, 2.81541896, 3.14311655, 3.22373554, 3.43033456, 3.50433385, 3.66794514, 3.462296, 3.59480959,
     3.56250726, 3.6209845,  3.63034523, 3.68238915, 3.69096892, 3.75560395, 3.83545191, 3.90419498]

plt.plot(x, y, '*')
plt.show()
f = interp1d(x, y, kind='cubic')
yinp = f(x)
plt.plot(x, yinp)
plt.show()

问题：拟合结果不满足单调递增要求。

2. scipy.interpolate.UnivariateSpline光滑样条拟合 #

代码实现：

import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import UnivariateSpline

x = np.array([0., 0.75, 1.8, 2.25, 3.75, 4.5, 6.45, 6.75, 
              7.5, 8.325, 10.875, 11.25, 12.525, 12.75, 15., 
              20.85, 21.])
y = np.array([2.83811035, 2.81541896, 3.14311655, 
              3.22373554, 3.43033456, 3.50433385, 
              3.66794514, 3.462296,   3.59480959, 
              3.56250726, 3.6209845,  3.63034523, 
              3.68238915, 3.69096892, 3.75560395, 
              3.83545191, 3.90419498])
              

spl = UnivariateSpline(x, y, s=0.05)
xs  = np.linspace(x.min(), x.max(), 100)
plt.plot(x, y, 'ro', ms=5)
plt.plot(xs, spl(xs), 'cyan', lw=1)
plt.show()

通过以下代码获取样条系数与节点：

print(spl.get_coeffs())
print(spl.get_knots())

k = 3
tck = splrep(x, y, k=k, s=0.09)
knots = tck[0]
coeffs = tck[1]
print('knot points=', knots)
print('coefficients=', coeffs)

问题：无法利用系数手动生成分段多项式函数，例如针对小数据集：

x = [0., 0.75, 1.8, 2.25]
y = [2.83811035, 2.81541896, 3.14311655, 3.22373554]

需要输出对应的分段多项式以理解样条结构。

解决方案 #

一、Scipy中实现单调递增光滑样条 #

1. 使用PchipInterpolator（推荐）

PCHIP（分段三次Hermite插值多项式）是保形插值方法，默认保证单调性，同时保持一阶连续可导的光滑性，效果接近Matlab的单调样条。

代码示例：

import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import PchipInterpolator

x = np.array([0., 0.75, 1.8, 2.25, 3.75, 4.5, 6.45, 6.75, 
              7.5, 8.325, 10.875, 11.25, 12.525, 12.75, 15., 
              20.85, 21.])
y = np.array([2.83811035, 2.81541896, 3.14311655, 
              3.22373554, 3.43033456, 3.50433385, 
              3.66794514, 3.462296,   3.59480959, 
              3.56250726, 3.6209845,  3.63034523, 
              3.68238915, 3.69096892, 3.75560395, 
              3.83545191, 3.90419498])

# 初始化PCHIP插值器，自动保证单调性
pchip = PchipInterpolator(x, y)
xs = np.linspace(x.min(), x.max(), 200)
plt.plot(x, y, 'ro', ms=5)
plt.plot(xs, pchip(xs), 'b-', lw=2)
plt.title('Monotonic PCHIP Spline')
plt.show()

2. 带约束的自定义光滑样条

若需要更灵活的光滑度控制，可通过scipy.optimize实现：以最小化样条曲率（保证光滑）为目标，同时约束每个分段的导数非负。这种方式实现复杂，适合对光滑度有特殊要求的场景。

二、手动生成分段多项式函数 #

1. 基于UnivariateSpline的分段多项式推导

三次样条在每个节点区间[t_i, t_{i+1}]的形式为：
$$S_i(x) = a_i + b_i(x-t_i) + c_i(x-t_i)^2 + d_i(x-t_i)^3$$
可通过区间内点反推系数：

import numpy as np
from scipy.interpolate import splrep, splev

# 小数据集示例
x = np.array([0., 0.75, 1.8, 2.25])
y = np.array([2.83811035, 2.81541896, 3.14311655, 3.22373554])

# 获取样条的tck元组（k=3为三次样条，s=0表示插值过所有点）
tck = splrep(x, y, k=3, s=0)
knots, coeffs, degree = tck

# 提取唯一节点区间
unique_knots = np.unique(knots)
segments = list(zip(unique_knots[:-1], unique_knots[1:]))

# 计算每个分段的多项式系数
for i, (t_start, t_end) in enumerate(segments):
    # 在区间内取4个点构建方程组
    x_segment = np.linspace(t_start, t_end, 4)
    y_segment = splev(x_segment, tck)
    dx = x_segment - t_start
    # 构建线性方程组求解a,b,c,d
    A = np.vstack([np.ones(4), dx, dx**2, dx**3]).T
    a, b, c, d = np.linalg.solve(A, y_segment)
    
    print(f"分段 {i+1}: [{t_start:.2f}, {t_end:.2f}]")
    print(f"多项式: S(x) = {a:.6f} + {b:.6f}(x-{t_start:.2f}) + {c:.6f}(x-{t_start:.2f})² + {d:.6f}(x-{t_start:.2f})³\n")

2. 基于PchipInterpolator的分段多项式推导

PCHIP本身基于Hermite插值，可直接通过节点导数构建分段多项式：

from scipy.interpolate import PchipInterpolator

# 小数据集示例
x = np.array([0., 0.75, 1.8, 2.25])
y = np.array([2.83811035, 2.81541896, 3.14311655, 3.22373554])

pchip = PchipInterpolator(x, y)
# 获取每个节点的导数
slopes = pchip.derivative()(x)

for i in range(len(x)-1):
    x0, x1 = x[i], x[i+1]
    y0, y1 = y[i], y[i+1]
    m0, m1 = slopes[i], slopes[i+1]
    
    dx = x1 - x0
    dy = y1 - y0
    # 计算三次Hermite多项式系数
    c0 = y0
    c1 = m0
    c2 = (3*dy/dx - 2*m0 - m1)/dx
    c3 = (m0 + m1 - 2*dy/dx)/(dx**2)
    
    print(f"分段 {i+1}: [{x0:.2f}, {x1:.2f}]")
    print(f"多项式: S(x) = {c0:.6f} + {c1:.6f}(x-{x0:.2f}) + {c2:.6f}(x-{x0:.2f})² + {c3:.6f}(x-{x0:.2f})³\n")

内容的提问来源于stack exchange，提问作者Natasha