最小二乘算法 C 语言实现

//
//  main.m
//  test
//
//  Created by Jack on 15/11/20.
//  Copyright © 2015年 宇之楓鷙. All rights reserved.
//

#import <stdio.h>
#import <stdlib.h>
#import <malloc/malloc.h>
#import <math.h>


Smooth(double *x,double *y,double *a,int n,int m,double *dt1,double *dt2,double *dt3)//(x,y,a,n,m,dt1,dt2,dt3 )
//double *x; /*实型一维数组，输入参数，存放节点的xi值*/
//double *y; /*实型一维数组，输入参数，存放节点的yi值*/
//double *a; /*双精度实型一维数组，长度为m。返回m一1次拟合多项式的m个系数*/
//int n; /*整型变量，输入参数，给定数据点的个数*/
//int m; /*整型变量，输入参数，拟合多项式的项数*/
//double *dt1; /*实型变量，输出参数，拟合多项式与数据点偏差的平方和*/
//double *dt2; /*实型变量，输出参数，拟合多项式与数据点偏差的绝对值之和*/
//double *dt3; /*实型变量，输出参数，拟合多项式与数据点偏差的绝对值最大值*/
{
    int i ,j ,k ;
    double *s,*t,*b,z,d1,p,c,d2,g,q,dt;
    /*分别为s ,t ,b分配存贮空间*/
    s = (double *)calloc(n,sizeof(double));
    if(s == NULL)
    {
        printf("内存分配失败\n");
        exit (0);
    }
    t = (double *)calloc(n,sizeof(double));
    if(t == NULL)
    {
        printf("内存分配失败\n");
        exit (0);
    }
    b = (double *)calloc(n,sizeof(double));
    if(b == NULL)
    {
        printf("内存分配失败\n");
        exit (0);
    }
    z = 0;
    for(i=1;i<=n;i++)
        z=z+x[i-1]/n; /*z为各个x的平均值*/
    b[0]=1;
    d1=n;
    p=0;
    c=0;
    for(i=1;i<=n;i++)
    {
        p=p+x[i-1]-z;
        c=c+y[i-1];
    }
    c=c/d1;
    p=p/d1;
    a[0]=c*b[0];
    if(m>1)
    {
        t[1]=1;
        t[0]=-p;
        d2=0;
        c=0;
        g=0;
        for(i=1;i<=n;i++)
        {
            q=x[i-1]-z-p;
            d2=d2+q*q;
            c=y[i-1]*q+c;
            g=(x[i-1]-z)*q*q+g;
        }
        c=c/d2;
        p=g/d2;
        q=d2/d1;
        d1=d2;
        a[1]=c*t[1];
        a[0]=c*t[0]+a[0];
    }
    for(j=3;j<=m;j++)
    {
        s[j-1]=t[j-2];
        s[j-2]=-p*t[j-2]+t[j-3];
        if(j>=4)
            for(k=j-2;k>=2;k--)
                s[k-1]=-p*t[k-1]+t[k-2]-q*b[k-1];
        s[0]=-p*t[0]-q*b[0];
        d2=0;
        c=0;
        g=0;
        for(i=1;i<=n;i++)
        {
            q=s[j-1];
            for(k=j-1;k>=1;k--)
                q=q*(x[i-1]-z)+s[k-1]; 
            d2=d2+q*q; 
            c=y[i-1]*q+c; 
            g=(x[i-1]-z)*q*q+g; 
        } 
        c=c/d2; 
        p=g/d2; 
        q=d2/d1; 
        d1=d2; 
        a[j-1]=c*s[j-1]; 
        t[j-1]=s[j-1]; 
        for(k=j-1;k>=1;k--) 
        { 
            a[k-1]=c*s[k-1]+a[k-1]; 
            b[k-1]=t[k-1]; 
            t[k-1]=s[k-1]; 
        } 
    } 
    *dt1=0; 
    *dt2=0; 
    *dt3=0; 
    for(i=1;i<=n;i++) 
    { 
        q=a[m-1]; 
        for(k=m-1;k>=1;k--) 
            q=q*(x[i-1]-z)+a[k-1]; 
        dt=q-y[i-1]; 
        if(fabs(dt)>*dt3) 
            *dt3=fabs(dt); 
        *dt1=*dt1+dt*dt; 
        *dt2=*dt2+fabs(dt); 
    } 
    /*释放存储空间*/ 
    free(s); 
    free(t); 
    free(b); 
    return(1); 
}

void main()
{
    int i ,n ,m ;
    double *x,*y,*a,dt1,dt2,dt3,b;
    n = 12;// 12个样点
    m = 4; //3次多项式拟合
    b = 0; //x的初值为0
    /*分别为x,y,a分配存贮空间*/
    x = (double *)calloc(n,sizeof(double));
    if(x == NULL)
    {
        printf("内存分配失败\n");
        exit (0);
    }
    y = (double *)calloc(n,sizeof(double));
    if(y == NULL)
    {
        printf("内存分配失败\n");
        exit (0);
    }
    a = (double *)calloc(n,sizeof(double));
    if(a == NULL)
    {
        printf("内存分配失败\n");
        exit (0);
    }
    for(i=1;i<=n;i++)
    {
        x[i-1]=b+(i-1)*5;
        /*每隔5取一个点，这样连续取12个点*/
    }
    y[0]=0;
    y[1]=1.27;
    y[2]=2.16;
    y[3]=2.86;
    y[4]=3.44;
    y[5]=3.87;
    y[6]=4.15;
    y[7]=4.37;
    y[8]=4.51;
    y[9]=4.58;
    y[10]=4.02;
    y[11]=4.64;
    /*x[i-1]点对应的y值是拟合已知值*/
    
    Smooth(x,y,a,n,m,&dt1,&dt2,&dt3); /*调用拟合函数*/
    for(i=1;i<=m;i++)
        printf("a[%d] = %.10f\n",(i-1),a[i-1]);
    printf("拟合多项式与数据点偏差的平方和为：\n");
    printf("%.10e\n",dt1);
    printf("拟合多项式与数据点偏差的绝对值之和为：\n");
    printf("%.10e\n",dt2);
    printf("拟合多项式与数据点偏差的绝对值最大值为：\n");
    printf("%.10e\n",dt3);
    free(x); /*释放存储空间*/
    free(y); /*释放存储空间*/
    free(a); /*释放存储空间*/
}