Python中基于splrep和splev的三次样条插值误差求解问题

我正在对包含x、y数据及y方向误差的数据集执行加权三次样条拟合，之后希望对曲线上其他点进行插值并获取对应误差。但在实现误差传播计算时遇到困难，当前代码生成的插值误差存在零点，结果不符合预期，请问我哪里出错了？

当前实现代码：

import sys
import numpy as np
from PySide6.QtWidgets import QApplication, QMainWindow, QVBoxLayout, QWidget
from matplotlib.backends.backend_qt5agg import FigureCanvasQTAgg as FigureCanvas
from matplotlib.figure import Figure
from scipy.interpolate import splrep, splev

class MplCanvas(FigureCanvas):
    def __init__(self, parent=None, width=5, height=4, dpi=100):
        fig = Figure(figsize=(width, height), dpi=dpi)
        self.axes = fig.add_subplot(111)
        super().__init__(fig)

class MainWindow(QMainWindow):
    def __init__(self):
        super().__init__()

        self.canvas = MplCanvas(self, width=5, height=4, dpi=100)
        layout = QVBoxLayout()
        layout.addWidget(self.canvas)
        container = QWidget()
        container.setLayout(layout)
        self.setCentralWidget(container)

        self.plot()

    def plot(self):
        # Generate random data with errors
        np.random.seed(0)
        x = np.linspace(0, 10, 10)
        y = np.sin(x) + np.random.normal(0, 0.1, len(x))
        y_err = np.random.normal(0.1, 0.02, len(x)) + 0.5

        # Fit a weighted cubic spline using splrep
        tck = splrep(x, y, w=1/y_err, k=3)

        # Interpolate values using splev
        x_interp = np.linspace(0, 10, 100)
        y_interp = splev(x_interp, tck)

        # Calculate the fitted values at the data points (for residuals)
        y_fit = splev(x, tck)

        # Calculate residuals
        residuals = y - y_fit

        # Calculate chi-squared and degrees of freedom
        chi2 = np.sum((residuals / y_err) ** 2)
        dof = len(x) - len(tck[1]) // 3 - 1

        # Covariance matrix for weighted least squares
        cov_matrix = np.diag(y_err ** 2) * chi2 / dof

        # Compute Jacobians at interpolated points
        jacobian = np.array([splev(xi, tck, der=1) for xi in x_interp])

        # Extract the diagonal elements of the covariance matrix (variances)
        variances = np.diag(cov_matrix)

        # Convert variances to column vector for broadcasting
        variances_column = variances[:, np.newaxis]

        # Propagate errors through the squared Jacobian and sum up
        jacobian_squared = jacobian ** 2
        weighted_jacobian = jacobian_squared * variances_column
        sum_weighted_jacobian = np.sum(weighted_jacobian, axis=0)

        # Calculate the final interpolated errors
        y_interp_err = np.sqrt(sum_weighted_jacobian)

        # Plot the results
        self.canvas.axes.errorbar(x, y, yerr=y_err, fmt='o', label='Data')
        self.canvas.axes.plot(x_interp, y_interp, label='Cubic Spline Fit')
        self.canvas.axes.fill_between(x_interp, y_interp - y_interp_err, y_interp + y_interp_err, alpha=0.2, label='Error')
        self.canvas.axes.legend()
        self.canvas.draw()

app = QApplication(sys.argv)
window = MainWindow()
window.show()
sys.exit(app.exec())

核心错误点 #

1. 雅可比矩阵计算完全错误
   你当前计算的是样条曲线的一阶导数（der=1），但误差传播需要的是插值点拟合值对原始数据点y值的偏导数——即每个插值点的拟合结果是原始y数据的线性组合，需要计算这个组合的权重，而非曲线的斜率。

2. 协方差矩阵构建逻辑错误
   加权最小二乘的协方差矩阵不是直接用np.diag(y_err ** 2) * chi2 / dof。正确的计算应该基于权重矩阵的逆，而非原始误差的对角矩阵缩放，这完全违背了加权最小二乘的协方差推导逻辑。

修正后的实现思路 #

1. 获取样条基函数矩阵
   三次样条的拟合值是基函数与系数的线性组合，利用BSpline.design_matrix可以直接得到原始数据点和插值点对应的基函数矩阵，进而推导插值结果对原始y数据的偏导数（雅可比矩阵）。

2. 正确计算加权协方差
   先构建权重矩阵$W = \text{diag}(1/y_{err})$，计算$(B^T W B)^{{-1}$（$B$为原始数据的基函数矩阵），再结合$\chi}2/dof$得到系数的协方差，最后通过雅可比矩阵完成误差传播。

修正后的代码 #

import sys
import numpy as np
from PySide6.QtWidgets import QApplication, QMainWindow, QVBoxLayout, QWidget
from matplotlib.backends.backend_qt5agg import FigureCanvasQTAgg as FigureCanvas
from matplotlib.figure import Figure
from scipy.interpolate import splrep, BSpline

class MplCanvas(FigureCanvas):
    def __init__(self, parent=None, width=5, height=4, dpi=100):
        fig = Figure(figsize=(width, height), dpi=dpi)
        self.axes = fig.add_subplot(111)
        super().__init__(fig)

class MainWindow(QMainWindow):
    def __init__(self):
        super().__init__()

        self.canvas = MplCanvas(self, width=5, height=4, dpi=100)
        layout = QVBoxLayout()
        layout.addWidget(self.canvas)
        container = QWidget()
        container.setLayout(layout)
        self.setCentralWidget(container)

        self.plot()

    def plot(self):
        # Generate random data with errors
        np.random.seed(0)
        x = np.linspace(0, 10, 10)
        y = np.sin(x) + np.random.normal(0, 0.1, len(x))
        y_err = np.random.normal(0.1, 0.02, len(x)) + 0.5

        # Fit weighted cubic spline
        w = 1 / y_err
        tck = splrep(x, y, w=w, k=3)
        t, c, k = tck
        spline = BSpline(t, c, k)

        # Interpolate values
        x_interp = np.linspace(0, 10, 100)
        y_interp = spline(x_interp)
        y_fit = spline(x)

        # Calculate residuals and reduced chi-squared
        residuals = y - y_fit
        chi2 = np.sum((residuals * w) ** 2)
        dof = len(x) - len(c)  # 自由度：数据点数量 - 样条系数数量
        chi2_red = chi2 / dof

        # 构建原始数据点和插值点的基函数矩阵
        B = spline.design_matrix(x, t, k)
        B_interp = spline.design_matrix(x_interp, t, k)

        # 权重矩阵与协方差计算
        W = np.diag(w)
        BTWB = B.T @ W @ B
        BTWB_inv = np.linalg.inv(BTWB)

        # 雅可比矩阵：插值结果对原始y数据的偏导数
        J = B_interp @ BTWB_inv @ B.T @ W

        # 误差传播：计算插值结果的协方差矩阵
        cov_y = np.diag(y_err ** 2)
        cov_interp = J @ cov_y @ J.T
        y_interp_err = np.sqrt(np.diag(cov_interp))

        # 绘图
        self.canvas.axes.errorbar(x, y, yerr=y_err, fmt='o', label='Data')
        self.canvas.axes.plot(x_interp, y_interp, label='Cubic Spline Fit')
        self.canvas.axes.fill_between(x_interp, y_interp - y_interp_err, y_interp + y_interp_err, alpha=0.2, label='Error')
        self.canvas.axes.legend()
        self.canvas.draw()

app = QApplication(sys.argv)
window = MainWindow()
window.show()
sys.exit(app.exec())

额外说明 #

修正后的误差计算不再出现零点：原始数据点附近的插值误差会与原始数据误差接近，远离数据点时误差会合理增大，符合误差传播的物理意义。

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