最小二乘法（高斯）

来源：互联网发布：CPDA数据分析师证书编辑：程序博客网时间：2024/05/22 10:34
# include<cstdio># include<cmath># define N 12   ///12个值# define M 4# include<cstring># include<iostream>using namespace std;///这是数据初始化。double x[N]= {0,5,10,15,20,25,30,35,40,45,50,55};double y[N]= {0,1.27,2.16,2.86,3.44,3.87,4.15,4.37,4.51,4.58,4.02,4.64};double xx[7];double yy[4];///这是要用数组，高斯消元要用double a[5][5]={    {0,0,0,0,0},    {0,0,0,0,0},    {0,0,0,0,0},    {0,0,0,0,0},    {0,0,0,0,0},};double b[5]={0,0,0,0,0};double nei(double t){    return b[1]*t+b[2]*t*t+b[3]*t*t*t;}int main(){    int n=3,k,i,j;    double sum,s,num;    ///这里是求x[i]、x[i]平方、、、、x[i]六次方的累和。就是正规表达式中的各项要求的x[i]部分。    for(i=1; i<=6; i++)    {        sum=0;        for(j=0; j<N; j++)        {            s=1;            for(k=1; k<=i; k++)                s=s*x[j];            sum=sum+s;        }        xx[i]=sum;    }    ///这里是看结果对不对。    for(i=1; i<=6; i++)    {        printf("xx[%d]=%lf\n",i,xx[i]);    }    ///这里是求正规表达式中的y[i]、、、、x[i]*x[i]*x[i]*y[i]的累和。就是正规表达式中的有y[i]的部分。    for(i=1; i<=4; i++)    {        sum=0;        for(j=0; j<12; j++)        {            s=1;            for(k=1; k<i; k++)                s=s*x[j];            sum=sum+s*y[j];        }        yy[i]=sum;    }    ///这里是看结果对不对。    for(i=1; i<=4; i++)    {        printf("yy[%d]=%lf\n",i,yy[i]);    }    ///这里是把上面求得的x[i],、、、、等放到数组里。    for(i=1; i<=3; i++)    {        b[i]=yy[i];        for(j=1; j<=3; j++)        {            a[i][j]=xx[i+j-1];        }    }     ///这里是看结果在数组放的位置对不对。    /*    for(i=1;i<=4;i++)    {        for(j=1;j<=4;j++)            printf("%lf  ", a[i][j]);        printf("\n");    }    */    /*    for(i=1;i<=4;i++)        printf("%lf  ",b[i]);    printf("\n");    for(i=1;i<=4;i++)    {        for(j=1;j<=4;j++)            printf("%lf  ", a[i][j]);        printf("\n");    }    */    ///后面是高斯消元的部分。    k=1;    do    {        for(i=k+1; i<=n; i++)            a[i][k]=a[i][k]/a[k][k];        for(i=k+1; i<=n; i++)        {            for(j=k+1; j<=n; j++)                a[i][j]=a[i][j]-a[i][k]*a[k][j];            b[i]=b[i]-a[i][k]*b[k];        }        if(k==n-1)            break;        k++;    }    while(1);    printf("\n");    for(i=1;i<=4;i++)    {        for(j=1;j<=4;j++)            printf("%lf  ", a[i][j]);        printf("\n");    }    b[n]=b[n]/a[n][n];    for(i=n-1; i>=1; i--)    {        s=0;        for(j=i+1; j<=n; j++)            s+=a[i][j]*b[j];        b[i]=(b[i]-s)/a[i][i];    }    ///这里的b[1],b[2],b[3]表示的是拟合函数中的a1,a2,a3.    for(i=1; i<=n; i++)        printf("b[%2d]=%lf\n",i,b[i]);    ///这里输出的是12个值的每一个误差。    for(i=0;i<12;i++)    {        printf("y[%d]=%lf   ",i,y[i]);        printf("%lf  ",nei(x[i]));        printf("error  %lf  ",y[i]-nei(x[i]));        printf("\n");    }    printf("\n");    ///这里求得是12个点平方和的总误差    num=0;    for(int i=0;i<12;i++)    {        num=(y[i]-nei(x[i]))*(y[i]-nei(x[i]));    }    printf("%lf ",num);    return 0;}
阅读全文
0 0