我有一堆数据，一般是表格 a, b, c, ..., 是的
其中 y = f(a, b, c...)
其中大多数是三个和四个变量，并且有 10k - 10M 条记录。我的一般假设是它们本质上是代数的，例如：
y = P1 a^E1 + P2 b^E2 + P3 c^E3
不幸的是，我上一次统计分析课是在 20 年前。获得 f 的良好近似值的最简单方法是什么？具有极小学习曲线的开源工具（即我可以在一小时左右得到一个像样的近似值的东西）将是理想的。谢谢！
【问题讨论】：
标签： statistics regression best-fit-curve
我花了一个多星期的时间试图做同样的事情。我尝试了一大堆优化的东西来微调系数，但基本上没有成功，然后我发现有一个封闭形式的解决方案，效果很好。
免责声明：我试图以固定的最大数量级拟合数据。如果您的 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 = (XTX)-1XTY。 B 然后成为多项式的所有系数。
对于多变量多项式回归，它的想法是一样的，只是现在你有一个巨大的多变量线性回归，其中每个回归量（你正在做回归的变量）是你的巨型多项式表达式的一个系数。
因此，如果您的输入数据如下所示：
| Inputs | Output |
|---|---|
| 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 |
你的 Y 矩阵很简单：
| y1 |
| y2 |
| ... |
| yN |
然后你做 B = (XTX)-1XTY 然后 B 等于
| c1 |
| c2 |
| c3 |
| c4 |
| c5 |
| c6 |
| c7 |
| c8 |
| c9 |
请注意，系数的总数将是 (o + 1)V，其中 o 是多项式的阶，V 是变量的数量，因此它增长得非常快。
如果你使用好的矩阵代码，那么我相信运行时复杂度将是 O(((o+1)V)3 + ((o + 1 )V)2N)，其中 V 是变量数，o 是多项式的阶数，N 是您拥有的数据输入数。起初这听起来很糟糕，但在大多数情况下，o 和 V 可能不会很高，所以这只是关于 N 的线性。
请注意，如果您正在编写自己的矩阵代码，那么确保您的反演代码使用类似LU decomposition 的内容很重要。如果您使用天真的反转方法（就像我一开始那样），那么 ((o+1)V)3 将变为 ((o+1)V )!，情况更糟。在我做出这个改变之前，我预测我的 5 阶 3 变量多项式将需要大约 400 谷歌千年才能完成。使用LU分解后，大约需要7秒。
这种方法要求 (XTX) 不是奇异矩阵（换句话说，它可以反转）。我的线性代数有点粗糙，所以我不知道会发生这种情况的所有情况，但我知道当输入变量之间存在完美的多重共线性时会发生这种情况。这意味着一个变量只是另一个因子乘以一个常数（例如，一个输入是完成一个项目的小时数，另一个是完成一个项目的美元，但美元只是基于小时费率乘以小时数)。
好消息是，当存在完美的多重共线性时，您就会知道。当你反转矩阵时，你最终会被零除或其他东西。
更大的问题是当你有不完美的多重共线性时。那时您有两个密切相关但不完全相关的变量（例如温度和高度，或速度和马赫数）。在这些情况下，这种方法在理论上仍然有效，但它变得非常敏感，以至于小的浮点错误可能会导致结果偏离。
然而，在我的观察中，结果要么非常好，要么非常糟糕，因此您可以为均方误差设置一些阈值，如果超过该阈值，则说“无法计算多项式”。
【讨论】：
如果它有用，这里有一个 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)
但是，如果您真的想了解正在发生的事情，您将不得不投入时间来扩展某些工具或编程环境的学习曲线 - 我真的认为没有任何方法可以解决这个问题。人们通常不会编写专门的工具来专门执行 3 项幂回归等操作。
【讨论】：
简答：没那么简单。考虑一种关于数据子集的非参数方法。
您需要决定 2 个主要问题 (1) 您是否真的关心函数的参数，即您的 P1、E1、...，或者您是否可以只估计平均函数 (2 ) 你真的需要估计所有数据的函数吗？
我要提到的第一件事是您指定的函数是非线性的（在要估计的参数中），所以普通的最小二乘法不起作用。假设您指定了一个线性函数。 10M 值仍然存在问题。线性回归可以使用 QR 分解以有效的方式执行，但您仍然需要 O(p * n^2) 算法，其中 p 是您尝试估计的参数数量。如果你想估计非线性平均函数，它会变得更糟。
要在如此庞大的数据集中估计任何内容，唯一的方法是使用子集来执行估计。基本上，您随机选择一个子集并使用它来估计函数。
如果您不关心参数值，而只想估计均值函数，则最好使用非参数估计技术。
希望这会有所帮助。
leif
【讨论】：
您知道要限制多项式的幂吗？
如果没有限制，那么您始终可以通过将 N 个点与具有 N 个系数的多项式进行匹配来获得 N 个点的精确匹配。为此，您将 N 个不同的点代入方程，产生 N 个方程和 N 个未知数（系数），然后您可以使用简单的高中代数或矩阵来求解未知数。
【讨论】：
数据拟合的基础包括假设解的一般形式，猜测常量的一些初始值，然后迭代以最小化猜测解的误差以找到特定解，通常在最小二乘意义上。
查看R 或Octave 以获取开源工具。它们都能够进行最小二乘分析，几个教程只需 Google 搜索即可。
编辑：用于估计二阶多项式系数的八度代码
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
【讨论】：
如果您猜测 f,[*] 的形式，您需要一个最小化器来找到最佳参数。 The tools Scottie T suggests 可以工作，ROOT 和许多其他人也可以。
如果您不知道 f 可能采用哪种形式，那么您确实有很大的麻烦。
[*]也就是说，你知道的
f = f(x,y,z,w,...;p1,p2,p3...)
其中ps 是参数，坐标是x、y...
【讨论】：