最小二乘法去趋势项

Paul George 13

已于 2024-03-22 14:19:27 修改

阅读量354

点赞数 5

文章标签：最小二乘法算法 stm32 c语言

于 2024-03-22 10:32:19 首次发布

本文链接：https://blog.csdn.net/weixin_49876887/article/details/136933265

版权

模型为 y = ax + b

这些公式是通过最小化误差的平方和来推导出来的。对于线性模型 $y = ax + b$ ，我们希望找到使得拟合的直线最小化实际数据点到该直线的距离的参数 a 和 b 。

考虑拟合直线 $y = ax + b$ 对于给定数据点 $(x_i, y_i)$ 的预测值 $\widehat{y_i}$ 与实际值 $y_i$ 之间的残差 $e_i = y_i - \hat{y_i}$ 。

拟合直线的预测值为： $\hat{y_i} = ax_i + b$

每个数据点 $(x_i, y_i)$ 到拟合直线的残差为： $e_i = y_i - (ax_i + b)$

我们的目标是找到使得残差的平方和最小的参数 a 和 b ，即最小化误差的平方和 S ：

$S = \sum_{i=1}^{n} e_i^2 = \sum_{i=1}^{n} (y_i - (ax_i + b))^2$

为了最小化 S ，我们可以对 a 和 b 分别求偏导数，并令其等于零，然后解出参数的值。

对 a 求偏导数：

$\frac{\partial S}{\partial a} = -2 \sum_{i=1}^{n} x_i(y_i - (ax_i + b))$

令上式等于零，并整理可得：

$\sum_{i=1}^{n} (x_i y_i - x_i(ax_i + b)) = 0$

$\sum_{i=1}^{n} (x_i y_i - a x_i^2 - b x_i) = 0$

$a \sum_{i=1}^{n} x_i^2 + b \sum_{i=1}^{n} x_i - \sum_{i=1}^{n} x_i y_i = 0$

$a \sum_{i=1}^{n} x_i^2 + b \sum_{i=1}^{n} x_i = \sum_{i=1}^{n} x_i y_i$

对 b 求偏导数：

$\frac{\partial S}{\partial b} = -2 \sum_{i=1}^{n} (y_i - (ax_i + b))$

令上式等于零，并整理可得：

$\sum_{i=1}^{n} (y_i - (ax_i + b)) = 0$

$\sum_{i=1}^{n} (y_i - ax_i - b) = 0$

$\sum_{i=1}^{n} y_i - a \sum_{i=1}^{n} x_i - nb = 0$

$a \sum_{i=1}^{n} x_i + nb = \sum_{i=1}^{n} y_i$

通过求解以上两个方程组，可以得到参数 $ a $ 和 $ b $ 的解析表达式。这些解就是我们需要的最小二乘估计值。

// 计算线性拟合的参数 a 和 b
void calculateTrend(int n, float x[], float y[], float *a, float *b) {
    float sum_x = 0.0, sum_y = 0.0, sum_xy = 0.0, sum_x2 = 0.0;

    // 计算各项和
    for (int i = 0; i < n; ++i) {
        sum_x += x[i];
        sum_y += y[i];
        sum_xy += x[i] * y[i];
        sum_x2 += x[i] * x[i];
    }

    // 计算参数 a 和 b
    *a = (n * sum_xy - sum_x * sum_y) / (n * sum_x2 - sum_x * sum_x);
    *b = (sum_x2 * sum_y - sum_x * sum_xy) / (n * sum_x2 - sum_x * sum_x);
}

// 去除趋势项
void removeTrend(int n, float x[], float y[], float a, float b, float detrended_y[]) {
    for (int i = 0; i < n; ++i) {
        detrended_y[i] = y[i] - (a * x[i] + b);
    }
}

模型为 y = a sin(b * pi * x)

要使用模型 $y = a \sin(b \pi x)$ 去除振动加速度数据的趋势项，可以通过最小二乘法来计算最佳的参数 $a$ 和 $b$ 。以下是这个过程的数学计算公式：

定义误差函数： 定义拟合误差为实际加速度数据点 $(x_i, y_i)$ 与模型预测值之间的残差的平方和，即： $S = \sum_{i=1}^{n} (y_i - \hat{y_i})^2 = \sum_{i=1}^{n} \left(y_i - a \sin(b \pi x_i)\right)^2$
最小化误差： 使用最小二乘法，通过最小化误差函数 $S$ 来找到使得误差最小的参数 $a$ 和 $b$ 。
计算参数 a 和 b ： 通过对误差函数求偏导数，并令其等于零，可以得到参数 a 和 b 的解析表达式。

首先，我们需要计算误差函数 S 对参数 a 和 b 的偏导数，并令其等于零，然后解出 a 和 b 的表达式。这将产生一个非线性方程组，需要使用数值方法来求解。

具体而言，误差函数对参数 a 和 b 的偏导数分别为：

对 a 求偏导数：

$\frac{\partial S}{\partial a} = -2 \sum_{i=1}^{n} \left(y_i - a \sin(b \pi x_i)\right) \sin(b \pi x_i)$

对 b 求偏导数：

$\frac{\partial S}{\partial b} = -2 \pi \sum_{i=1}^{n} \left(y_i - a \sin(b \pi x_i)\right) a x_i \cos(b \pi x_i)$

然后，我们令上述偏导数等于零，并解出 a 和 b 的值，以最小化误差函数 S 。

// 计算拟合误差
double calculateError(int n, double x[], double y[], double a, double b) {
    double error = 0.0;
    for (int i = 0; i < n; ++i) {
        double predicted_y = a * sin(b * M_PI * x[i]);
        error += pow(y[i] - predicted_y, 2);
    }
    return error;
}

// 使用最小二乘法拟合数据，计算参数 a 和 b
void fitData(int n, double x[], double y[], double *a, double *b) {
    double best_error = INFINITY;
    double best_a, best_b;

    // 尝试不同的 a 和 b 的值，并选择使得拟合误差最小的参数
    for (double a_val = 0.1; a_val <= 10.0; a_val += 0.1) {
        for (double b_val = 0.1; b_val <= 10.0; b_val += 0.1) {
            double error = calculateError(n, x, y, a_val, b_val);
            if (error < best_error) {
                best_error = error;
                best_a = a_val;
                best_b = b_val;
            }
        }
    }

    *a = best_a;
    *b = best_b;
}

// 去除趋势项
void removeTrend(int n, double x[], double y[], double a, double b, double detrended_y[]) {
    for (int i = 0; i < n; ++i) {
        detrended_y[i] = y[i] - (a * sin(b * M_PI * x[i]));
    }
}