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

前言

因为被数学作业逼疯了所以写了个这个小东西（强哥暗杀名单） 原理为最小二乘法。 (hat{b}=frac{sumlimits^n_{i=1}x_iy_i-noverline{x}overline{y}}{sumlimits^n_{i=1}x_i^2-noverline{x}^2}) (hat{a}=overline{y}-hat{b}overline{x})

使用方法

第一行输入数据组数。 第二行输入$x$的值，用空格隔开。 第三行输入对应的$y$的值，用空格隔开。 数据组数不能超过10000。 注：在实际计算中可能会因为$hat(保留的不同位数而导致)hat$的细小误差。默认保留六位小数。 UPD：然后我发现有的题的$hat$竟然能和给的数据差4...所以精度可能还是会有一些问题，理性使用吧。 UPD：增加对残差平方和$sumlimits_^nhat^{2$和相关指数$R}2$的计算。 (sumlimits_{i=1}^nhat{e}^2=sumlimits_{i=1}^n(y_i-hat{y}_i)^2) (R^2=1-frac{sumlimits_{i=1}^n(y_i-hat{y}_i)^2}{sumlimits^n_{i=1}(y_i-overline{y})^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
a=%lf
e^2=%lf
R^2=%lf
",b,a,sum1,R2);

        system("pause");
	return 0;
}

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

UPD：增加了对二次函数模型和指数函数模型转化的线性回归方程的计算： (hat{y}=c_1x^2+c_2) (hat{y}=c_3e^{c_4x})

二次函数模型

#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
a=%lf
e^2=%lf
R^2=%lf
",b,a,sum1,R2);

        system("pause");
	return 0;
}

指数函数模型

• 令$z=ln y$，本代码输出对应$z=hatx+hat(的)hat(和)hat(。则对应回归方程为)hat=e^{hatx+hat}$
• 注：本代码$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
a=%lf
e^2=%lf
R^2=%lf",b,a,sum1,R2);

        system("pause");
	return 0;
}

下载

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