最小二乘法求线性回归方程

前言

因为被数学作业逼疯了所以写了个这个小东西（强哥暗杀名单）
原理为最小二乘法。

使用方法

第一行输入数据组数。
第二行输入\(x\)的值，用空格隔开。
第三行输入对应的\(y\)的值，用空格隔开。
数据组数不能超过10000。
注：在实际计算中可能会因为\(\hat{b}\)保留的不同位数而导致\(\hat{a}\)的细小误差。默认保留六位小数。
UPD：然后我发现有的题的\(\hat{a}\)竟然能和给的数据差4...所以精度可能还是会有一些问题，理性使用吧。
UPD：增加对残差平方和\(\sum\limits_{i=1}^n\hat{e}^2\)和相关指数\(R^2\)的计算。

UPD：添加system("pause");，方便查看结果。

样例输入

5
15.0 25.8 30.0 36.6 44.4
39.4 42.9 42.9 43.1 49.2

样例输出

b=0.291046
a=34.663848
e^2=8.426560
R^2=0.832073

代码

#include <bits/stdc++.h>
using namespace std;
const int maxn=1e4+10;
int n;
double x[maxn],y[maxn];
double a,b,R2;

double cal(double k){
    return b*k+a;
}

int main(){
    scanf("%d",&n);

    double avex=0,avey=0;
    for(int i=1;i<=n;i++){
        scanf("%lf",&x[i]);
        avex+=x[i];
    }
    for(int i=1;i<=n;i++){
        scanf("%lf",&y[i]);
        avey+=y[i];
    }

    avex/=n;avey/=n;

    double sum1=0,sum2=0;
    for(int i=1;i<=n;i++){
        sum1+=x[i]*y[i];
        sum2+=x[i]*x[i];
    }

    b=(sum1-n*avex*avey)/(sum2-n*avex*avex);
    a=avey-b*avex;

    sum1=0,sum2=0;
    for(int i=1;i<=n;i++){
        sum1+=(y[i]-cal(x[i]))*(y[i]-cal(x[i]));
        sum2+=(y[i]-avey)*(y[i]-avey);
    }

    R2=1-sum1/sum2;

    printf("b=%lf\na=%lf\ne^2=%lf\nR^2=%lf\n",b,a,sum1,R2);

        system("pause");
    return 0;
}

函数模型转化的线性回归方程

UPD：增加了对二次函数模型和指数函数模型转化的线性回归方程的计算：

二次函数模型

#include <bits/stdc++.h>
using namespace std;
const int maxn=1e4+10;
int n;
double x[maxn],y[maxn];
double a,b,R2;

double cal(double k){
    return b*k+a;
}

int main(){
    scanf("%d",&n);

    double avex=0,avey=0;
    for(int i=1;i<=n;i++){
        scanf("%lf",&x[i]);
        x[i]*=x[i];//其实就多了这一句...
        avex+=x[i];
    }
    for(int i=1;i<=n;i++){
        scanf("%lf",&y[i]);
        avey+=y[i];
    }

    avex/=n;avey/=n;

    double sum1=0,sum2=0;
    for(int i=1;i<=n;i++){
        sum1+=x[i]*y[i];
        sum2+=x[i]*x[i];
    }

    b=(sum1-n*avex*avey)/(sum2-n*avex*avex);
    a=avey-b*avex;

    sum1=0,sum2=0;
    for(int i=1;i<=n;i++){
        sum1+=(y[i]-cal(x[i]))*(y[i]-cal(x[i]));
        sum2+=(y[i]-avey)*(y[i]-avey);
    }

    R2=1-sum1/sum2;

    printf("b=%lf\na=%lf\ne^2=%lf\nR^2=%lf\n",b,a,sum1,R2);

        system("pause");
    return 0;
}

指数函数模型

• 令\(z=ln\ y\)，本代码输出对应\(z=\hat{b}x+\hat{a}\)的\(\hat{a}\)和\(\hat{b}\)。则对应回归方程为\(\hat{y}=e^{\hat{b}x+\hat{a}}\)
• 注：本代码\(e\)的不同取值也可能导致微小误差，本代码取2.718281。日常计算通常取2.7，可能导致\([-0.01,0.01]\)的误差。

#include <bits/stdc++.h>
using namespace std;
const int maxn=1e4+10;
int n;
double x[maxn],y[maxn],z[maxn];
double a,b,R2;

double cal(double k){
    return pow(2.718281,b*k+a);
}

int main(){
    scanf("%d",&n);

    double avex=0,avez=0;
    for(int i=1;i<=n;i++){
        scanf("%lf",&x[i]);
        avex+=x[i];
    }
    for(int i=1;i<=n;i++){
        scanf("%lf",&y[i]);
        z[i]=log(y[i]);
        avez+=z[i];
    }

    avex/=n;avez/=n;

    double sum1=0,sum2=0;
    for(int i=1;i<=n;i++){
        sum1+=x[i]*z[i];
        sum2+=x[i]*x[i];
    }

    b=(sum1-n*avex*avez)/(sum2-n*avex*avex);
    a=avez-b*avex;

    sum1=0,sum2=0;
    for(int i=1;i<=n;i++){
        sum1+=(y[i]-cal(x[i]))*(y[i]-cal(x[i]));
        sum2+=(y[i]-avez)*(y[i]-avez);
    }

    R2=1-sum1/sum2;

    printf("b=%lf\na=%lf\ne^2=%lf\nR^2=%lf",b,a,sum1,R2);

        system("pause");
    return 0;
}

下载

蓝奏云地址
密码:203m
考虑到编译不了的...就帮大家编译一下。
仅限Windows系统。
（如果你觉得还可以就点个推荐吧orz）