使用scipy.interpolate.splrep拟合线性样条出现大量NaN系数求助

我有一组数据，手动可轻松拟合带单个内部节点的线性样条，但使用scipy.interpolate.splrep拟合时，无论设置何种平滑参数都无法正常工作，系数中出现大量NaN。编写迭代尝试不同平滑参数的代码直到抛出RuntimeWarning，结果仍存在大量NaN；为初始点设置高权重也无效。手动指定节点（如t=[75.0]）时拟合可成功，但需要无需人工干预的可重复自动化方案，不限于scipy库。

相关代码：

import math
import itertools
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
from scipy.interpolate import splrep, splev


def is_sorted(a):
    return np.all(a[:-1] <= a[1:])


def fit_spline(x, y):
    assert(is_sorted(x))
    
    best_n_nans = 2*len(x)
    for log_s in itertools.count(0.01, 0.0001):
        try:
            spl = splrep(x, y, s=math.exp(log_s), k=1)
        except RuntimeWarning:
            break
        
        n_nans = sum(math.isnan(t) for t in spl[1])
        if n_nans < best_n_nans:
            best_fit = spl
            best_n_nans = n_nans
            print(f"Found fit with {n_nans} nans")
            
        if len(spl[0]) <= 5:
            break
    return best_fit

matplotlib.use("Qt5Agg")
plot = plt.figure()
x = np.array([  0. ,   0. ,   0. ,   0. ,   0. ,   0.5,   1.2,   2.3,   2.6,
         3. ,   3.1,   3.1,   3.1,   5. ,   5.1,   5.2,   5.6,   6.3,
         6.5,   6.9,   9.6,   9.8,  10. ,  10.3,  10.7,  11.3,  11.3,
        12.4,  12.4,  12.4,  12.5,  12.5,  13.1,  14.2,  15.4,  15.7,
        16.1,  17.7,  18.5,  18.6,  20.5,  21.3,  22.7,  23. ,  24.3,
        26.3,  26.3,  26.6,  27. ,  27.7,  28.7,  31. ,  31. ,  31.1,
        31.1,  31.6,  32.7,  34.4,  34.5,  35.7,  35.7,  38.7,  39.8,
        43.9,  43.9,  46.9,  51. ,  52.6,  52.6,  53.3,  53.7,  54.5,
        57.7,  61.2,  62.5,  65. ,  66. ,  68.1,  70.8,  74.3,  74.6,
        76.3,  83.1,  85.7,  87.5,  90.8,  91. ,  92.2,  94.3,  94.6,
       101.1, 104.8, 109.6, 113.3, 114.8, 119. , 120.5, 121.1, 126.1,
       133.2, 136.2, 137.1, 138.5, 141.4, 144.6, 145.1, 145.3, 150.5,
       160.7, 167.4, 168.8, 170.5, 170.5, 173.3, 177. , 181.3, 185.1,
       193. , 200.3, 201.8, 204.3, 208.4, 209.9, 223.7, 230.9, 240.3,
       243.9, 250.6, 263.1])
y = np.array([ 0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,
        0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,
        0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,
        0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  7,  6,
        0,  0,  0,  9,  0,  0,  0,  0,  0,  0,  9,  1,  0,  3,  0,  9,  2,
       10,  0, 14, 19,  7, 10, 18, 22, 31, 12, 29, 16, 20,  2, 38, 26, 28,
       26, 25, 17, 35, 20, 29, 39, 51, 29, 27, 43, 38, 40, 47, 41, 39, 54,
       44, 53, 65, 50, 59, 51, 62, 69, 74, 66])

spl = fit_spline(x, y)
y2 = splev(x, spl)
plt.plot(x, y, 'o', x, y2)
plot.show()

拟合得到的系数存在大量NaN：

array([nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan,
       nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan,
       nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan,
       nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, 35., 38.,
       40., 47., 41., 39., 54., 44., 53., 65., 50., 59., 51., 62., 69.,
       74., 66.,  0.,  0.])

方案1：基于数据特征自动找拐点，手动构造线性样条 #

你的数据前半段y全为0，后半段y开始上升，拐点非常明显。可以直接通过检测第一个非零y值对应的x位置作为内部节点，然后分两段做线性拟合：

import numpy as np
import matplotlib.pyplot as plt

# 原始数据
x = np.array([  0. ,   0. ,   0. ,   0. ,   0. ,   0.5,   1.2,   2.3,   2.6,
         3. ,   3.1,   3.1,   3.1,   5. ,   5.1,   5.2,   5.6,   6.3,
         6.5,   6.9,   9.6,   9.8,  10. ,  10.3,  10.7,  11.3,  11.3,
        12.4,  12.4,  12.4,  12.5,  12.5,  13.1,  14.2,  15.4,  15.7,
        16.1,  17.7,  18.5,  18.6,  20.5,  21.3,  22.7,  23. ,  24.3,
        26.3,  26.3,  26.6,  27. ,  27.7,  28.7,  31. ,  31. ,  31.1,
        31.1,  31.6,  32.7,  34.4,  34.5,  35.7,  35.7,  38.7,  39.8,
        43.9,  43.9,  46.9,  51. ,  52.6,  52.6,  53.3,  53.7,  54.5,
        57.7,  61.2,  62.5,  65. ,  66. ,  68.1,  70.8,  74.3,  74.6,
        76.3,  83.1,  85.7,  87.5,  90.8,  91. ,  92.2,  94.3,  94.6,
       101.1, 104.8, 109.6, 113.3, 114.8, 119. , 120.5, 121.1, 126.1,
       133.2, 136.2, 137.1, 138.5, 141.4, 144.6, 145.1, 145.3, 150.5,
       160.7, 167.4, 168.8, 170.5, 170.5, 173.3, 177. , 181.3, 185.1,
       193. , 200.3, 201.8, 204.3, 208.4, 209.9, 223.7, 230.9, 240.3,
       243.9, 250.6, 263.1])
y = np.array([ 0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,
        0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,
        0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,
        0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  7,  6,
        0,  0,  0,  9,  0,  0,  0,  0,  0,  0,  9,  1,  0,  3,  0,  9,  2,
       10,  0, 14, 19,  7, 10, 18, 22, 31, 12, 29, 16, 20,  2, 38, 26, 28,
       26, 25, 17, 35, 20, 29, 39, 51, 29, 27, 43, 38, 40, 47, 41, 39, 54,
       44, 53, 65, 50, 59, 51, 62, 69, 74, 66])

# 自动定位拐点：第一个y>0的位置对应的x值
non_zero_idx = np.argmax(y > 0)
knot_x = x[non_zero_idx]

# 拆分数据为两段
x_segment1 = x[x <= knot_x]
y_segment1 = y[x <= knot_x]
x_segment2 = x[x >= knot_x]
y_segment2 = y[x >= knot_x]

# 分别拟合线性模型
coeffs1 = np.polyfit(x_segment1, y_segment1, 1)
coeffs2 = np.polyfit(x_segment2, y_segment2, 1)

# 生成拟合曲线
x_fit = np.linspace(x.min(), x.max(), 500)
y_fit = np.where(x_fit <= knot_x,
                 np.polyval(coeffs1, x_fit),
                 np.polyval(coeffs2, x_fit))

# 可视化
plt.figure()
plt.plot(x, y, 'o', label='原始数据')
plt.plot(x_fit, y_fit, '-', linewidth=2, label='自动拟合线性样条')
plt.legend()
plt.show()

这个方案完全自动化，利用数据的特征直接定位节点，拟合结果和手动操作一致，不会出现NaN问题。

方案2：使用scipy.interpolate.LSQUnivariateSpline约束节点数量 #

如果偏好使用scipy的样条工具，LSQUnivariateSpline可以指定固定数量的内部节点，避免splrep的自动节点选择导致的NaN问题。通过一阶差分自动定位初始节点，再让算法优化到最优位置：

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

# 原始数据
x = np.array([  0. ,   0. ,   0. ,   0. ,   0. ,   0.5,   1.2,   2.3,   2.6,
         3. ,   3.1,   3.1,   3.1,   5. ,   5.1,   5.2,   5.6,   6.3,
         6.5,   6.9,   9.6,   9.8,  10. ,  10.3,