简单的多维曲线拟合

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

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

\n

免责声明：我试图用固定的最大数量级来拟合数据。如果您的 E1、E2 等值没有限制，那么这对您不起作用。

\n

现在我已经花时间学习这些东西了，我实际上发现，如果我理解某些答案，它们会给出很好的提示。距离上一次统计和线性代数课也已经有一段时间了。

\n

因此，如果还有其他人缺乏线性代数知识，这就是我所做的。

\n

尽管这不是您想要拟合的线性函数，但事实证明这仍然是一个线性回归问题。维基百科有一篇关于线性回归的非常好的文章。我建议慢慢阅读：https://en.wikipedia.org/wiki/Linear_regression# :~:text =In%20statistics%2C%20linear%20regression%20is,as%20dependent%20and%20independent%20variables）。它还链接了许多其他相关的优秀文章。

\n

如果您不知道如何使用矩阵解决简单（单变量）线性回归问题，请花一些时间学习如何做到这一点。

\n

一旦您学习了如何进行简单线性回归，就可以尝试一些多元线性回归。基本上，要进行多变量线性回归，您需要创建一个 X 矩阵，其中每个输入数据项都有一行，每行包含该数据条目的所有变量值（加上最后一列中使用的 1）多项式末尾的常量值（称为截距））。然后创建一个 Y 矩阵，该矩阵是单列，每个数据项占一行。然后你求解 B = (X ^T X) ^-1 X ^T Y。然后 B 就成为多项式的所有系数。

\n

对于多变量多项式回归，其想法相同，现在您有一个巨大的多变量线性回归，其中每个回归量（您正在执行回归的变量）都是您的巨型多项式表达式的系数。

\n

因此，如果您的输入数据如下所示：

\n

\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n \n

输入	输出
a1、b1、	y1
a2、b2、	y2
...	...
aN, bN,	yN

\n

\n

您想要拟合形式为 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 矩阵将如下所示：

\n

\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n \n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n \n\n\n\n\n\n\n\n\n\n\n

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

\n

\n

你的 Y 矩阵很简单：

\n

\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n

y1
y2
...
yN

\n

\n

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

\n

\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n \n\n\n\n\n\n\n\n\n\n\n

c1
c2
c3
c4
c5
c6
c7
c8
c9

\n

\n

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

\n

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

\n

请注意，如果您正在编写自己的矩阵代码，那么确保您的反演代码使用诸如LU 分解之类的东西非常重要。如果您使用 na\xc3\xafve 反转方法（就像我一开始所做的那样），那么 ((o+1) ^V ) ³就会变成 ((o+1) ^V )!，这更糟糕。在进行更改之前，我预测我的 5 阶 3 变量多项式大约需要 400 google millennia 才能完成。使用LU分解后，大约需要7秒。

\n

另一项免责声明

\n

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

\n

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

\n

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

\n

然而，根据我的观察，结果要么非常好，要么非常糟糕，因此您可以为均方误差设置一些阈值，如果超过该阈值，则说“无法计算多项式”。

\n

Answer 3

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

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