最小二乘法

#include <stdio.h>
#include <stdlib.h>
#include "math.h"
void polyfit(int n, double *x, double *y, int poly_n, double p[]);
void gauss_solve(int n,double A[],double x[],double b[]);
double xy(double x,double p[]);
/*==================polyfit(n,x,y,poly_n,a)===================*/
/*=======拟合y=a0+a1*x+a2*x^2+……+apoly_n*x^poly_n========*/
/*=====n是数据个数 xy是数据值 poly_n是多项式的项数======*/
/*===返回a0,a1,a2,……a[poly_n]，系数比项数多一（常数项）=====*/
/*
由x得出y

*/
void polyfit(int n,double x[],double y[],int poly_n,double p[])
{
int i,j;
double *tempx,*tempy,*sumxx,*sumxy,*ata;

tempx = (double *)calloc(n , sizeof(double));
sumxx = (double *)calloc((poly_n*2+1) , sizeof(double));
tempy = (double *)calloc(n , sizeof(double));
sumxy = (double *)calloc((poly_n+1) , sizeof(double));
ata = (double *)calloc( (poly_n+1)*(poly_n+1) , sizeof(double) );
for (i=0;i<n;i++)
{
tempx[i]=1;
tempy[i]=y[i];
}
for (i=0;i<2*poly_n+1;i++)
{
for (sumxx[i]=0,j=0;j<n;j++)
{
sumxx[i]+=tempx[j];
tempx[j]*=x[j];
}
}
for (i=0;i<poly_n+1;i++)
{
for (sumxy[i]=0,j=0;j<n;j++)
{
sumxy[i]+=tempy[j];
tempy[j]*=x[j];
}
}
for (i=0;i<poly_n+1;i++)
{
for (j=0;j<poly_n+1;j++)
{
ata[i*(poly_n+1)+j]=sumxx[i+j];
}
}
gauss_solve(poly_n+1,ata,p,sumxy);

free(tempx);
free(sumxx);
free(tempy);
free(sumxy);
free(ata);
}
/*============================================================
高斯消元法计算得到 n 次多项式的系数
n: 系数的个数
ata: 线性矩阵
sumxy: 线性方程组的Y值
p: 返回拟合的结果
============================================================*/
void gauss_solve(int n,double A[],double x[],double b[])
{
int i,j,k,r;
double max;
for (k=0;k<n-1;k++)
{
max=fabs(A[k*n+k]); // find maxmum
r=k;
for (i=k+1;i<n-1;i++)
{
if (max<fabs(A[i*n+i]))
{
max=fabs(A[i*n+i]);
r=i;
}
}
if (r!=k)
{
for (i=0;i<n;i++) //change array:A[k]&A[r]
{
max=A[k*n+i];
A[k*n+i]=A[r*n+i];
A[r*n+i]=max;
}
max=b[k]; //change array:b[k]&b[r]
b[k]=b[r];
b[r]=max;
}

for (i=k+1;i<n;i++)
{
for (j=k+1;j<n;j++)
A[i*n+j]-=A[i*n+k]*A[k*n+j]/A[k*n+k];
b[i]-=A[i*n+k]*b[k]/A[k*n+k];
}
}

for (i=n-1;i>=0;x[i]/=A[i*n+i],i--)
{
for (j=i+1,x[i]=b[i];j<n;j++)
x[i]-=A[i*n+j]*x[j];
}
}

/*************************************************
数据解析,可用于仪器校准，因为场地，环境等不同，有测量误差，进行自动校准，然后测量数据
例如，纸张计算中，比赛前，一张纸一张纸累积侧过去，然后得到校准参数，用于比赛中的测量
*************************************************/
double xy(double x,double p[])
{
double y;
y=p[0]+x*p[1]+x*x*p[2]+x*x*x*p[3]+x*x*x*x*p[4]+x*x*x*x*x*p[5];
return y;
}

void main()
{
int i, sizenum;
double P[6];
int dimension = 5; //5次多项式拟合,即5阶
// 要拟合的数据
double xx[]= {0.995119, 2.001185, 2.999068, 4.001035, 4.999859, 6.004461, 6.999335,
7.999433, 9.002257, 10.003888, 11.004076, 12.001602, 13.003390, 14.001623, 15.003034,
16.002561, 17.003010, 18.003897, 19.002563, 20.003530};
double yy[] = {-7.620000, -2.460000, 108.760000, 625.020000, 2170.500000, 5814.580000,
13191.840000, 26622.060000, 49230.220000, 85066.500000, 139226.280000, 217970.140000, 328843.860000,
480798.420000, 684310.000000, 951499.980000, 1296254.940000, 1734346.660000, 2283552.120000, 2963773.500000};

sizenum = sizeof(xx)/ sizeof(xx[0]); // 拟合数据的维数
polyfit(sizenum, xx, yy, dimension, P);

printf("拟合系数, 按升序排列如下:\n");
for (i=0;i<dimension+1;i++) //这里是升序排列，Matlab是降序排列
{
printf("P[%d]=%lf\n",i,P[i]);
}

printf("解得=%lf\n",xy(20.003530,P));
}

 

 还有种办法，查表，根据每个阶段对应实际值，从而运用算法得到对应值