Как получить ошибки при интерполяции кубическим сплайном (Python; splrep, splev) ⇐ Python

[Anonymous]

Я выполняю взвешенный кубический сплайн, аппроксимирующий некоторые данные 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())

При этом создаются точки, в которых ошибки равны нулю, например:
[img]https://i. sstatic.net/cWGWTVCg.png[/img]

что выглядит неправильно. Что я делаю не так?

Подробнее здесь: https://stackoverflow.com/questions/790 ... lrep-splev