如何使用 numy linalg lstsq 来拟合具有相同斜率但截距不同的两个数据集？

Question

我正在尝试进行加权最小二乘拟合，并遇到了numpy.linalg.lstsq。我需要拟合加权最小二乘。因此，以下工作：

# Generate some synthetic data from the model.
N = 50
x = np.sort(10 * np.random.rand(N))
yerr = 0.1 + 0.5 * np.random.rand(N)
y = 10.0 * x + 15
y += yerr * np.random.randn(N)
#do the fitting
err = 1/yerr**2
W = np.sqrt(np.diag(err))
x = x.flatten()
y = y.flatten()
A = np.vstack([x, np.ones(len(x))]).T
xw = np.dot(W,A)
yw = np.dot(W,y)
m, b = np.linalg.lstsq(xw, yw)[0]

这给了我最合适的斜率和截距。现在，假设我有两个具有相同斜率但截距不同的数据集？我将如何进行联合拟合以获得最佳拟合斜率加上两个截距。我仍然需要加权最小二乘版本。对于未加权的情况，我发现以下工作：

(m,b1,b2),_,_,_ = np.linalg.lstsq(np.stack([np.concatenate((x1,x2)),
                                        np.concatenate([np.ones(len(x1)),np.zeros(len(x2))]),
                                        np.concatenate([np.zeros(len(x1)),np.ones(len(x2))])]).T, 
                              np.concatenate((y1,y2)))

Answer 1

首先我重写你的第一种方法，因为我认为它可以写得更清楚

weights = 1 / yerr
m, b = np.linalg.lstsq(np.c_[weights * x, weights], weights * y, rcond=None)[0]

为了适应 2 个数据集，您可以堆叠 2 个数组，但将矩阵的一些元素设为 0。

np.random.seed(12)
N = 3
x = np.sort(10 * np.random.rand(N))
yerr = 0.1 + 0.5 * np.random.rand(N)
y = 10.0 * x + 15
y += yerr * np.random.randn(N)

M = 2
x1 = np.sort(10 * np.random.rand(M))
yerr1 = 0.1 * 0.5 * np.random.rand(M)
y1 = 10.0 * x1 + 25
y1 += yerr1 * np.random.randn(M)
#do the fitting
weights = 1 / yerr
weights1 = 1 / yerr1
first_column = np.r_[weights * x, weights1 * x1]
second_column = np.r_[weights, [0] * x1.size]
third_column = np.r_[[0] * x.size, weights1]
a = np.c_[first_column, second_column, third_column]
print(a)
# [[  4.20211437   2.72576342   0.        ]
#  [ 24.54293941   9.32075195   0.        ]
#  [ 13.22997409   1.78771428   0.        ]
#  [126.37829241   0.          26.03711851]
#  [686.96961895   0.         124.44253391]]
c = np.r_[weights * y, weights1 * y1]
print(c)
# [  83.66073785  383.70595203  159.12058215 1914.59065915 9981.85549321]
m, b1, b2 = np.linalg.lstsq(a, c, rcond=None)[0]
print(m, b1, b2)
# 10.012202998026055 14.841412336510793 24.941219918240172

编辑

如果你想要不同的斜率和一个截距，你可以这样做。可能最好掌握单斜率 2 截距情况的总体思路。看看数组a：你从权重和c 构造它，所以现在它是未加权的问题。您尝试找到这样的vector = [slope, intercept1, intercept2] 和a @ vector = c（通过最小化差异的平方和来尽可能地）。通过在a 中置零，我们使其可分离：矩阵的上半部分a 变化slope 和intercept1，a 的下半部分变化slope 和intercept2。类似于 vector = [slope1, slope2, intercept] 的 2 个斜坡案例。

first_column = np.r_[weights * x, [0] * x1.size]
second_column = np.r_[[0] * x.size, weights1 * x1]
third_column = np.r_[weights, weights1]