资讯

历史

科技

环境与自然

成长

游戏

财经

文学与艺术

美食

健康

家居

文化

情感

汽车

三农

军事

旅行

运动

教育

生活

星座命理

如何用C语言实现最小二乘法

创作时间:

作者:

@小白创作中心

如何用C语言实现最小二乘法

引用

来源

https://docs.pingcode.com/baike/1285634

如何用C语言实现最小二乘法

核心观点：最小二乘法是一种用于求解数据拟合问题的数学方法、C语言可以通过矩阵运算和线性代数实现最小二乘法、通过编写C代码实现最小二乘法需要理解线性代数的基本概念。下面将详细解释如何用C语言实现最小二乘法。

最小二乘法是一种用于求解数据拟合问题的数学方法，它通过最小化数据点与拟合线之间的平方误差来找到最佳拟合线。要在C语言中实现最小二乘法，首先需要理解线性代数的基本概念，例如矩阵和向量运算。然后，可以通过编写C代码来实现这些运算，最终实现最小二乘法。

一、最小二乘法的原理

最小二乘法是通过最小化误差平方和来找到拟合数据的最佳曲线或直线。假设有一组数据点 $(x_i, y_i)$，拟合直线的方程为 $y = ax + b$。目标是找到系数 $a$ 和 $b$，使得所有数据点到拟合直线的距离平方和最小。

1、误差平方和公式

误差平方和（SSE）的公式为：

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

通过对 $a$ 和 $b$ 求导并设导数为零，可以得到一组线性方程：

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

2、求解方程组

通过求解上述方程组，可以得到 $a$ 和 $b$ 的值。使用矩阵表示，可以将其写成：

$$
A \cdot \beta = Y
$$

其中，

$$
A = \begin{pmatrix}
\sum_{i=1}^{n} x_i & n \
\sum_{i=1}^{n} x_i^2 & \sum_{i=1}^{n} x_i
\end{pmatrix}, \quad
\beta = \begin{pmatrix}
a \
b
\end{pmatrix}, \quad
Y = \begin{pmatrix}
\sum_{i=1}^{n} y_i \
\sum_{i=1}^{n} x_i y_i
\end{pmatrix}
$$

通过矩阵运算，可以求解出 $\beta$，即 $a$ 和 $b$。

二、C语言实现最小二乘法

1、引入必要的库和定义数据结构

首先，需要引入必要的库文件，并定义用于存储数据点和结果的结构体。

#include <stdio.h>
#include <stdlib.h>

typedef struct {
    double x;
    double y;
} DataPoint;

typedef struct {
    double a;
    double b;
} LinearFitResult;

2、实现矩阵运算函数

实现基本的矩阵运算函数，例如矩阵乘法和矩阵求逆。

void matrixMultiply(double mat1[2][2], double mat2[2][1], double res[2][1]) {
    for (int i = 0; i < 2; i++) {
        res[i][0] = 0;
        for (int j = 0; j < 2; j++) {
            res[i][0] += mat1[i][j] * mat2[j][0];
        }
    }
}

int matrixInverse(double mat[2][2], double inv[2][2]) {
    double det = mat[0][0] * mat[1][1] - mat[0][1] * mat[1][0];
    if (det == 0) {
        return -1;  // Matrix is not invertible
    }
    double invDet = 1.0 / det;
    inv[0][0] = mat[1][1] * invDet;
    inv[0][1] = -mat[0][1] * invDet;
    inv[1][0] = -mat[1][0] * invDet;
    inv[1][1] = mat[0][0] * invDet;
    return 0;
}

3、实现最小二乘法函数

实现用于执行最小二乘法的函数，计算出系数 $a$ 和 $b$。

LinearFitResult leastSquaresFit(DataPoint* data, int n) {
    LinearFitResult result;
    double sumX = 0, sumY = 0, sumXY = 0, sumX2 = 0;
    for (int i = 0; i < n; i++) {
        sumX += data[i].x;
        sumY += data[i].y;
        sumXY += data[i].x * data[i].y;
        sumX2 += data[i].x * data[i].x;
    }
    double A[2][2] = {{sumX2, sumX}, {sumX, n}};
    double Y[2][1] = {{sumXY}, {sumY}};
    double A_inv[2][2];
    if (matrixInverse(A, A_inv) != 0) {
        fprintf(stderr, "Matrix inversion failed.\n");
        exit(EXIT_FAILURE);
    }
    double B[2][1];
    matrixMultiply(A_inv, Y, B);
    result.a = B[0][0];
    result.b = B[1][0];
    return result;
}

4、测试最小二乘法函数

编写一个简单的主函数来测试最小二乘法函数。

int main() {
    DataPoint data[] = {
        {1, 1},
        {2, 2},
        {3, 2},
        {4, 4},
        {5, 5}
    };
    int n = sizeof(data) / sizeof(data[0]);
    LinearFitResult result = leastSquaresFit(data, n);
    printf("a: %f, b: %f\n", result.a, result.b);
    return 0;
}

三、优化和扩展

1、处理更大规模的数据集

对于更大规模的数据集，可以使用动态数组和内存管理来存储数据点，并优化矩阵运算的效率。

DataPoint* readDataFromFile(const char* filename, int* n) {
    FILE* file = fopen(filename, "r");
    if (!file) {
        perror("Failed to open file");
        exit(EXIT_FAILURE);
    }
    fscanf(file, "%d", n);
    DataPoint* data = (DataPoint*)malloc((*n) * sizeof(DataPoint));
    if (!data) {
        perror("Failed to allocate memory");
        exit(EXIT_FAILURE);
    }
    for (int i = 0; i < *n; i++) {
        fscanf(file, "%lf %lf", &data[i].x, &data[i].y);
    }
    fclose(file);
    return data;
}

2、支持多元线性回归

扩展最小二乘法以支持多元线性回归，处理多个自变量的情况。

typedef struct {
    double* x;
    double y;
} MultiVarDataPoint;

typedef struct {
    double* coefficients;
    int size;
} MultiVarLinearFitResult;

MultiVarLinearFitResult multiVarLeastSquaresFit(MultiVarDataPoint* data, int n, int m) {
    MultiVarLinearFitResult result;
    result.size = m + 1;
    result.coefficients = (double*)malloc(result.size * sizeof(double));
    if (!result.coefficients) {
        perror("Failed to allocate memory");
        exit(EXIT_FAILURE);
    }
    // Implement the multi-variable least squares fitting algorithm here
    // ...
    return result;
}