二次线性回归方程二次线性回归法

转载

墨染青衫 2024-08-05 12:11:27

文章标签 二次线性回归方程 python 机器学习最小二乘法线性回归 文章分类 机器学习人工智能

文章目录

要求
原理
代码

矩阵的建立
矩阵各种函数

求逆

求转置
矩阵相乘
输出矩阵

不足

完整代码

Cpp文件
python文件
traindata.txt

报告形式

要求

二次线性回归方程二次线性回归法_python

这里采用的是最简单的最小二乘的方法来找到这个多元问题的线性模型，然后根据此模型来进行测试分类。

原理

二次线性回归方程二次线性回归法_最小二乘法_02

代码

根据流程我给出代码片段，后面贴上完整代码的链接

矩阵的建立

因为考虑到代码的通用性，不能够限制输入矩阵的大小，因此这里采用vector的方式来作为建立矩阵类型

//矩阵类型
#define MATRIX vector<vector<float>>
MATRIX L

但是需要注意的是这种类型不能够像平常使用vector那样去用它，必须要经过初始化才能够使用其可变长的特性，在初始化的时候，只有第0行是可以使用的，其它的必须再次resize才能够使用

说明：

宏定义含义：

二次线性回归方程二次线性回归法_机器学习_03

比如：

二次线性回归方程二次线性回归法_机器学习_04

或者

增加一行：

二次线性回归方程二次线性回归法_线性回归_05

矩阵各种函数

求逆

//计算每一行每一列的每个元素所对应的余子式，组成A*
//按第一行展开计算|A|
float getA(MATRIX arcs, int n)
{
    if (n == 1)
    {
        return arcs[0][0];
    }
    float ans = 0;
    MATRIX temp;
    temp.resize(C(arcs));
    for (int i = 0; i < C(arcs); i++)
        temp[i].resize(C(arcs));
    int i, j, k;
    for (i = 0; i < n; i++)
    {
        for (j = 0; j < n - 1; j++)
        {
            for (k = 0; k < n - 1; k++)
            {
                temp[j][k] = arcs[j + 1][(k >= i) ? k + 1 : k];
            }
        }
        float t = getA(temp, n - 1);
        if (i % 2 == 0)
        {
            ans += arcs[0][i] * t;
        }
        else
        {
            ans -= arcs[0][i] * t;
        }
    }
    return ans;
}
void getAStart(MATRIX arcs, int n, MATRIX &ans)
{
    if (n == 1)
    {
        ans.resize(1);
        ans[0].resize(1);
        ans[0][0] = 1;
        return;
    }
    int i, j, k, t;
    MATRIX temp;
    temp.resize(C(arcs));
    for (int i = 0; i < C(arcs); i++)
        temp[i].resize(C(arcs));
    for (i = 0; i < n; i++)
    {
        for (j = 0; j < n; j++)
        {
            for (k = 0; k < n - 1; k++)
            {
                for (t = 0; t < n - 1; t++)
                {
                    temp[k][t] = arcs[k >= i ? k + 1 : k][t >= j ? t + 1 : t];
                }
            }
            ans[j][i] = getA(temp, n - 1);
            if ((i + j) % 2 == 1)
            {
                ans[j][i] = -ans[j][i];
            }
        }
    }
}
//得到给定矩阵src的逆矩阵保存到des中。
bool GetMatrixInverse(MATRIX src, int n, MATRIX &des)
{
    des.resize(C(src));
    for (int i = 0; i < C(src); i++)
        des[i].resize(C(src));
    float flag = getA(src, n);
    MATRIX t;
    t.resize(C(src));
    for (int i = 0; i < C(src); i++)
        t[i].resize(C(src));
    if (flag == 0)
    {
        return false;
    }
    else
    {
        getAStart(src, n, t);
        for (int i = 0; i < n; i++)
        {
            for (int j = 0; j < n; j++)
            {
                des[i][j] = t[i][j] / flag;
            }
        }
    }
    return true;
}

求转置

/**
 * @brief 获取矩阵的转置矩阵 
 * 
 * @param src 原矩阵
 * @param des 保存转置矩阵
 */
void GetXT(MATRIX src, MATRIX &des)
{
    des.resize(C(src));
    for (int i = 0; i < C(src); i++)
        des[i].resize(R(src));
    for (int i = 0; i < R(src); i++)
    {
        for (int j = 0; j < C(src); j++)
        {
            des[j][i] = src[i][j];
        }
    }
}

矩阵相乘

/**
 * @brief 获取两矩阵相乘
 * 
 * @param arg1 相乘的矩阵1
 * @param arg2 相乘的矩阵2
 * @param des 保存相乘结果
 */
void GetMuti(MATRIX arg1, MATRIX arg2, MATRIX &des)
{
    des.resize(R(arg1));
    for (int i = 0; i < R(arg1); i++)
        des[i].resize(C(arg2));
    for (int i = 0; i < R(des); i++)
    {
        for (int j = 0; j < C(des); j++)
        {
            float _data = 0;
            for (int k = 0; k < C(arg1); k++)
            {
                _data += arg1[i][k] * arg2[k][j];
            }
            des[i][j] = _data;
        }
    }
}

输出矩阵

void printMAT(MATRIX _data)
{
    for (int i = 0; i < R(_data); i++)
    {
        for (int j = 0; j < C(_data); j++)
        {
            printf("%.3f ", _data[i][j]);
        }
        cout << endl;
    }
}

不足

没有考虑到X^TX不是满秩或者正定矩阵

完整代码

Cpp文件

// #define TEST
#define RUN
/**
 * @file male.cpp
 * @author willpower
 * @brief 基于最小二乘法的多元线性回归
 * @version 0.1
 * @date 2021-04-06
 * 
 * @copyright Copyright (c) 2021
 * 
 */
#include <iostream>
#include <string>
#include <vector>
#include <cmath>
#include <time.h>
#include <Windows.h>
using namespace std;
//训练数据路径
const string filePath = "./traindata.txt";
//存放结果路径
const string savePath = "./savedata.txt";
//获取行
#define R(x) x.size()
//获取列
#define C(x) x[0].size()
//数据行和列
#define N 17
#define M 3
//5次
#define TESTCNT 5
//2折
#define FOLDCNT 2
//矩阵类型
#define MATRIX vector<vector<float>>
// #define VER(x) 1 / (1.0 + exp(-(x)))
MATRIX trainData;
MATRIX X;
MATRIX Y;
//最后的线性模型参数
float A1 = 0, A2 = 0, B = 0;
//计算每一行每一列的每个元素所对应的余子式，组成A*
//按第一行展开计算|A|
float getA(MATRIX arcs, int n)
{
    if (n == 1)
    {
        return arcs[0][0];
    }
    float ans = 0;
    MATRIX temp;
    temp.resize(C(arcs));
    for (int i = 0; i < C(arcs); i++)
        temp[i].resize(C(arcs));
    int i, j, k;
    for (i = 0; i < n; i++)
    {
        for (j = 0; j < n - 1; j++)
        {
            for (k = 0; k < n - 1; k++)
            {
                temp[j][k] = arcs[j + 1][(k >= i) ? k + 1 : k];
            }
        }
        float t = getA(temp, n - 1);
        if (i % 2 == 0)
        {
            ans += arcs[0][i] * t;
        }
        else
        {
            ans -= arcs[0][i] * t;
        }
    }
    return ans;
}
void getAStart(MATRIX arcs, int n, MATRIX &ans)
{
    if (n == 1)
    {
        ans.resize(1);
        ans[0].resize(1);
        ans[0][0] = 1;
        return;
    }
    int i, j, k, t;
    MATRIX temp;
    temp.resize(C(arcs));
    for (int i = 0; i < C(arcs); i++)
        temp[i].resize(C(arcs));
    for (i = 0; i < n; i++)
    {
        for (j = 0; j < n; j++)
        {
            for (k = 0; k < n - 1; k++)
            {
                for (t = 0; t < n - 1; t++)
                {
                    temp[k][t] = arcs[k >= i ? k + 1 : k][t >= j ? t + 1 : t];
                }
            }
            ans[j][i] = getA(temp, n - 1);
            if ((i + j) % 2 == 1)
            {
                ans[j][i] = -ans[j][i];
            }
        }
    }
}
//得到给定矩阵src的逆矩阵保存到des中。
bool GetMatrixInverse(MATRIX src, int n, MATRIX &des)
{
    des.resize(C(src));
    for (int i = 0; i < C(src); i++)
        des[i].resize(C(src));
    float flag = getA(src, n);
    MATRIX t;
    t.resize(C(src));
    for (int i = 0; i < C(src); i++)
        t[i].resize(C(src));
    if (flag == 0)
    {
        return false;
    }
    else
    {
        getAStart(src, n, t);
        for (int i = 0; i < n; i++)
        {
            for (int j = 0; j < n; j++)
            {
                des[i][j] = t[i][j] / flag;
            }
        }
    }
    return true;
}
/**
 * @brief 将任意大小矩阵输出
 * 
 * @param _data 输入的矩阵
 */
void printMAT(MATRIX _data)
{
    for (int i = 0; i < R(_data); i++)
    {
        for (int j = 0; j < C(_data); j++)
        {
            printf("%.3f ", _data[i][j]);
        }
        cout << endl;
    }
}
/**
 * @brief 读取数据
 */
void ReadData()
{
    for (int i = 0; i < N; i++)
    {
        //每多一行，给矩阵加一行
        trainData.resize(i + 1);
        for (int j = 0; j < M; j++)
        {
            //每多一列，给矩阵加一列
            trainData[i].resize(j + 1);
            cin >> trainData[i][j];
        }
    }
#ifdef TEST
    printMAT(trainData);
#endif
}
//
/**
 * @brief 获取矩阵的转置矩阵 
 * 
 * @param src 原矩阵
 * @param des 保存转置矩阵
 */
void GetXT(MATRIX src, MATRIX &des)
{
    des.resize(C(src));
    for (int i = 0; i < C(src); i++)
        des[i].resize(R(src));
    for (int i = 0; i < R(src); i++)
    {
        for (int j = 0; j < C(src); j++)
        {
            des[j][i] = src[i][j];
        }
    }
}
/**
 * @brief 获取两矩阵相乘
 * 
 * @param arg1 相乘的矩阵1
 * @param arg2 相乘的矩阵2
 * @param des 保存相乘结果
 */
void GetMuti(MATRIX arg1, MATRIX arg2, MATRIX &des)
{
    des.resize(R(arg1));
    for (int i = 0; i < R(arg1); i++)
        des[i].resize(C(arg2));
    for (int i = 0; i < R(des); i++)
    {
        for (int j = 0; j < C(des); j++)
        {
            float _data = 0;
            for (int k = 0; k < C(arg1); k++)
            {
                _data += arg1[i][k] * arg2[k][j];
            }
            des[i][j] = _data;
        }
    }
}
/**
 * @brief 计算最小二乘的参数 
 * 
 * @param arg1 XTX的逆
 * @param arg2 XT
 * @param arg3 Y
 * @param _a1 线性模型参数
 * @param _a2 
 * @param _b 
 */
void CalcArg(MATRIX arg1, MATRIX arg2, MATRIX arg3, float &_a1, float &_a2, float &_b)
{
    MATRIX res1;
    MATRIX res2;
    GetMuti(arg1, arg2, res1);
    GetMuti(res1, arg3, res2);
    _a1 = res2[0][0];
    _a2 = res2[1][0];
    _b = res2[2][0];
}
/**
 * @brief 输出参数
 * 
 * @param a1 线性参数1
 * @param a2 线性参数2
 * @param b 线性参数3
 */
void PrintArg(float &a1, float &a2, float &b)
{
// #ifdef RUN
    cout << "a1:" << a1 << " a2:" << a2 << " b:" << b << endl;
// #endif
}
/**
 * @brief 训练数据，并且进行测试
 * 
 * @param _trainData 训练的数据
 * @param _testData 测试数据
 * @param score 得分
 * @param enTrain 是否为训练
 */
void testRun(MATRIX _trainData, MATRIX _testData, float score, bool enTrain)
{
    float a1, a2, b;
    if (enTrain)
    {
        MATRIX XTX;
        MATRIX XTX1;
        MATRIX XT;
        GetXT(_trainData, XT);
        GetMuti(XT, _trainData, XTX);
        GetMatrixInverse(XTX, C(XTX), XTX1);
        CalcArg(XTX1, XT, Y, a1, a2, b);
        A1 += a1;
        A2 += a2;
        B += b;
        PrintArg(a1, a2, b);
    }
    else
    {
        a1 = A1 / (TESTCNT * FOLDCNT);
        a2 = A2 / (TESTCNT * FOLDCNT);
        b = B / (TESTCNT * FOLDCNT);
        PrintArg(A1, A2, B);
    }
    float scores = 0;
    for (int i = 0; i < R(_testData); i++)
    {
        if (((a1 * _testData[i][0] + a2 * _testData[i][1] + b) > 0.5) && (_testData[i][2] == 1.0))
            scores += 1;
        else if (((a1 * _testData[i][0] + a2 * _testData[i][1] + b) < 0.5) && (_testData[i][2] == 0))
        {
            scores += 1;
        }
#ifdef RUN
        cout << ((a1 * _testData[i][0]) + (a2 * _testData[i][1]) + b) << endl;
#endif
    }
    score = scores / R(_testData);
// #ifdef RUN
    cout << "score:" << score << endl;
    cout << "************************************" << endl;
// #endif
    //延时1s保证随机
    Sleep(1000);
}
/**
 * @brief 将训练数据矩阵做出处理
 * 
 * @param _src 输入的训练数据
 * @param _des 处理后数据存放的位置
 */
void GenerateMAT(MATRIX _src, MATRIX &_des)
{
    _des.resize(R(_src));
    for (int i = 0; i < R(_src); i++)
        _des[i].resize(C(_src));
    for (int i = 0; i < R(_src); i++)
    {
        for (int j = 0; j < C(_src) - 1; j++)
        {
            _des[i][j] = _src[i][j];
        }
    }
    for (int i = 0; i < R(_des); i++)
    {
        _des[i][C(_des) - 1] = 1;
    }
#ifdef TEST
    cout << "src marix" << endl;
    printMAT(_src);
    cout << "des marix" << endl;
    printMAT(_des);
#endif
}
int main()
{
    //重定向输入为训练数据文件
    freopen(filePath.c_str(), "rb", stdin);
    //将数据全部取出
    ReadData();
    float score;
    //5次2折验证
    for (int cnt = 0; cnt < TESTCNT; cnt++)
    {
        //以时间作为种子，保证随机将数据对折
        srand((unsigned)time(NULL));
        //其中X为训练数据，testX为测试数据
        MATRIX X, testX, trainX;
        //将数据随机划分
        for (int i = 0; i < R(trainData); i++)
        {
            if ((rand() % 2) && (R(testX) < (N / 2)))
            {
                //增加行数
                testX.resize(R(testX) + 1);
                //因为列数一定，因此采用M初始化
                testX[R(testX) - 1].resize(M);
                for (int j = 0; j < M; j++)
                {
                    testX[R(testX) - 1][j] = trainData[i][j];
                }
                continue;
            }
            if (R(trainX) < (N - (N / 2)))
            {
                trainX.resize(R(trainX) + 1);
                Y.resize(R(Y) + 1);
                for (int j = 0; j < M; j++)
                {
                    trainX[R(trainX) - 1].resize(j + 1);
                    trainX[R(trainX) - 1][j] = trainData[i][j];
                    if ((M - 1) == j)
                    {
                        Y[R(Y) - 1].resize(1);
                        Y[R(Y) - 1][0] = trainData[i][j];
                    }
                }
                continue;
            }
        }
        #ifdef TEST
        cout<<"train data "<<endl;
        printMAT(trainX);
        cout<<"test data "<<endl;
        printMAT(testX);
        cout<<"********************"<<endl;
        #endif
        GenerateMAT(trainX, X);
#ifdef RUN
        printMAT(testX);
#endif
        testRun(X, testX, score, 1);
        /**
         @brief 此处将二折的训练集和测试集互换
         * 
         */
        GenerateMAT(testX, X);
#ifdef RUN
        printMAT(trainX);
#endif
        //训练数据并且进行测试
        testRun(X, trainX, score, 1);
    }
#ifdef RUN
    printMAT(trainData);
#endif
    //将训练后参数的平均值作为最后的参数，然后将整个训练数据作为测试数据进行测试
    testRun(X, trainData, score, 0);
    freopen(savePath.c_str(), "wb", stdout);
    cout<<A1<<" "<<A2<<" "<<B<<endl;
    return 0;
}

python文件

# encoding:utf-8
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from mpl_toolkits import mplot3d
import sys
PATH='./traindata.txt'
ARG = "./savedata.txt"
data = []
with open(PATH) as f:
    for line in f:
        data.append(line[:-1].split(' '))
f.close
arg = []
with open(ARG) as f:
    for line in f:
        arg.append(line[:-1].split(' '))
f.close
for i in range(len(arg[0])):
    arg[0][i] = float(arg[0][i])
print(data)
x1,y1,x2,y2 = [],[],[],[]
for da in data:
    if(da[2] == '1'):
        x1.append(float(da[0]))
        y1.append(float(da[1]))
    else:
        x2.append(float(da[0]))
        y2.append(float(da[1]))
x = np.arange(0, 0.8, 0.01)
y = (0.5-(arg[0][0]*x+arg[0][2]))/arg[0][1]
plt.rcParams['font.sans-serif'] = ['SimHei']  # 用来正常显示中文标签
plt.scatter(x1,y1,s=100,c='g',label='好瓜')
plt.scatter(x2,y2,s=100,c='r',label='坏瓜')
plt.xlabel("密度")
plt.ylabel("含糖量")
plt.legend()#显示点注释
plt.plot(x,y)
print(arg)
plt.plot()
plt.show()

traindata.txt

0.697 0.460 1
0.774 0.376 1
0.634 0.264 1
0.608 0.318 1
0.556 0.215 1
0.403 0.237 1
0.481 0.149 1
0.437 0.211 1
0.666 0.091 0
0.243 0.267 0
0.245 0.057 0
0.343 0.099 0
0.639 0.161 0
0.657 0.198 0
0.360 0.370 0
0.593 0.042 0
0.719 0.103 0