statistics - 简单的多维曲线拟合

Question

我有一堆数据，一般是a,b,c,...,y

其中 y = f(a, b, c...)

其中大多数是三四个变量，并且有 10k - 10M 记录。我的一般假设是它们本质上是代数的，例如：

y = P1 a^E1 + P2 b^E2 + P3 c^E3

不幸的是，我上一次统计分析课是在 20 年前。获得 f 的良好近似值的最简单方法是什么？具有极小学习曲线的开源工具（即我可以在一小时左右得到一个体面的近似值的工具）将是理想的。谢谢！

Answer 1

如果它有用，这里有一个 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 项幂回归之类的操作。

Answer 2

数据拟合的基础包括假设解的一般形式，猜测常量的一些初始值，然后迭代以最小化猜测解的误差以找到特定解，通常在最小二乘意义上。

查看R或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

Answer 3

您知道要限制多项式的幂吗？

如果没有限制，那么您始终可以通过将其与具有 N 个系数的多项式进行匹配来获得 N 个点的精确匹配。为此，您将 N 个不同的点插入方程，产生 N 个方程和 N 个未知数（系数），然后您可以使用简单的高中代数或矩阵来求解未知数。

Answer 4

我花了一个多星期的时间试图做同样的事情。我尝试了一大堆优化方法来微调系数，但基本上没有成功，然后我发现有一个封闭形式的解决方案，而且效果很好。

免责声明：我试图以固定的最大数量级拟合数据。如果您的 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 = c1 a^2 b^2 + c2 a^2 b + c3 a^2 + c4 a b^2 + c5 a b + c6 a + c7 b^2形式的二阶多项式+ c8 b + 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	一个N^2	aN*bN^2	aN*bN	一个	bN^2	氮化硼	1

您的 Y 矩阵很简单：

y1
y2
...
yN

然后你做 B = (X ^T X) ^-1 X ^T Y 然后 B 将等于

c1
c2
c3
c4
c5
c6
c7
c8
c9

请注意，系数的总数将是 (o + 1) ^V，其中 o 是多项式的阶，V 是变量的数量，因此它增长得非常快。

如果您使用的是好的矩阵代码，那么我相信运行时复杂度将为 O(((o+1) ^V ) ³ + ((o + 1) ^V ) ² N)，其中 V 是变量的数量，o 是多项式的阶数，N 是您拥有的数据输入数。最初这听起来很糟糕，但在大多数情况下，o 和 V 可能不会很高，所以这只是关于 N 的线性关系。

请注意，如果您正在编写自己的矩阵代码，那么确保您的反演代码使用类似于LU 分解的东西很重要。如果你使用一种朴素的反转方法（就像我一开始所做的那样），那么 ((o+1) ^V ) ³就会变成 ((o+1) ^V )!，这就更糟了。在我做出这个改变之前，我预测我的 5 阶 3 变量多项式将需要大约 400 谷歌千年才能完成。使用LU分解后，大约需要7秒。

另一个免责声明

这种方法要求 (X ^T X) 不是奇异矩阵（换句话说，它可以反转）。我的线性代数有点粗糙，所以我不知道会发生这种情况的所有情况，但我知道当输入变量之间存在完美的多重共线性时会发生这种情况。这意味着一个变量只是另一个因子乘以一个常数（例如，一个输入是完成一个项目的小时数，另一个是完成一个项目的美元，但美元只是基于小时费率乘以小时数）。

好消息是，当存在完美的多重共线性时，你会知道的。当你反转矩阵时，你最终会被零除或其他东西。

更大的问题是当你有不完美的多重共线性时。那时您有两个密切相关但不完全相关的变量（例如温度和高度，或速度和马赫数）。在这些情况下，这种方法在理论上仍然有效，但它变得非常敏感，以至于小的浮点错误可能会导致结果偏离。

然而，在我的观察中，结果要么非常好，要么非常糟糕，所以你可以为你的均方误差设置一些阈值，如果超过了这个阈值，那么就说“无法计算多项式”。

Answer 5

如果您猜测 f,[*] 的形式，您需要一个最小化器来找到最佳参数。Scottie T 建议的工具以及ROOT和许多其他工具都可以使用。

如果您不知道 f 可能采用哪种形式，那么您确实陷入了大麻烦。

---

[*] 也就是说，你知道

f = f(x,y,z,w,...;p1,p2,p3...)

其中ps 是参数，坐标是x, y...

Answer 6

简短的回答：这不是那么简单。考虑一种关于数据子集的非参数方法。

您需要决定两个主要问题 (1) 您是否真的关心函数的参数，即您的 P1、E1、...，或者您是否可以只估计平均函数 (2)？真的需要估计所有数据的函数吗？

我要提到的第一件事是您指定的函数是非线性的（在要估计的参数中），所以普通的最小二乘法不起作用。假设您指定了一个线性函数。10M 值仍然存在问题。可以使用 QR 分解以有效的方式执行线性回归，但您仍然需要 O(p * n^2) 算法，其中 p 是您尝试估计的参数数量。如果你想估计非线性平均函数，它会变得更糟。

在如此大的数据集中估计任何东西的唯一方法是使用一个子集来执行估计。基本上，您随机选择一个子集并使用它来估计函数。

如果您不关心参数值，而只想估计平均函数，那么使用非参数估计技术可能会更好。

希望这会有所帮助。

莱夫