用两条线拟合曲线的二次算法

Question

问题是找到实值二维曲线（由点集给出）与由两条线组成的折线的最佳拟合。

蛮力方法是为曲线的每个点找到“左”和“右”线性拟合，然后选择误差最小的对。我可以在遍历曲线点的同时增量计算两个线性拟合，但我找不到增量计算误差的方法。因此，这种方法会产生二次复杂度。

问题是是否有一种算法可以提供次二次复杂度？

第二个问题是这种算法是否有一个方便的 C++ 库？

---

编辑 对于单行拟合，有公式：

m = (Σxiyi - ΣxiΣyi/N) / (Σxi2 - (Σxi)2/N)
b = Σyi/N - m * Σxi/N

其中m 是斜率，b 是直线的偏移量。 拥有这样的拟合误差公式可以最好地解决问题。

Answer 1

免责声明：我不想弄清楚如何在 C++ 中执行此操作，因此我将使用 Python (numpy) 表示法。这些概念是完全可以转换的，因此您可以毫不费力地翻译回您选择的语言。

假设您有一对数组，x 和 y，包含数据点，并且 x 是单调递增的。还假设您将始终选择一个分区点，该分区点在每个分区中至少留下两个元素，因此方程是可解的。

现在您可以计算一些相关数量：

N = len(x)

sum_x_left = x[0]
sum_x2_left = x[0] * x[0]
sum_y_left = y[0]
sum_y2_left = y[0] * y[0]
sum_xy_left = x[0] * y[0]

sum_x_right = x[1:].sum()
sum_x2_right = (x[1:] * x[1:]).sum()
sum_y_right = y[1:].sum()
sum_y2_right = (y[1:] * y[1:]).sum()
sum_xy_right = (x[1:] * y[1:]).sum()

我们需要这些量（它们是O(N) 来初始化）的原因是您可以直接使用它们来计算一些众所周知的线性回归参数公式。例如，y = m * x + b 的最佳 m 和 b 由下式给出

μ_x = Σx_i/N μ_y = Σy_i/N m = Σ(x_i - μ_x)(y_i - μ_y) / Σ(xi - μ_x)² b = μ_y - m * μ_x

误差平方和由下式给出

e = Σ(y_i - m * x_i - b)²

这些可以使用简单的代数展开为：

m = (Σx_iy_i - Σx_iΣy_i/N) / (Σx_i ² - (Σx_i)²/N) b = Σy_i/N - m * Σx_i/N e = Σy_i² + m² * Σx_i² + N * b² - 2 * m * Σx_iy_i - 2 * b * Σy_i + 2 * m * b * Σx_i

因此，您可以遍历所有可能性并记录最小的e：

for p in range(1, N - 3):
    # shift sums: O(1)
    sum_x_left += x[p]
    sum_x2_left += x[p] * x[p]
    sum_y_left += y[p]
    sum_y2_left += y[p] * y[p]
    sum_xy_left += x[p] * y[p]

    sum_x_right -= x[p]
    sum_x2_right -= x[p] * x[p]
    sum_y_right -= y[p]
    sum_y2_right -= y[p] * y[p]
    sum_xy_right -= x[p] * y[p]

    # compute err: O(1)
    n_left = p + 1
    slope_left = (sum_xy_left - sum_x_left * sum_y_left * n_left) / (sum_x2_left - sum_x_left * sum_x_left / n_left)
    intercept_left = sum_y_left / n_left - slope_left * sum_x_left / n_left
    err_left = sum_y2_left + slope_left * slope_left * sum_x2_left + n_left * intercept_left * intercept_left - 2 * (slope_left * sum_xy_left + intercept_left * sum_y_left - slope_left * intercept_left * sum_x_left)

    n_right = N - n_left
    slope_right = (sum_xy_right - sum_x_right * sum_y_right * n_right) / (sum_x2_right - sum_x_right * sum_x_right / n_right)
    intercept_right = sum_y_right / n_right - slope_right * sum_x_right / n_right
    err_right = sum_y2_right + slope_right * slope_right * sum_x2_right + n_right * intercept_right * intercept_right - 2 * (slope_right * sum_xy_right + intercept_right * sum_y_right - slope_right * intercept_right * sum_x_right)

    err = err_left + err_right
    if p == 1 || err < err_min
        err_min = err
        n_min_left = n_left
        n_min_right = n_right
        slope_min_left = slope_left
        slope_min_right = slope_right
        intercept_min_left = intercept_left
        intercept_min_right = intercept_right

您可能还可以进行其他简化，但这足以拥有O(n) 算法。

Answer 2

如果它有助于这是我用于这种事情的C代码。它对疯狂的物理学家说的很少。

首先关闭，公式。如果通过某些点拟合y ^：x-> a * x + b线，则错误给出：

E = Sum{ sqr(y[i]-y^(x[i])) }/ N = Vy - Cxy*Cxy/Vx
where 
Vx is the variance of the xs
Vy that of the ys 
Cxy the covariance of the as and the ys

下面的代码使用包含均值，差异，协方差和计数的结构。

函数moms_acc_pt（）添加新点时更新这些。函数moms_line（）返回该行的A和B，以及如上所述的错误。返回的fmax（0，）是在近乎完美的符合胜度的情况下，舍入错误可以发送（数学上非负）结果负。

虽然可以具有从瞬间删除一个点的函数，但我发现决定通过乘坐副本来添加一个点，累积副本中的点，恢复副本中的点，恢复副本的函数更容易处理。保持副本的副本是最适合的点，以及其他

的原件

typedef struct
{   int n;      // number points
    double  xbar,ybar;  // means of x,y
    double  Vx, Vy;     // variances of x,y
    double  Cxy;        // covariance of x,y
}   momentsT;

// update the moments to include the point x,y
void    moms_acc_pt( momentsT* M, double x, double y)
{   M->n += 1;
double  f = 1.0/M->n;
double  dx = x-M->xbar;
double  dy = y-M->ybar;
    M->xbar += f*dx;
    M->ybar += f*dy;
double  g = 1.0 - f;
    M->Vx   = g*(M->Vx  + f*dx*dx);
    M->Cxy  = g*(M->Cxy + f*dx*dy);
    M->Vy   = g*(M->Vy  + f*dy*dy);
}

// return the moments for the combination of A and B (assumed disjoint)
momentsT    moms_combine( const momentsT* A, const momentsT* B)
{
momentsT    C;
    C.n = A->n + B->n;
double  alpha = (double)A->n/(double)C.n;
double  beta = (double)B->n/(double)C.n;
    C.xbar = alpha*A->xbar + beta*B->xbar;
    C.ybar = alpha*A->ybar + beta*B->ybar;
double  dx = A->xbar - B->xbar;
double  dy = A->ybar - B->ybar;
    C.Vx = alpha*A->Vx + beta*B->Vx + alpha*beta*dx*dx;
    C.Cxy= alpha*A->Cxy+ beta*B->Cxy+ alpha*beta*dx*dy;
    C.Vy = alpha*A->Vy + beta*B->Vy + alpha*beta*dy*dy;
    return C;
}

// line is y^ : x -> a*x + b; return Sum{ sqr( y[i] - y^(x[i])) }/N
double  moms_line( momentsT* M, double* a, double *b)
{   *a = M->Cxy/M->Vx;
    *b = M->ybar - *a*M->xbar;
    return fmax( 0.0, M->Vy - *a*M->Cxy);
}