我有一堆数据，通常是以这种形式呈现：
a, b, c, ..., y
其中y = f(a, b, c...)
它们大多数都是三个或四个变量，并且有10k-10M的记录。我的基本假设是它们具有代数性质，例如:
y = P1 a^E1 + P2 b^E2 + P3 c^E3
不幸的是，我上一次统计分析课已经是20年前的事了。有什么简单的方法可以很好地近似f吗？最好是开源工具且学习曲线非常低（即我可以在一个小时左右获得不错的近似值）。谢谢！
简单的多维曲线拟合
23
- user64258
3
6关于标题，“多维曲线拟合”有何简单之处？ :-) - Chris Ballance
http://www.prz.rzeszow.pl/~janand/Theory_of_LSM.pdf - jfs
1正交距离回归可以用于这个问题。 - jfs
6个回答
14
如果有用的话，这里是一个Numpy/Scipy (Python)模板来实现你想要的:
from numpy import array
from scipy.optimize import leastsq
def __residual(params, y, a, b, c):
p0, e0, p1, e1, p2, e2 = params
return p0 * a ** e0 + p1 * b ** e1 + p2 * c ** e2 - y
# load a, b, c
# guess initial values for p0, e0, p1, e1, p2, e2
p_opt = leastsq(__residual, array([p0, e0, p1, e1, p2, e2]), args=(y, a, b, c))
print 'y = %f a^%f + %f b^%f %f c^%f' % map(float, p_opt)
如果你真的想理解正在发生的事情，那么你需要投入时间去学习一些工具或者编程环境-我真的认为没有任何绕过这个步骤的方法。人们通常不会专门为做类似于三项幂回归之类的事情编写专用工具。
- David
2
1如果a，b，c没有无限精度（最小二乘法假定坐标有无限精度），scipy.odr（正交距离回归）可能会很有用。 - jfs
这个函数肯定需要一些样本输出来最小化，即给定一组a、b、c值的一些样本y值？ - Brendan
5
我花了一周多的时间试图做基本相同的事情。我尝试了很多优化方法来微调系数，但基本上没有成功，然后我发现有一个闭式解决方案，它非常有效。
免责声明：我试图拟合数据，并限制最大数量级。如果您的E1、E2等值没有限制，那么这对您无效。
现在我花时间学习这些东西，我实际上看到，如果我理解了它们，其中一些答案会给出好的提示。这也是自从我的上一次统计和线性代数课以来过了一段时间。
因此，如果还有其他缺乏线性代数知识的人，请看看我做了什么。
尽管您正在尝试拟合的不是线性函数，但它仍然是一个线性回归问题。维基百科有一篇关于线性回归的非常好的文章。我建议您慢慢阅读：https://en.wikipedia.org/wiki/Linear_regression#:~:text=In%20statistics%2C%20linear%20regression%20is,as%20dependent%20and%20independent%20variables）。它还链接了许多其他相关的好文章。
如果您不知道如何使用矩阵做一个简单（单变量）的线性回归问题，请花些时间学习如何做。
一旦学会了简单的线性回归，就可以尝试一些多变量线性回归。基本上，要进行多元线性回归，您需要创建一个X矩阵，其中每个输入数据项都有一行，并且每行包含该数据条目的所有变量值（最后一列为1，用于表示多项式末端的常数（称为拦截器））。然后，您创建一个Y矩阵，它是一个单列，每个数据项都有一行。然后解决B =（X T X） -1 X T Y。 B然后成为您的多项式的所有系数。
对于多元多项式回归，想法相同，现在您拥有一个巨大的多元线性回归，其中每个回归器（要进行回归的变量）都是巨大多项式表达式的系数。
如果您的输入数据如下：
您的Y矩阵简单地表示为：
然后你做B = (XTX)-1XTY，然后B将等于：
请注意，系数的总数将是(o + 1)V，其中o是多项式的阶数，V是变量的数量，因此它增长得非常快。
如果您使用良好的矩阵代码，那么我相信运行时复杂度将是O(((o+1)V)3 + ((o+1)V)2N)，其中V是变量数，o是多项式的次数，N是输入数据的数量。一开始听起来很糟糕，但在大多数情况下，o和V可能不会很高，因此它就变成了相对于N的线性复杂度。
请注意，如果您正在编写自己的矩阵代码，则重要的是确保您的求逆代码使用类似LU分解的东西。如果您使用天真的求逆方法（就像我一开始做的那样），那么((o+1)V)3变成了((o+1)V)！，这更糟糕。在进行修改之前，我预计我的五阶三元多项式需要大约400亿年才能完成。使用LU分解后，只需要大约7秒钟。
另一个免责声明
这种方法要求(XTX)不是奇异矩阵（换句话说，可以被反转）。我的线性代数有点生疏，所以我不知道所有会出现这种情况的案例，但我知道当输入变量之间存在完美的多重共线性时会发生这种情况。这意味着一个变量只是另一个因素乘以一个常数（例如，一个输入是完成项目所需的小时数，另一个是完成项目所需的美元，但美元只是基于每小时费率乘以小时数）。
好消息是，当存在完美的多重共线性时，你会知道。当你反转矩阵时，你会得到除以零或其他错误信息。
更大的问题是当你有不完美的多重共线性时。这是指你有两个密切相关但不完全相关的变量（如温度和海拔，或速度和马赫数）。在这些情况下，这种方法在理论上仍然有效，但它变得非常敏感，小的浮点误差就可能导致结果完全偏离预期。
根据我的观察，结果要么很好，要么很差，因此你可以设置一些均方误差的阈值，如果超过了这个阈值，那么就说“无法计算多项式”。
免责声明：我试图拟合数据，并限制最大数量级。如果您的E1、E2等值没有限制，那么这对您无效。
现在我花时间学习这些东西，我实际上看到，如果我理解了它们，其中一些答案会给出好的提示。这也是自从我的上一次统计和线性代数课以来过了一段时间。
因此，如果还有其他缺乏线性代数知识的人，请看看我做了什么。
尽管您正在尝试拟合的不是线性函数，但它仍然是一个线性回归问题。维基百科有一篇关于线性回归的非常好的文章。我建议您慢慢阅读：https://en.wikipedia.org/wiki/Linear_regression#:~:text=In%20statistics%2C%20linear%20regression%20is,as%20dependent%20and%20independent%20variables）。它还链接了许多其他相关的好文章。
如果您不知道如何使用矩阵做一个简单（单变量）的线性回归问题，请花些时间学习如何做。
一旦学会了简单的线性回归，就可以尝试一些多变量线性回归。基本上，要进行多元线性回归，您需要创建一个X矩阵，其中每个输入数据项都有一行，并且每行包含该数据条目的所有变量值（最后一列为1，用于表示多项式末端的常数（称为拦截器））。然后，您创建一个Y矩阵，它是一个单列，每个数据项都有一行。然后解决B =（X T X） -1 X T Y。 B然后成为您的多项式的所有系数。
对于多元多项式回归，想法相同，现在您拥有一个巨大的多元线性回归，其中每个回归器（要进行回归的变量）都是巨大多项式表达式的系数。
如果您的输入数据如下：
| 输入 | 输出 |
|---|---|
| a1，b1， | y1 |
| a2，b2， | y2 |
| ... | ... |
| aN，bN， | yN |
如果你想要拟合一个二次多项式，形如 y = c1a^2b^2 + c2a^2b + c3a^2 + c4ab^2 + c5ab + c6a + c7b^2 + c8b + c9，那么你的 X 矩阵将会是：
| a1^2*b1^2 | a1^2*b1 | a1^2 | a1*b1^2 | a1*b1 | a1 | b1^2 | b1 | 1 |
| a2^2*b2^2 | a2^2*b2 | a2^2 | a2*b1^2 | a2*b2 | a2 | b2^2 | b2 | 1 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... |
| aN^2*bN^2 | aN^2*bN | aN^2 | aN*bN^2 | aN*bN | aN | bN^2 | bN | 1 |
| y1 |
| y2 |
| ... |
| yN |
| c1 |
| c2 |
| c3 |
| c4 |
| c5 |
| c6 |
| c7 |
| c8 |
| c9 |
如果您使用良好的矩阵代码，那么我相信运行时复杂度将是O(((o+1)V)3 + ((o+1)V)2N)，其中V是变量数，o是多项式的次数，N是输入数据的数量。一开始听起来很糟糕，但在大多数情况下，o和V可能不会很高，因此它就变成了相对于N的线性复杂度。
请注意，如果您正在编写自己的矩阵代码，则重要的是确保您的求逆代码使用类似LU分解的东西。如果您使用天真的求逆方法（就像我一开始做的那样），那么((o+1)V)3变成了((o+1)V)！，这更糟糕。在进行修改之前，我预计我的五阶三元多项式需要大约400亿年才能完成。使用LU分解后，只需要大约7秒钟。
另一个免责声明
这种方法要求(XTX)不是奇异矩阵（换句话说，可以被反转）。我的线性代数有点生疏，所以我不知道所有会出现这种情况的案例，但我知道当输入变量之间存在完美的多重共线性时会发生这种情况。这意味着一个变量只是另一个因素乘以一个常数（例如，一个输入是完成项目所需的小时数，另一个是完成项目所需的美元，但美元只是基于每小时费率乘以小时数）。
好消息是，当存在完美的多重共线性时，你会知道。当你反转矩阵时，你会得到除以零或其他错误信息。
更大的问题是当你有不完美的多重共线性时。这是指你有两个密切相关但不完全相关的变量（如温度和海拔，或速度和马赫数）。在这些情况下，这种方法在理论上仍然有效，但它变得非常敏感，小的浮点误差就可能导致结果完全偏离预期。
根据我的观察，结果要么很好，要么很差，因此你可以设置一些均方误差的阈值，如果超过了这个阈值，那么就说“无法计算多项式”。
- NateW
3
数据拟合的基本原理是假设解的一般形式，猜测常数的初始值，然后迭代以最小化猜测解的误差，找到一个特定的解，通常采用最小二乘法。
请查看开源工具R或Octave。它们都能进行最小二乘分析，并有多个教程可供搜索。 编辑： 用于估计二次多项式系数的Octave代码
在我的机器上，我得到：
请查看开源工具R或Octave。它们都能进行最小二乘分析，并有多个教程可供搜索。 编辑： 用于估计二次多项式系数的Octave代码
x = 0:0.1:10;
y = 5.*x.^2 + 4.*x + 3;
% Add noise to y data
y = y + randn(size(y))*0.1;
% Estimate coefficients of polynomial
p = polyfit(x,y,2)
在我的机器上，我得到：
ans =
5.0886 3.9050 2.9577
- Scottie T
6
谢谢，我明白了...这就是为什么我说“学习曲线非常低”！这两种语言都是出色的通用统计语言，但是它们的学习曲线相对较陡峭（依我之见）。 - user64258
我明白了。我认为，对于简单函数而言，使用任何一种工具都不应该需要太长时间来适应，甚至可以使用Python或Perl来完成这个任务。 - Scottie T
我认为它们相对简单（我在问题中添加了详细信息），而且我已经在谷歌上花了大约一个小时，这就是为什么我来这里的原因;-) - user64258
不幸的是，polyfit 只适用于单值函数 f(x)。OP 特别提到了非线性多维拟合，我怀疑 Octave 默认情况下不支持此功能。 - Kena
我认为你不可能会比那个Octave代码更简单（或者用几乎相同的语法使用Python中的Numpy/Scipy-请参考http://www.scipy.org）。 - David Z
但它不是多维的，对于一维情况才适用。 - tkarahan
1
你知道你想要将多项式限制到哪个幂次吗？
如果没有限制，那么你总是可以通过匹配具有N个系数的多项式来获得N个点的精确匹配。为此，您将N个不同的点插入方程式中，得到N个方程和N个未知数（系数），然后您可以使用简单的高中代数或矩阵来解决这些未知数。
如果没有限制，那么你总是可以通过匹配具有N个系数的多项式来获得N个点的精确匹配。为此，您将N个不同的点插入方程式中，得到N个方程和N个未知数（系数），然后您可以使用简单的高中代数或矩阵来解决这些未知数。
- mbeckish
1
+1，我在某个地方读到稀疏网格数据可以使用比规则网格数据更少的节点实现相同的多项式精度。您知道这是如何实现的吗？ - owari
0
如果你对 f 的形式有猜测，[*] 你需要一个最小化器来找到最优参数。Scottie T 建议的工具 可以使用，ROOT 也可以，还有许多其他工具。
如果你不知道 f 可能采取什么形式，那你真的很麻烦。
[*] 意思是，你知道
f = f(x,y,z,w,...;p1,p2,p3...)
其中p是参数，而坐标是x、y等等。
- dmckee
0
简短回答：这并不简单。考虑在数据子集上使用非参数方法。
有两个主要问题需要决定：（1）您是否真的关心函数的参数，即P1、E1等，或者只是估计平均函数；（2）您是否真的需要在所有数据上估计函数？
首先要提到的是，您指定的函数是非线性的（在要估计的参数中），因此普通最小二乘法无法工作。假设您指定了一个线性函数。您仍然会遇到10M个值的问题。可以使用QR分解以有效的方式执行线性回归，但您仍然面临O(p * n^2)算法的问题，其中p是您尝试估计的参数数量。如果您想估计非线性平均函数，情况会变得更糟。
您唯一能够在如此大的数据集中估计任何东西的方法是使用子集来执行估计。基本上，您随机选择一个子集并使用它来估计函数。
如果您不关心参数值，只想估计平均函数，那么使用非参数估计技术可能会更好。
希望这能有所帮助。
Leif
有两个主要问题需要决定：（1）您是否真的关心函数的参数，即P1、E1等，或者只是估计平均函数；（2）您是否真的需要在所有数据上估计函数？
首先要提到的是，您指定的函数是非线性的（在要估计的参数中），因此普通最小二乘法无法工作。假设您指定了一个线性函数。您仍然会遇到10M个值的问题。可以使用QR分解以有效的方式执行线性回归，但您仍然面临O(p * n^2)算法的问题，其中p是您尝试估计的参数数量。如果您想估计非线性平均函数，情况会变得更糟。
您唯一能够在如此大的数据集中估计任何东西的方法是使用子集来执行估计。基本上，您随机选择一个子集并使用它来估计函数。
如果您不关心参数值，只想估计平均函数，那么使用非参数估计技术可能会更好。
希望这能有所帮助。
Leif
- leif
网页内容由stack overflow 提供, 点击上面的可以查看英文原文，
原文链接
原文链接