python - 使用 numpy polyfit 在 python 中使用具有多个变量的 polyfit 将数据拟合到曲线

Question

有没有办法计算两个变量的多项式模型的参数。它们是独立的，因此:

z = a + bx + cx^2 + dy + ey^2

有人告诉我，您可以为此使用 numpy.polyfit，但它只能支持两个变量，而不是我需要的三个变量。我的数据目前存储在三个 numpy 数组中，这样数组中每条数据的索引都与其他变量中该索引处的数据相关联。即

Y = [1 2 3 4 5]
X = [3 5 7 9 11]
Z = [2 4 6 2 6]

1 与 3 和 2 关联；2 与 5 和 4 关联，依此类推。我该如何解决这个问题？

Answer 1

polyfit 假定一个变量。但是，您要做的是求解一般线性方程组。

用线性代数表达问题

你有一个等式:

z = a + bx + cx^2 + dy + ey^2

以及观察 x、y 和 z 的 5 个点。这给了我们 5 个等式:

z1 = a + bx1 + cx1^2 + dy1 + ey1^2
z2 = a + bx2 + cx2^2 + dy2 + ey2^2
z3 = a + bx3 + cx3^2 + dy3 + ey3^2
z4 = a + bx4 + cx4^2 + dy4 + ey4^2
z5 = a + bx5 + cx5^2 + dy5 + ey5^2

最容易将其视为线性代数问题。我们可以将方程组重写为矩阵乘法:

|z1|   |1  x1  x1^2  y1  y1^2|   |a|
|z2|   |1  x2  x2^2  y2  y2^2|   |b|
|z3| = |1  x3  x3^2  y3  y3^2| x |c|
|z4|   |1  x4  x4^2  y4  y4^2|   |d|
|z5|   |1  x5  x5^2  y5  y5^2|   |e|

我们知道 x1..5、y1..5 和 z1..5，但是 a、b、 c、d、e 是未知数。我们分别称这些矩阵为 B、A 和 x:

B = A x

我们想要求解 x，它是我们的 a、b、c、d、e 参数的矩阵。

---

奇异矩阵

但是，我们有一个问题。您为 x、y 和 z 提供的确切数字会生成 A 的奇异矩阵。换句话说，行不是独立的。我们实际上已经将同一个等式输入了两次。事实上，在这种情况下，我们只有 3 个等式。其他两个只是前三个的组合。

无法使用您提供的 X、Y、Z 数据求解方程组。

考虑到这一点，让我们将问题更改为使用 5 个随机 x、y、z 点。

精确解

在这种特定情况下，我们正好有 5 个未知数和 5 个方程。因此，我们可以准确地解决这个问题(例如使用 np.linalg.solve)。这被称为“偶数确定问题”。

import numpy as np
np.random.seed(1)

# Each array with have 5 random points
x, y, z = np.random.random((3, 5))

# A is going to look like:
#  |1  x1  x1^2  y1  y1^2|
#  |1  x2  x2^2  y2  y2^2|
#  |1  x3  x3^2  y3  y3^2|
#  |1  x4  x4^2  y4  y4^2|
#  |1  x5  x5^2  y5  y5^2|
A = np.column_stack([np.ones(len(x)), x, x**2, y, y**2])

# And "B" will just be our "z" array
B = z

# Now we can solve the system of equations Ax = B
a, b, c, d, e = np.linalg.solve(A, B)

---

超过 5 个观察结果

但是，假设我们有 10 个或 100 个观测值。在这种情况下，我们就会遇到一个超定问题。我们无法精确求解，而是需要使用最小二乘拟合。

在那种情况下，您仍然会用矩阵乘法来表达事物并求解 Ax = B。但是，A 不会是方阵。因此，我们需要使用不同的工具来解决问题。对于 numpy，它是 np.linalg.lstsq 而不是 np.linalg.solve:

我将很快对此进行详细说明(可能会有点)，但目前:

import numpy as np
np.random.seed(1)

# Each array with have 20 random points this time
x, y, z = np.random.random((3, 20))

# We're solving Ax = B
A = np.column_stack([np.ones(len(x)), x, x**2, y, y**2])
B = z

# Solve the system of equations.
result, _, _, _ = np.linalg.lstsq(A, B)
a, b, c, d, e = result

删除条款

如果您想从等式中删除 bx 和 dy 项，您可以将它们从 A 中删除:

import numpy as np
np.random.seed(1)

x, y, z = np.random.random((3, 20))

# Note that we've remove the `x` and `y` terms.
# We're solving `z = a + cx^2 + ey^2`
A = np.column_stack([np.ones(len(x)), x**2, y**2])
B = z

# Solve the system of equations.
result, _, _, _ = np.linalg.lstsq(A, B)
a, c, e = result