最小二乘法(c语言实现线性,matlab进行拟合)及相关系数的求解

现在给定n个点，（x1,y1),(x2,y2),(x3,y3),(x4,y4),(x5,y5)..(xn,yn).

现在希望得到一条最好的曲线（也就是求一个函数关系式~~）

能尽可能的描述这n个点（不一定所有点都经过，但是总的拟合最小）

现在探讨什么叫总的拟合误差最小：

为了方便，我们考虑最简单的线性模型。

1:

n

∑（yi-f(xi)) 

i=1

但是考虑到有些点在线上面，有些点在线下面，正负会相互抵消，所以我们很容易想到另一种方法。

2:

n

∑（yi-f(xi)) *（yi- f(xi))

i=1

然后令f(x)=a+bx;

可以得到公式

n

∑((yi-(a+bxi))^2)

i=1
上式分别对a、b求偏导

整理后得到方程组

aN+b∑xi=∑yi

a∑xi+b∑(xi*xi)=∑xi*∑yi

解上述方程组（上面乘∑xi,下面式子乘N）便可求得直线参数a和b的最佳估计值。



说句题外话，以上的步骤只需要掌握求导的知识，手算时间不超过2分钟
当然，我们得到的拟合结果必然有好有坏，比如对于线性模型，直线的效果自然最好，直接得到y=kx+b的系数k,b，但是如果不是点不都在直线上，效果肯定不咋样。

所以引出最小二乘法的相关系数



当rxy接近1时说明拟合效果不错，而当接近0时，说明结果很差，基本无关。

至于用程序实现就很简单了，都是基础的运算。

#include <stdio.h>
#include <math.h>
#define N 10
double calculatesumx(double x[])
{
int i;
double sum=0;
for (i=0;i<N;i++)
sum=sum+x[i];
return sum;
}
double calculatesquare(double x[])
{
int i;
double sum=0;
for (i=0;i<N;i++)
sum=sum+x[i]*x[i];
return sum;
}
double calculatesumxy(double x[],double y[])
{
int i,j,n;
double sum=0;
for (i=0;i<N;i++)
{
sum=sum+x[i]*y[i];
}
return sum;
}
double sumaverage(double x[],double y[])
{
int i;
double averagex=0,averagey=0,sum=0;
for (i=0;i<N;i++)
{
averagex=averagex+x[i];
}
averagex=averagex/N;
for (i=0;i<N;i++)
{
averagey=averagey+x[i];
}
averagey=averagey/N;
for (i=0;i<N;i++)
{
sum=sum+(x[i]-averagex)*(y[i]-averagey);
}
return sum;
}
double squareaverage(double x[])
{
int i;
double sum=0,average=0;
for (i=0;i<N;i++)
average=average+x[i];
average=average/N;
for (i=0;i<N;i++)
{
sum=sum+(x[i]-average)*(x[i]-average);
}
return sqrt(sum);
}
int main()
{
double x[10] = {1,2,3,4,5,6,7,8,9,0};
double y[10] = {3,4,5,6,7,8,9,10,11,2};
int  i,j;
double sumx,sumy,sumxy,sumxx,averagesumxy,squareaveragex,squareaveragey,r,k,b;
sumx=calculatesumx(x);
sumy=calculatesumx(y);
sumxy=calculatesumxy(x,y);
sumxx=calculatesquare(x);
b=(sumxx*sumy-sumx*sumxy)/(N*sumxx-sumx*sumx);
k=(N*sumxy-sumx*sumy)/(N*sumxx-sumx*sumx);
printf("%.4f %.4f\n",k,b);//求得最小二乘法下的k,b，接下来评估拟合程度.
averagesumxy=sumaverage(x,y);
squareaveragex=squareaverage(x);
squareaveragey=squareaverage(y);
r=averagesumxy/(squareaveragex*squareaveragey);
printf("拟合程度为:%.4f\n",r);
return 0;

}

以上都是线性的拟合，然而实际上用到的更多的是更高次幂的拟合，所以就得用matlab啦。

在matlab中数据拟合主要有两种，一种是多项式的拟合，也就是告诉你一些点，求这些点拟合的曲线

x=[1 2 3 4 5 6 7 8 9];
y=[9 7 6 3 -1 2 5 7 20];
p=polyfit(x,y,3);%3是最大次幂，polyfit是多项式的个个系数
xi=0:2:10;%为了画图显示的范围，可根据数据范围自行调整。
yi=polyval(p,xi);%计算多项式的值
plot(xi,yi,x,y,'r*');%画出图形
%(x,y)是需要拟合的点的横纵坐标，3是多项式最大次数

另一种是指定函数的拟合，使用确实频繁。

它一般是用来求系数的，上一个是我不关心你用什么函数，只要在这个次幂下尽可能拟合就可以。

这个是告诉你具体的关系式，根据、点的坐标求出关系式的系数。

syms t
x=[1;2];%注意这里的数据都必须为列向量形式
y=[3;4];
f=fittype('a*t+b','independent','t','coefficients',{'a','b'});%这里t是告诉别人t是变量，而a,b都是待求系数
cfun=fit(x,y,f)
xi=0:1:20;
yi=cfun(xi);
plot(x,y,'r*',xi,yi,'b-');