C语言实现一元线性回归分析和多元线性回归分析

环能jvav大师

于 2024-07-15 08:15:00 发布

阅读量265

点赞数 6

文章标签： c语言线性回归回归机器学习算法

本文链接：https://blog.csdn.net/m0_73500130/article/details/140417382

版权

数学理论部分

线性归回

线性回归（Linear Regression）用于确定两种或两种以上变量间相互依赖的定量关系，是一种平时很常用的统计分析方法，它利用称为线性回归方程的最小平方函数对一个或多个自变量（特征值）和因变量（目标值）之间关系进行建模，这种函数是一个或多个称为回归系数的模型参数的线性组合。
线性回归又分为一元线性回归和多元线性回归：

一元线性回归分析：只包括一个自变量和一个因变量，且二者的关系可用一条直线近似表示。
在多元线性回归分析中，包括两个或两个以上的自变量，且因变量和自变量之间是线性关系。
多元线性回归分析：模型的形式通常为： $\beta_0 + \beta_1 x_1 + \beta_2 x_2 + \cdots + \beta_n x_n + \epsilon$
$y$ :因变量（响应变量），
$\beta_{0}, \beta_{1}, \ldots, \beta_{n}$ :回归系数，
$x_1, x_2, \ldots, x_n$ :自变量（预测变量），
$\epsilon$ :误差项。

最小二乘法

大部分的线性回归模型使用最小二乘法，最小二乘法基于最小化误差的平方和来求解模型参数，在自变量不是多重共线性的情况下，最小二乘法求解的模型参数具有唯一性，在给定数据集下，通过最小二乘法得到的模型是确定的。同时最小二乘法在正态分布假设下可以用`极大似然估计（MLE）``解释，也可以证明解是最优线性无偏估计，通过特征检测等小技巧，最小二乘法还可以用到处理非线性关系的数据的场景下。

在可视化方面看，最小二乘法试图找到一条直线，使得所有样本点到这条直线的欧氏距离之和最小，这种优化目标符合许多实际问题的需求。但是嘞，当数据中存在异常值的时候，最小二乘法会受到较大影响，它并非是万能的。

这次使用到的数据的最小二乘法一元线性回归趋势线如下：

在这里插入图片描述

准备数据

这次的数据和算法使用的数据都是从书籍《EXCEL机器学习》(作者：周红)中摘抄。
一元线性回归模型使用的温度与冰激凌销售额表格CSV文件如下(ice_cream.csv)：
数据一列为华氏温度(Temperature)，一列为冰激凌销售额(Ice_Cream_Sale)

Temperature,Ice_Cream_Sale
91,89.8
87,90.2
86,81.1
88,83
92.8,90.9
95.2,119
93.3,94.9
97.7,132.4

多元线性回归模型使用的温度、游客数和晴天数与冰激凌销售额表格CSV文件如下(ice_cream2.csv)：
数据由上一份表格修改，一列为华氏温度(Temperature)，一列为冰激凌销售额(Ice_Cream_Sale),一列为游客数(Number_of_SunnyDays)，一列为晴天数(Ice_Cream_Sale)

Temperature,Number_of_Tourists,Number_of_SunnyDays,Ice_Cream_Sale
91,998,4,89.8
87,1256,7,90.2
86,791,6,81.1
88,705,5,83
92.8,1089,3,90.9
95.2,1135,6,119
93.3,1076,4,94.9
97.7,1198,7,132.4

计算过程

一元线性回归方程

误差的平方和公式：
$E=\Sigma^{n}_{i=1}{(y_{i}-(mx_{i}+b))^{2}}$
这个公式用于计算实际值 $y_{i}$ 和预测值 $mx_{i}+b$ 之间的平方差的总和，其中 $n$ 是数据点的数量， $m$ 是斜率， $b$ 是截距。
将误差函数对 $m$ 和 $b$ 的偏导数设为零：
$\frac{\partial E}{\partial x} = 0 且 \frac{\partial E}{\partial y} = 0$
为了找到使误差最小的 $m$ 和 $b$ ，需要将误差函数对 $m$ 和 $b$ 的偏导数设为零。
斜率计算公式：
$\frac{n \Sigma^{n}_{i=1}x_{i}y_{i}-(\Sigma^{n}_{i=1}x_{i})(\Sigma^{n}_{i=1}y_{1})}{n(\Sigma^{n}_{i=1}x_{i}^{2})-(\Sigma^{n}_{i=1}x_{i})^{2}}$
使用了数据点的 $x_{i}$ 和 $y_{i}$ 值的总和以及它们的乘积的总和来计算斜率 $m$ 。
截距计算公式:
$b=\frac{(\Sigma^{n}_{i=1}{x_{i}^{2}})(\Sigma^{n}_{i=1}{y_{i}})-(\Sigma^{n}_{i=1}{x_{i}})(\Sigma^{n}_{i=1}{x_{i}y_{i}})}{n(\Sigma^{n}_{i=1}{x_{i}^{2}})(\Sigma^{n}_{i=1}{x_{i}})^{2}}$
使用数据点的 $x_{i}$ 和 $y_{i}$ 值的总和以及它们的乘积的总和来计算截距 $b$

多元线性回归方程

设计矩阵（Design Matrix）
设计矩阵 $X$ 用于将模型中的自变量组织起来，通常包括一列常数项（用于估计截距 $\beta_0$ ）：
$\begin{bmatrix} 1 & x_{11} & x_{12} & \cdots & x_{1n} \\ 1 & x_{21} & x_{22} & \cdots & x_{2n} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & x_{m1} & x_{m2} & \cdots & x_{mn} \end{bmatrix}$
$m$ 是观测值的数量，
$n$ 是自变量的数量（不包括常数项）。
因变量矩阵 $Y$
因变量矩阵 $Y$ 是一个列向量，包含了所有观测值的因变量值：
$\begin{bmatrix} y_1 \\ y_2 \\ \vdots \\ y_m \end{bmatrix}$
最小二乘法求解回归系数
回归系数 $\beta$ 可以通过最小二乘法求解，具体是通过正规方程（Normal Equations）：
$\beta = (X^T X)^{-1} X^T Y$
$X^T$ :设计矩阵 $X$ 的转置，
$X^T X)^{-1}$ 是 $X^T X$ 的逆矩阵。
高斯-约旦消元法（Gauss-Jordan Elimination）
1.在每一个循环过程中寻找到主元，并通过行变换将主元移动到矩阵的主对角线上。
2.将主元所在的行内的所有元素除以主元，使得主元化为1。
3.观察主元所在的列上的其他元素，通过行变换将它们所在的行减去主元所在的行乘以一定的倍数，使得主元所在的列内、除主元外的其他元素化为0。
4.循环重复上述步骤，直到处理完所有列。
5.最终，当系数矩阵被变换为单位矩阵时，增广矩阵的右侧部分会直接给出方程组的解。

C语言实现

一元线性回归分析

在这里使用了<stdio.h>头文件中预定义的结构体类型FILE,这个结构体用于处理文件输入输出，随后定义了一个指针fp指向了FILE，随后使用fopen标准函数打开csv文件。随后定义变量用于存储数据和累积和并定义缓冲区用于读取文件中的每一行。随后根据公式来计算斜率和截距。

#include <stdio.h>  
  
int main() {  
    FILE *fp = fopen("ice_cream.csv", "r");  
    if (fp == NULL) {  
        perror("没打开");  
        return -1;  
    }  
  
    double temp, sale, sumX = 0, sumY = 0, sumXY = 0, sumX2 = 0;  
    int n = 0;  
    char buffer[1024];  
  
    while (fgets(buffer, 1024, fp)) {  
        if (sscanf(buffer, "%lf,%lf", &temp, &sale) == 2) {  
            sumX += temp;  
            sumY += sale;  
            sumXY += temp * sale;  
            sumX2 += temp * temp;  
            n++;  
        }  
    }  
  
    fclose(fp);  
  
    double m = (n * sumXY - sumX * sumY) / (n * sumX2 - sumX * sumX);  
    double b = (sumX2 * sumY - sumX * sumXY) / (n * sumX2 - sumX * sumX);  
  
    printf("斜率m= %lf\n", m);  
    printf("y轴上的截距b= %lf\n", b);  
  
    return 0;  
}

运行结果如下：
在这里插入图片描述
根据运行结果和公式 $y = m x + b$ 可以得到预估模型:
$3.8394 * 预计的温度 - 253.17 = 预测的销售额$

多元线性回归分析

使用C语言进行多元线性回归存在一点点难度，大概步骤如下：

1.先使用#define函数定义两个全局变量MAX_DATA_POINTS和NUM_FEATURES，分别是最大数据点和元素数量
2.定义数据结构DataPoint用于存储数据点和Matrix用于存储矩阵
3.定义矩阵操作函数:
createMatrix：根据给定的行数和列数创建一个新的矩阵，并返回指向该矩阵的指针。
destroyMatrix：销毁一个矩阵，释放其占用的内存
printMatrix：打印矩阵的内容
matrixMultiply：计算两个矩阵的乘积，并返回结果矩阵
transposeMatrix：计算并返回矩阵的转置
4.定义数据读取函数与矩阵构建函数：
readData：从CSV文件中读取数据，并存储到DataPoint数组中，同时忽略文件的第一行。
buildDesignMatrix：根据读取的数据点构建一个设计矩阵，该矩阵包括常数项（y轴截距）和指定的特征（温度、游客数、晴天数）。
5.定义多元线性回归分析函数：
solveGaussJordan：使用高斯-约旦消元法求解线性方程组，先构建增广矩阵[A|B]，然后执行高斯-约旦消元，最后提取解向量。
calculateParameters：这是进行多元线性回归分析的核心函数。它首先计算设计矩阵X的转置XT，然后计算XT*X和XT*Y，最后调用solveGaussJordan函数求解回归参数。
6.定义主函数main:
在随后的主函数中声明DataPoint数组和变量来存储数据点和数据大小 -> 调用readData函数从CSV文件中读取数据 -> 调用buildDesignMatrix函数构建设计矩阵X,构建因变量矩阵Y，其中包含每个数据点的冰激凌销售额 -> 调用calculateParameters函数计算回归参数。最后即可打印回归参数（y轴截距、温度、游客数、晴天数的系数），随后清理分配的内存，释放矩阵占用的空间。

#include <stdio.h>
#include <stdlib.h>
#nclude <string.h>
#include <math.h>

#define MAX_DATA_POINTS 100
#define NUM_FEATURES 3 // 3分别是温度、游客数、晴天数

// 存储数据点
typedef struct {
    double temperature;
    double tourists;
    double sunny_days;
    double ice_cream_sale;
} DataPoint;

// 矩阵类型
typedef struct {
    double** elements;
    int rows;
    int cols;
} Matrix;

// 声明函数
Matrix* createMatrix(int rows, int cols);
vid destroyMatrix(Matrix* matrix);
void printMatrix(Matrix* matrix);
Matrix* matrixMultiply(Matrix* A, Matrix* B);
Matrix* transposeMatrix(Matrix* A);
void readData(const char* filename, DataPoint* data, int* dataSize);
Matrix* buildDesignMatrix(DataPoint* data, int dataSize, int numFeatures);
Matrix* calculateParameters(Matrix* X, Matrix* Y);

// 创建矩阵
Matrix* createMatrix(int rows, int cols) {
    Matrix* m = (Matrix*)malloc(sizeof(Matrix));
    m->rows = rows;
    m->cols = cols;
    m->elements = (double**)malloc(rows * sizeof(double*));
    for (int i = 0; i < rows; i++) {
        m->elements[i] = (double*)malloc(cols * sizeof(double));
    }
    return m;
}

// 销毁矩阵
void destroyMatrix(Matrix* matrix) {
    for (int i = 0; i < matrix->rows; i++) {
        free(matrix->elements[i]);
    }
    free(matrix->elements);
    free(matrix);
}

// 打印矩阵
void printMatrix(Matrix* matrix) {
    for (int i = 0; i < matrix->rows; i++) {
        for (int j = 0; j < matrix->cols; j++) {
            printf("%f ", matrix->elements[i][j]);
        }
        printf("\n");
    }
}

// 矩阵乘法
Matrix* matrixMultiply(Matrix* A, Matrix* B) {
    if (A->cols != B->rows) {
        printf("矩阵乘法中维度不匹配\n");
        return NULL;
    }
    
    Matrix* C = createMatrix(A->rows, B->cols);
    for (int i = 0; i < A->rows; i++) {
        for (int j = 0; j < B->cols; j++) {
            C->elements[i][j] = 0;
            for (int k = 0; k < A->cols; k++) {
                C->elements[i][j] += A->elements[i][k] * B->elements[k][j];
            }
        }
    }
    return C;
}

// 矩阵转置
Matrix* transposeMatrix(Matrix* A) {
    Matrix* AT = createMatrix(A->cols, A->rows);
    for (int i = 0; i < A->rows; i++) {
        for (int j = 0; j < A->cols; j++) {
            AT->elements[j][i] = A->elements[i][j];
        }
    }
    return AT;
}

// 从文件中读取数据
void readData(const char* filename, DataPoint* data, int* dataSize) {
    FILE* fp = fopen(filename, "r");
    if (!fp) {
        printf("无法打开文件 %s\n", filename);
        exit(1);
    }

    char line[256];
    fgets(line, sizeof(line), fp);
    while (fgets(line, sizeof(line), fp) != NULL && *dataSize < MAX_DATA_POINTS) {
        sscanf(line, "%lf,%lf,%lf,%lf",
               &(data[*dataSize].temperature),
               &(data[*dataSize].tourists),
               &(data[*dataSize].sunny_days),
               &(data[*dataSize].ice_cream_sale));
        (*dataSize)++;
    }
    fclose(fp);
}

// 构建设计矩阵
Matrix* buildDesignMatrix(DataPoint* data, int dataSize, int numFeatures) {
    Matrix* X = createMatrix(dataSize, numFeatures + 1);
    for (int i = 0; i < dataSize; i++) {
        X->elements[i][0] = 1; // 常数项
        X->elements[i][1] = data[i].temperature;
        X->elements[i][2] = data[i].tourists;
        X->elements[i][3] = data[i].sunny_days;
    }
    return X;
}

Matrix* solveGaussJordan(Matrix* A, Matrix* B) {  
    int n = A->rows;  
    int m = B->cols;  
    if (n != A->cols) {  
        printf("矩阵A必须是方阵\n");  
        return NULL;  
    }  
    if (n != B->rows) {  
        printf("矩阵A和B的行数必须相等\n");  
        return NULL;  
    }  
  
    // 创建增广矩阵 [A|B]  
    Matrix* AB = createMatrix(n, n + m);  
    for (int i = 0; i < n; i++) {  
        for (int j = 0; j < n; j++) {  
            AB->elements[i][j] = A->elements[i][j];  
        }  
        for (int j = 0; j < m; j++) {  
            AB->elements[i][n + j] = B->elements[i][j];  
        }  
    }  
  
    // 高斯-约旦消元法  
    for (int k = 0; k < n; k++) {  
        // 查找主元  
        int maxRow = k;  
        for (int i = k + 1; i < n; i++) {  
            if (fabs(AB->elements[i][k]) > fabs(AB->elements[maxRow][k])) {  
                maxRow = i;  
            }  
        }  
  
        // 如果主元为0则无法继续  
        if (fabs(AB->elements[maxRow][k]) < 1e-9) {  
            printf("警告: 矩阵接近奇异，可能无法求解\n");  
            destroyMatrix(AB);  
            return NULL;  
        }  
  
        // 交换行  
        if (maxRow != k) {  
            for (int j = 0; j < n + m; j++) {  
                double temp = AB->elements[k][j];  
                AB->elements[k][j] = AB->elements[maxRow][j];  
                AB->elements[maxRow][j] = temp;  
            }  
        }  
  
        // 归一化当前行  
        double pivot = AB->elements[k][k];  
        for (int j = k; j < n + m; j++) {  
            AB->elements[k][j] /= pivot;  
        }  
  
        // 消元  
        for (int i = 0; i < n; i++) {  
            if (i != k) {  
                double factor = AB->elements[i][k];  
                for (int j = k; j < n + m; j++) {  
                    AB->elements[i][j] -= factor * AB->elements[k][j];  
                }  
            }  
        }  
    }  
  
    // 提取解向量  
    Matrix* solution = createMatrix(n, m);  
    for (int i = 0; i < n; i++) {  
        for (int j = 0; j < m; j++) {  
            solution->elements[i][j] = AB->elements[i][n + j];  
        }  
    }  
  
    // 清理增广矩阵  
    destroyMatrix(AB);  
  
    return solution;  
}
// 使用高斯消元法
Matrix* calculateParameters(Matrix* X, Matrix* Y) {
    Matrix* XT = transposeMatrix(X);
    Matrix* XT_X = matrixMultiply(XT, X);
    Matrix* XT_Y = matrixMultiply(XT, Y);
    Matrix* parameters = solveGaussJordan(XT_X, XT_Y);
    destroyMatrix(XT_X);
    destroyMatrix(XT_Y);
    destroyMatrix(XT);
    return parameters;
}

int main() {  
    DataPoint data[MAX_DATA_POINTS];  
    int dataSize = 0;  
    readData("ice_cream2.csv", data, &dataSize);  
    // 构建设计矩阵 
    Matrix* X = buildDesignMatrix(data, dataSize, NUM_FEATURES);
    // 构建因变量矩阵 
    Matrix* Y = createMatrix(dataSize, 1);  
    for (int i = 0; i < dataSize; i++) {  
        Y->elements[i][0] = data[i].ice_cream_sale;  
    }  
  
    // 计算并打印回归参数  
    Matrix* parameters = calculateParameters(X, Y);  
    printf("回归参数如下:\n");  
    printf("y轴截距: %f\n", parameters->elements[0][0]); 
    printf("温度: %f\n", parameters->elements[1][0]); 
    printf("游客数: %f\n", parameters->elements[2][0]);
    printf("晴天数: %f\n", parameters->elements[3][0]); 

    // 清理分配的内存  
    destroyMatrix(X);  
    destroyMatrix(Y);  
    destroyMatrix(parameters);  
  
    return 0;  
}

运行结果如下：
在这里插入图片描述

环能jvav大师

关注

6
点赞
踩
7

收藏

觉得还不错? 一键收藏
0
评论
C语言实现一元线性回归分析和多元线性回归分析

使用最小二乘法进行一元或多元线性回归分析，优点是模型简单直观、计算效率高且广泛应用于多领域，能够解释自变量对因变量的影响。然而，它也存在对非线性关系拟合效果差、对异常值敏感、假设条件严格等缺点，尤其是在多元线性回归中还需注意多重共线性问题，且仅能表明相关性而非因果关系。
复制链接

扫一扫