Python 中的加权最小二乘法

Question

我正在尝试对我的数据进行加权最小二乘拟合。当 x=0 和 x>0.2 时的值应给予较大的权重。 0时的值

我的数据是：

filtered_x = array([0.        , 0.03807991, 0.05077321, 0.06346652, 0.07615982,
       0.08885313, 0.10154643, 0.11423973, 0.12693304, 0.13962634,
       0.15231964, 0.16501295, 0.17770625, 0.19039955, 0.20309286,
       0.21578616, 0.22847947, 0.24117277, 0.25386607, 0.26655938,
       0.27925268, 0.29194598, 0.30463929, 0.31733259, 0.33002589,
       0.3427192 ])

filtered_y = array([0.        , 1.53989397, 2.04460628, 4.18043213, 2.97621482,
       2.82642339, 2.98335023, 2.98964836, 2.12218901, 1.42801972,
       1.25930683, 0.71644077, 0.48220866, 0.21165985, 0.24756609,
       0.21123179, 0.57344999, 0.49362762, 0.20282767, 0.50321186,
       0.50347165, 0.74259408, 0.48493783, 0.81785588, 0.54543666,
       0.53218838])

我已经有以下代码：

# fitted function
def fpow(x, a, b, c):
    return a*x**b+c

# custom weight function
def custom_weights(x):
    weights = np.ones_like(x)  # Default weight for all points
    weights[x == 0] = 0.1     # Low Error for x = 0
    weights[x > 0.2] = 0.1   # Low Error for x > 0.15
    weights[(0 < x) & (x <= 0.2)] = 0.5 # Medium error for 0 < x <= 0.2
    return weights
# Weight 0 and higher wavenumbers heavier

# initial guess
pars0 = (0.4, 0.4, 1)

# perform fitting with custom weights
popt, pcov = curve_fit(fpow, filtered_x, filtered_y, absolute_sigma=True, p0=pars0, sigma=custom_weights(filtered_x), maxfev =5000)

# parameters
a_opt = popt[0]
b_opt = popt[1]
c_opt = popt[2]

# plot data
plt.errorbar(filtered_x, filtered_y, yerr=0, fmt =".", color='black', label ='data', zorder=1, markersize=5)

# creating x interval to include in y fit
x_interval = np.linspace(0, max(filtered_x), 1000)
y_fit = fpow( x_interval , *popt)
plt.plot( x_interval, y_fit, color = "red", label = "Weighted fit", zorder=2 , linewidth=3)

plt.grid(True)
plt.ylabel("U [m/s]")
plt.xlabel("Wavenumber [rad/m]")
plt.title("LS Weighted Fit of Current")
plt.legend()
plt.show()

我得到以下内容

它还给了我一个运行时错误

<ipython-input-747-a3ac7b8135e0>:2: RuntimeWarning: divide by zero encountered in power
  return a*x**b+c

您可以看到出了问题，对于较大的 x 值拟合效果不好。我尝试删除数据中的 (0,0) 点，但是，我得到了更糟糕的结果，如下所示：

Answer 1

这里有一个关于这种

RuntimeError

的问题，答案建议对数据使用非常小的值而不是零。

filtered_x = array([1e-8, 0.03807991, 0.05077321, 0.06346652, 0.07615982,
                    0.08885313, 0.10154643, 0.11423973, 0.12693304, 0.13962634,
                    0.15231964, 0.16501295, 0.17770625, 0.19039955, 0.20309286,
                    0.21578616, 0.22847947, 0.24117277, 0.25386607, 0.26655938,
                    0.27925268, 0.29194598, 0.30463929, 0.31733259, 0.33002589,
                    0.3427192])

filtered_y = array([1e-8, 1.53989397, 2.04460628, 4.18043213, 2.97621482,
                    2.82642339, 2.98335023, 2.98964836, 2.12218901, 1.42801972,
                    1.25930683, 0.71644077, 0.48220866, 0.21165985, 0.24756609,
                    0.21123179, 0.57344999, 0.49362762, 0.20282767, 0.50321186,
                    0.50347165, 0.74259408, 0.48493783, 0.81785588, 0.54543666,
                    0.53218838])

注意第一个位置的

1e-8

。这样就消除了警告。 它还会在一定程度上改善贴合度。

但请查看另一个问题，因为它建议了更多方法来改善贴合度。

编辑：

我忘了提及，我更改了您的自定义权重函数，因此它使用

np.isclose

。检查浮点数是否相等很容易出错，而且通常不是您想要的。

np.isclose

检查数字是否处于要比较的值的给定 epsilon 环境中。 stdlib 中有一个类似的函数

math.isclose

。

# custom weight function
def custom_weights(x):
    weights = np.ones_like(x)  # Default weight for all points
    weights[np.isclose(x, 0, atol=1e-7)] = 0.1  # Low Error for x = 0
    weights[x > 0.2] = 0.1  # Low Error for x > 0.15
    weights[(0 < x) & (x <= 0.2)] = 0.5  # Medium error for 0 < x <= 0.2
    return weights