拟合由两条直线组成的曲线的亚二次算法

Question

问题是找到一组点表示的实值2D曲线与由两条线组成的折线最佳匹配。

暴力解法是为曲线的每个点找到“左”和“右”线性拟合，并选择误差最小的一对。我可以在迭代曲线的点时逐步计算出这两个线性适合度，但我找不到一种逐步计算误差的方法。因此，这种方法会导致二次复杂度。

问题是是否存在可以提供次二次复杂度的算法？

第二个问题是是否有方便的C++库可用于此类算法？

---

编辑 对于使用单条直线进行拟合，有以下公式：

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

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

- Vahagn

---

@MadPhysicist 你的意思是在没有实际线路的情况下从线路计算有效值吗？ - Vahagn

有点像。所有的量m、b、rms等都可以用sum(x)、sum(x^2)、N来表示。 - Mad Physicist

如果你愿意，我可以为你解方程。你介意我使用Python标记吗？这应该很容易转换成C++。 - Mad Physicist

展示一下你所做的事情：“我可以在迭代曲线点时逐步计算出两个线性拟合。” - chux - Reinstate Monica

如果 y 取决于 x 并且最佳拟合将误差从 (x，y) 减小到 (x，curve(x))，则所使用的公式是有意义的。 如果 x 和 y 是独立的，则需要不同的误差计算方法。 - chux - Reinstate Monica

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 = Σxi/N
μy = Σyi/N
m = Σ(xi - μx)(yi - μy) / Σ(xi - μx)2
b = μy - m * μx

平方误差的总和由以下公式给出：

e = Σ(yi - m * xi - b)2

这些可以用简单的代数展开为以下公式：

m = (Σxiyi - ΣxiΣyi/N) / (Σxi2 - (Σxi)2/N)
b = Σyi/N - m * Σxi/N
e = Σyi2 + m2 * Σxi2 + N * b2 - 2 * m * Σxiyi - 2 * b * Σyi + 2 * m * b * Σxi

因此，您可以遍历所有可能性并记录最小的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) 算法。

- Mad Physicist

---

公式中的最后三项必须有因子2。您能否请修正答案？ - Vahagn

Answer 2

如果有帮助的话，这里是我用于此类事情的一些C代码。它对Mad Physicist所说的内容没有太大的补充。

首先，一个公式。如果你通过一些点拟合一条直线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,)是为了防止近乎完美拟合的情况下舍入误差导致结果变为负数（在数学上是非负数）。

虽然有可能编写一个从momentsT中删除点的函数，但我发现通过复制momentsT并将点累加到副本中、获取线并将副本保留在适合该点的一侧、原始momentsT保留在另一侧更容易处理。

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);
}

- dmuir