Python中非均匀数据卷积拟合的插值权重优化问题

问题背景 #

我需要对非均匀分布的时间分辨吸收数据（数据点呈指数分布）进行曲线拟合，但拟合过程涉及卷积。我的理解是：卷积要求数组具有均匀间距，因此必须先对数据进行插值——这个理解是否正确？

插值后出现的问题是：原数据点稀少的区域会生成大量插值点，导致这些区域在拟合中的权重被过度放大。

目前的思路是用Python的curve_fit设置权重来抵消失衡，但不确定具体实现方式，想知道这类问题有没有标准解决方案。

补充说明与尝试 #

由于完整代码过大，以下用线性函数示例演示（实际问题需卷积，因此必须插值）：

• 示例中x值前半部分间距大，后半部分间距小
• 分别做了三种拟合：无插值拟合、插值后拟合、插值加权重拟合

权重设置思路 #

对数据点密度过高的区域降低权重：若原区域仅1个数据点，插值后生成10个点，则设置sigma = 插值步长/原区域步长，期望还原无插值时的拟合结果，但实际未实现。

拟合结果 #

• 仅插值的拟合结果误差过小（因新增大量邻近点，符合预期）
• 加权拟合结果虽不总是接近无插值拟合，但误差处于同一数量级

示例代码 #

import numpy as np
import matplotlib.pyplot as plt

def interpolate_data(x_data, y_data, stepsize=0.5):
    """
    Interpolates given data and returns interpolated values at x-steps of 0.5.

    Parameters:
    - x_data: Array of x values (sorted in ascending order).
    - y_data: Array of corresponding y values.

    Returns:
    - Interpolated x values at steps of stepsize.
    - Interpolated y values at corresponding x values.
    """
    # Define the target x values with steps of stepsize
    x_interp = np.arange(np.min(x_data), np.max(x_data) + stepsize, stepsize)

    # Interpolate y values corresponding to the x_interp values
    y_interp = np.interp(x_interp, x_data, y_data)

    return x_interp, y_interp

def linear(x,a,b):
    return a*x+b

################################## Generate data ###############################################
noise_amplitude = 10
# Parameters for the linear function
a = 3.5
b = 2.0
# Define step sizes
x_step1 = 1
x_step2 = 0.1
# Length of each half
half_length = 10
# First half with step size x_step1
x_first_half = np.arange(0, half_length, x_step1)
# Second half with step size x_step2
x_second_half = np.arange(half_length, 2*half_length, x_step2)
# Combine both halves
x = np.concatenate((x_first_half, x_second_half))
# print("Generated x array:")
# print(x)

# Generate y values without noise
y_true = linear(x, a, b)
# Add normally distributed noise to y values
noise = np.random.normal(0, noise_amplitude, size=x.size)
y_data = y_true + noise

interpolation_stepsize = x_step2
x_ip, y_ip =interpolate_data(x,y_data,stepsize=interpolation_stepsize)
########################################## Fit ############################################
%matplotlib inline
from scipy.optimize import curve_fit
###################################### No interpolation ###############################
popt, pcov = curve_fit(linear, x, y_data)
# Extract the fitted parameters and their errors (standard deviations)
a_fit, b_fit = popt
std_err_a, std_err_b = np.sqrt(np.diag(pcov))
# Print the fitted parameters and their errors 
print(f"Fitted parameters (No interpolation):")
print(f"  a = {a_fit:.3f} ± {std_err_a:.3f}")
print(f"  b = {b_fit:.3f} ± {std_err_b:.3f}")
###################################### With interpolation ###############################
# Perform curve fitting using curve_fit
popt_ip, pcov_ip = curve_fit(linear, x_ip, y_ip)
# Extract the fitted parameters and their errors (standard deviations)
a_fit, b_fit = popt_ip
std_err_a, std_err_b = np.sqrt(np.diag(pcov_ip))
# print the fitted parameters and their errors
print(f"Fitted parameters (With interpolation):")
print(f"  a = {a_fit:.3f} ± {std_err_a:.3f}")
print(f"  b = {b_fit:.3f} ± {std_err_b:.3f}")
###################################### With interpolation and weights ###############################
# Attempt 1: weigh down / weigh up by the factor interpolation_stepsize/x_region. So that regions that have to many points 
# due to interpolation have less weight and vice versa.
## Construct weight array
# First half with step size x_step1
w_first_half = np.ones(int(len(x_ip)/2))*interpolation_stepsize/x_step1
# Second half with step size x_step2
w_second_half =     np.ones(int(len(x_ip)/2))*interpolation_stepsize/x_step2
# Combine both halves
weigth_array = np.concatenate((w_first_half,w_second_half))
###############################################
# Perform curve fitting using curve_fit
popt_ip_weighted, pcov_ip_weighted = curve_fit(linear, x_ip, y_ip,sigma=1/weigth_array)

# Extract the fitted parameters and their errors (standard deviations)
a_fit, b_fit = popt_ip_weighted
std_err_a, std_err_b = np.sqrt(np.diag(pcov_ip_weighted))

# Print the fitted parameters and their errors
print(f"Fitted parameters (With interpolation and weights):")
print(f"  a = {a_fit:.3f} ± {std_err_a:.3f}")
print(f"  b = {b_fit:.3f} ± {std_err_b:.3f}")

################################ Plot the data and the fitted models ######################################


plt.scatter(x, y_data, label='Data with noise')
plt.scatter(x_ip, y_ip, marker="x", s=10,label='Interpolated Data with noise')

plt.plot(x, y_true, color='k', linestyle='--', label='True linear model')
plt.plot(x, linear(x, *popt), color='green', label='Fitted linear model')
plt.plot(x_ip, linear(x_ip, *popt_ip), color='red', label='Fitted linear model of interpolated data')
plt.plot(x_ip, linear(x_ip, *popt_ip_weighted), color='blue', label='Fitted linear model of weighted interpolated data')

plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.title('Curve fitting with linear model')
plt.grid(True)
plt.show()

示例输出 #

无插值拟合参数：
  a = 3.164 ± 0.232
  b = 6.627 ± 3.382
插值后拟合参数：
  a = 2.873 ± 0.136
  b = 9.113 ± 1.566
插值加权重拟合参数：
  a = 3.279 ± 0.222
  b = 5.053 ± 3.364

拟合结果图 #

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