如何计算三次样条插值后数据集各新元素的标准差？

问题背景 #

我已通过以下代码对数据集完成三次样条插值：

import numpy as np
from scipy import interpolate

y = [-2.35771490e-04, -2.04672372e-04,  5.40550849e-06, -3.26933980e-04,
       -4.31306957e-04, -4.06928662e-04, -7.50889655e-04, -1.03669983e-03,
       -5.28464040e-04, -5.01935670e-04, -8.09657243e-04, -4.21891416e-04,
       -1.35441615e-04, -3.46855608e-04, -2.90791620e-04, -8.19357356e-05,
       -2.58070932e-04, -1.58262175e-04,  1.32105126e-04,  2.42572582e-05,
        9.88706879e-05]

x = [175., 165., 155., 145., 135., 125., 115., 105.,  95.,  85.,  75.,
    65.,  55.,  45.,  35.,  25.,  15.,   5.,  -5., -15., -25.]

f = interpolate.interp1d(x, y, kind='cubic')
x_new = np.arange(-25,175, 5)
y_interp = f(x_new)

现希望计算y_interp中每个新元素的标准差，原始y元素的误差如下：

err = [0.00012969551729667857, 0.00014864077922340332, 0.00010732361332456688, 8.95365810198269e-05, 8.559972679579093e-05, 0.00010669277818690999, 6.880710200343582e-05, 9.249378243164652e-05, 9.885179214527947e-05, 8.352246348207366e-05, 6.586538652949116e-05, 7.428127049298688e-05, 6.506191339534108e-05, 7.166116284348538e-05, 6.379032268503384e-05, 7.101008284147194e-05, 6.0593495968463165e-05, 4.699920194672463e-05, 5.854938570487613e-05, 7.140629409894438e-05, 9.66241238445066e-05]

解决方案 #

可以计算插值结果的标准差，核心逻辑是利用三次样条的线性性质：插值结果是原始样本点的线性组合，通过计算每个插值点对原始样本的权重，结合原始样本的误差推导插值结果的标准差。

步骤1：改用CubicSpline获取插值权重 #

interp1d没有直接返回权重的接口，改用scipy.interpolate.CubicSpline可以更方便地计算插值权重：

from scipy.interpolate import CubicSpline
import numpy as np

# 原始数据
y = [-2.35771490e-04, -2.04672372e-04,  5.40550849e-06, -3.26933980e-04,
        -4.31306957e-04, -4.06928662e-04, -7.50889655e-04, -1.03669983e-03,
        -5.28464040e-04, -5.01935670e-04, -8.09657243e-04, -4.21891416e-04,
        -1.35441615e-04, -3.46855608e-04, -2.90791620e-04, -8.19357356e-05,
        -2.58070932e-04, -1.58262175e-04,  1.32105126e-04,  2.42572582e-05,
         9.88706879e-05]
x = [175., 165., 155., 145., 135., 125., 115., 105.,  95.,  85.,  75.,
        65.,  55.,  45.,  35.,  25.,  15.,   5.,  -5., -15., -25.]
err = [0.00012969551729667857, 0.00014864077922340332, 0.00010732361332456688, 8.95365810198269e-05, 8.559972679579093e-05, 0.00010669277818690999, 6.880710200343582e-05, 9.249378243164652e-05, 9.885179214527947e-05, 8.352246348207366e-05, 6.586538652949116e-05, 7.428127049298688e-05, 6.506191339534108e-05, 7.166116284348538e-05, 6.379032268503384e-05, 7.101008284147194e-05, 6.0593495968463165e-05, 4.699920194672463e-05, 5.854938570487613e-05, 7.140629409894438e-05, 9.66241238445066e-05]

# 构造三次样条（默认自然边界，和interp1d的cubic一致）
cs = CubicSpline(x, y)
x_new = np.arange(-25, 175, 5)
y_interp = cs(x_new)

步骤2：计算插值权重矩阵 #

通过对单位向量做插值，得到每个插值点对原始样本的权重：

def get_interpolation_weights(x, x_new):
    n = len(x)
    # 对每个单位向量执行插值，得到权重矩阵
    cs = CubicSpline(x, np.eye(n))
    return cs(x_new)

# 获取权重矩阵：每行对应一个x_new点，每列对应一个原始x点的权重
weights = get_interpolation_weights(x, x_new)

步骤3：计算插值结果的标准差 #

假设原始err是每个y元素的标准差，且样本误差独立，那么插值结果的标准差计算公式为：
$\text{std}(y_{\text{interp}}) = \sqrt{\sum_{i=1}^n w_i^2 \times \text{err}_i^2}$

代码实现：

# 计算原始y的方差（标准差的平方）
var_y = np.array(err) ** 2
# 计算每个插值点的方差
var_interp = np.sum(weights ** 2 * var_y, axis=1)
# 得到插值结果的标准差
std_interp = np.sqrt(var_interp)

注意事项 #

• 上述方法基于原始样本误差独立的假设，如果误差之间存在相关性，需要用协方差矩阵$\Sigma$计算，公式为$\text{Var}(y_{\text{interp}}) = A \Sigma A^T$，其中$A$是权重矩阵。
• CubicSpline默认使用自然边界（端点二阶导数为0），和interp1d(kind='cubic')的边界条件一致，确保插值结果匹配。

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