0

我想对 2D 数据点进行代数曲线拟合，但由于各种原因 - 一次在内存中存储大量样本数据是不可能的，并且迭代所有这些数据是一个昂贵的过程。

（这样做的原因是实际上我需要根据我从磁盘读取的千兆字节数据同时拟合数千条曲线，因此这很慢）。

请注意，多项式系数的数量将受到限制（可能是 5-10），因此完全拟合的可能性极小，但这没关系，因为我试图在具有大量随机噪声的数据中找到潜在模式。我了解如何使用遗传算法将曲线拟合到数据集，但这需要多次通过样本数据，因此对我的应用程序不实用。

有没有办法用单次数据拟合曲线，其中必须从样本到样本保持的状态是最小的？

我应该补充一点，数据的性质是这些点可能位于 X 轴上介于 0.0 和 1.0 之间的任何位置，但 Y 值将始终为 1.0 或 0.0。

因此，在 Java 中，我正在寻找具有以下接口的类：

public interface CurveFit {
   public void addData(double x, double y);
   public List<Double> getBestFit(); // Returns the polynomial coefficients
}

实现这一点的类必须不需要在其实例字段中保留太多数据，即使对于数百万个数据点也不超过一千字节。这意味着您不能只存储数据，因为您稍后会对其进行多次传递。

编辑：有人建议在一次通过中找到最佳曲线可能是不可能的，但是不需要最佳拟合，就像我们可以在一次通过中获得它一样接近。

一种方法的基本原理可能是，如果我们有一种从曲线开始的方法，然后有一种方法对其进行修改，使其在新数据点进入时稍微接近新数据点——实际上是一种梯度下降的形式。希望有足够的数据（并且数据会很丰富），我们得到一个很好的曲线。也许这会激发某人找到解决方案。

4

8 回答 8

2

是的，这是一个投影。为了

y = X beta + error

其中小写项是向量，X 是矩阵，你有解向量

\hat{beta} = inverse(X'X) X' y

根据OLS页面。您几乎不想直接计算它，而是使用 LR、QR 或 SVD 分解。统计文献中有大量参考文献。

如果您的问题只有一个参数（因此 x 也是一个向量），那么这将简化为 y 和 x 之间的叉积的总和。

于 2009-11-11T17:52:02.893 回答
2

如果您不介意得到一条直线“曲线”，那么对于任何数量的数据，您只需要六个变量。这是我即将出版的书中的源代码；我相信您可以弄清楚 DataPoint 类是如何工作的：

插值.h：

#ifndef __INTERPOLATION_H
#define __INTERPOLATION_H

#include "DataPoint.h"

class Interpolation
{
private:
  int m_count;
  double m_sumX;
  double m_sumXX;  /* sum of X*X */
  double m_sumXY;  /* sum of X*Y */
  double m_sumY;
  double m_sumYY;  /* sum of Y*Y */

public:
  Interpolation();

  void addData(const DataPoint& dp);

  double slope() const;
  double intercept() const;

  double interpolate(double x) const;
  double correlate() const;
};

#endif // __INTERPOLATION_H

插值.cpp：

#include <cmath>

#include "Interpolation.h"

Interpolation::Interpolation()
{
  m_count = 0;
  m_sumX = 0.0;
  m_sumXX = 0.0;
  m_sumXY = 0.0;
  m_sumY = 0.0;
  m_sumYY = 0.0;
}

void Interpolation::addData(const DataPoint& dp)
{
  m_count++;
  m_sumX += dp.getX();
  m_sumXX += dp.getX() * dp.getX();
  m_sumXY += dp.getX() * dp.getY();
  m_sumY += dp.getY();
  m_sumYY += dp.getY() * dp.getY();
}

double Interpolation::slope() const
{
  return (m_sumXY - (m_sumX * m_sumY / m_count)) /
    (m_sumXX - (m_sumX * m_sumX / m_count));
}

double Interpolation::intercept() const
{
  return (m_sumY / m_count) - slope() * (m_sumX / m_count);
}


double Interpolation::interpolate(double X) const
{
  return intercept() + slope() * X;
}


double Interpolation::correlate() const
{
  return m_sumXY / sqrt(m_sumXX * m_sumYY);
}
于 2009-11-13T06:22:22.890 回答
0

您是否限制多项式系数的数量（即拟合多项式中 x 的最大幂）？

如果不是，那么您不需要“最佳拟合”算法 - 您始终可以将 N 个数据点完全拟合到 N 个系数的多项式。

只需使用矩阵来求解 N 个未知数的 N 个联立方程（多项式的 N 个系数）。

如果您限制系数的最大数量，您的最大值是多少？

按照您的评论和编辑：

你想要的是一个低通滤波器来滤除噪声，而不是用多项式拟合噪声。

于 2009-11-13T19:31:13.373 回答
0

为什么不使用一些固定大小的环形缓冲区（例如最后 1000 个点），并使用标准的基于 QR 分解的最小二乘法拟合缓冲数据？一旦缓冲区填满，每次获得新点时，您都​​会替换最旧的点并重新拟合。这样你就有了一个有限的工作集，它仍然有一些数据局部性，没有实时流（无记忆）处理的所有挑战。

于 2009-11-13T19:44:21.970 回答
0

鉴于您的数据的性质：

这些点可能位于 X 轴上介于 0.0 和 1.0 之间的任何位置，但 Y 值将始终为 1.0 或 0.0。

那么你甚至不需要一次通过，因为这两行将完全通过每个点：

  • X = [0.0 ... 1.0]，Y = 0.0
  • X = [0.0 ... 1.0]，Y = 1.0

两条短线段，单位长度，每个点落在一条线上。

诚然，在单程中为任意点找到良好曲线拟合的算法很有趣，但是（根据您的问题），这不是您所需要的。

于 2009-11-13T20:07:23.693 回答
0

您需要超定线性系统的解决方案。流行的方法是正规方程（通常不推荐）、QR 分解和奇异值分解 (SVD)。维基百科有不错的解释，Trefethen 和 Bau非常好。您的选择：

  1. 通过正规方程的核外实现。这需要行数多于列数的产品A'A（A因此结果非常小）。该矩阵A完全由样本位置定义，因此您不必存储它，因此计算A'A相当便宜（如果您不需要为节点位置命中内存，则非常便宜）。计算完成后，您可以一次性A'A通过输入数据获得解决方案，但该方法可能不稳定。

  2. 实施核外 QR 分解。经典的 Gram-Schmidt 将是最快的，但您必须注意稳定性。

  3. 使用分布式内存在内核中执行此操作（如果您有可用的硬件）。像 PLAPACK 和 SCALAPACK 这样的库可以做到这一点，性能应该比 1 好得多。并行可扩展性不是很好，但如果它是一个你甚至可以考虑串行做的问题大小，那就没问题了。

  4. 使用迭代方法计算 SVD。根据您系统的光谱特性（可能在预处理之后），这可能会很快收敛并且不需要存储矩阵（在您的情况下，矩阵有 5-10 列，每列都是输入数据的大小。一个好的库因为这是SLEPc，你只需要找到 Vandermonde 矩阵与向量的乘积（所以你只需要存储样本位置）。这是非常可并行扩展的。

于 2009-11-13T20:13:17.097 回答
0

假设您不知道哪个点应该属于哪个曲线，那么像霍夫变换之类的东西可能会提供您需要的东西。

霍夫变换是一种允许您识别数据集中结构的技术。一种用途是用于计算机视觉，它可以轻松识别视野内的线条和边界。

这种情况的优点：

  • 每个点只需要考虑一次
  • 您不需要为每条候选线保留一个数据结构，只需一个（复杂的、多维的）结构
  • 每一行的处理都很简单
  • 您可以在任何时候停下来并输出一组好的匹配项
  • 您永远不会丢弃任何数据，因此它不依赖于任何偶然的引用位置
  • 您可以在准确性和内存要求之间进行权衡
  • 不限于完全匹配，但也会突出显示部分匹配。

一种方法

要找到三次拟合，您将构建一个 4 维霍夫空间，将每个数据点投影到其中。霍夫空间内的热点将为您提供通过这些点的三次参数。

于 2009-11-13T20:18:45.103 回答
0

我相信我根据此代码的修改版本找到了我自己的问题的答案。对于那些感兴趣的人，我的 Java 代码在这里。

于 2009-11-13T23:15:22.240 回答