main.c

这是一个关于逐步回归求多元线性回归的一个VC例程

#include "math.h"
  void sqt3(n,k,x,f1,f2,eps,xx,b,v,s,dt,ye,yr,r)
  int n,k;
  double f1,f2,eps,x[],xx[],b[],v[],s[],dt[],ye[],yr[],r[];
  { int i,j,ii,m,imi,imx,l,it;
    double z,phi,sd,vmi,vmx,q,fmi,fmx;
    m=n+1; q=0.0;
    for (j=0; j<=n; j++)
      { z=0.0;
        for (i=0; i<=k-1; i++)
          z=z+x[i*m+j]/k;
        xx[j]=z;
      }
    for (i=0; i<=n; i++)
      for (j=0; j<=i; j++)
        { z=0.0;
          for (ii=0; ii<=k-1; ii++)
            z=z+(x[ii*m+i]-xx[i])*(x[ii*m+j]-xx[j]);
          r[i*m+j]=z;
        }
    for (i=0; i<=n; i++)
      ye[i]=sqrt(r[i*m+i]);
    for (i=0; i<=n; i++)
      for (j=0; j<=i; j++)
        { r[i*m+j]=r[i*m+j]/(ye[i]*ye[j]);
          r[j*m+i]=r[i*m+j];
        }
    phi=k-1.0;
    sd=ye[n]/sqrt(k-1.0);
    it=1;
    while (it==1)
      { it=0;
        vmi=1.0e+35; vmx=0.0;
        imi=-1; imx=-1;
        for (i=0; i<=n; i++)
          { v[i]=0.0; b[i]=0.0; s[i]=0.0;}
        for (i=0; i<=n-1; i++)
          if (r[i*m+i]>=eps)
            { v[i]=r[i*m+n]*r[n*m+i]/r[i*m+i];
              if (v[i]>=0.0)
                { if (v[i]>vmx)
                    { vmx=v[i]; imx=i;}
                }
              else
                { b[i]=r[i*m+n]*ye[n]/ye[i];
                  s[i]=sqrt(r[i*m+i])*sd/ye[i];
                  if (fabs(v[i])<vmi)
                    { vmi=fabs(v[i]); imi=i;}
                }
            }
        if (phi!=n-1.0)
          { z=0.0;
            for (i=0; i<=n-1; i++)
              z=z+b[i]*xx[i];
            b[n]=xx[n]-z; s[n]=sd; v[n]=q;
          }
        else
          { b[n]=xx[n]; s[n]=sd;}
        fmi=vmi*phi/r[n*m+n];
        fmx=(phi-1.0)*vmx/(r[n*m+n]-vmx);
        if ((fmi<f2)||(fmx>=f1))
          { if (fmi<f2)
              { phi=phi+1.0; l=imi;}
            else
              { phi=phi-1.0; l=imx;}
            for (i=0; i<=n; i++)
              if (i!=l)
                for (j=0; j<=n; j++)
                  if (j!=l)
                    r[i*m+j]=r[i*m+j]-
                             (r[l*m+j]/r[l*m+l])*r[i*m+l];
            for (j=0; j<=n; j++)
              if (j!=l)
                r[l*m+j]=r[l*m+j]/r[l*m+l];
            for (i=0; i<=n; i++)
              if (i!=l)
                r[i*m+l]=-r[i*m+l]/r[l*m+l];
            r[l*m+l]=1.0/r[l*m+l];
            q=r[n*m+n]*ye[n]*ye[n];
            sd=sqrt(r[n*m+n]/phi)*ye[n];
            dt[0]=sqrt(1.0-r[n*m+n]);
            dt[1]=(phi*(1.0-r[n*m+n]))/
          ((k-phi-1.0)*r[n*m+n]);
            it=1;
          }
      }
    for (i=0; i<=k-1; i++)
      { z=0.0;
        for (j=0; j<=n-1; j++)
          z=z+b[j]*x[i*m+j];
        ye[i]=b[n]+z; yr[i]=x[i*m+n]-ye[i];
      }
    return;
  }
#include "stdio.h"
// #include "11sqt3.c"
void  main()
  { int i,j,k;
    double eps,xx[5],b[5],v[5],s[5],ye[13],yr[13];
    double r[5][5],dt[2];
    static double x[13][5]={ 
                  {7.0,26.0,6.0,60.0,78.5},
                  {1.0,29.0,15.0,52.0,74.3},
                  {11.0,56.0,8.0,20.0,104.3},
                  {11.0,31.0,8.0,47.0,87.6},
                  {7.0,52.0,6.0,33.0,95.9},
                  {11.0,55.0,9.0,22.0,109.2},
                  {3.0,71.0,17.0,6.0,102.7},
                  {1.0,31.0,22.0,44.0,72.5},
                  {2.0,54.0,18.0,22.0,93.1},
                  {21.0,47.0,4.0,26.0,115.9},
                  {1.0,40.0,23.0,34.0,83.8},
                  {11.0,66.0,9.0,12.0,113.3},
                  {10.0,68.0,8.0,12.0,109.4}};
    static double f1[3]={1.46,4.75,9.33};
    static double f2[3]={1.45,4.67,9.07};
    eps=1.0e-30;
    for (k=0; k<=2; k++)
      { sqt3(4,13,x,f1[k],f2[k],eps,xx,b,v,s,dt,ye,yr,r);
        printf("\n");
        printf("f1=%e  f2=%e\n",f1[k],f2[k]);
        printf("\n");
        printf("original x(i) and y values:\n");
        for (i=0; i<=12; i++)
          { for (j=0; j<=3; j++)
              printf("x(%d)=%6.2f ",j,x[i][j]);
            printf("y=%6.2f\n",x[i][4]);
          }
        printf("\n");
        printf("mean of x(i) and y:\n");
        for (i=0; i<=3; i++)
          printf("x(%d)=%6.3f ",i,xx[i]);
        printf("y=%6.3f\n",xx[4]);
        printf("\n");
        printf("regression coeffi b(i):\n");
        for (i=0; i<=4; i++)
          printf("b(%d)=%9.6f ",i,b[i]);
        printf("\n\n");
        printf("standard partial sum of square of \n");
        printf("regression for x(i) and sum of \n");
        printf("square of residuals:\n");
        for (i=0; i<=3; i++)
          printf("v(%d)=%e ",i,v[i]);
        printf("\n");
        printf("q=%e\n",v[4]);
        printf("\n");
        printf("standard deviation of regression  \n");
        printf("coeffi and regression equation:\n");
        for (i=0; i<=3; i++)
          printf("s(%d)=%e ",i,s[i]);
        printf("\n");
        printf("s=%e\n",s[4]);
        printf("\n");
        printf("multi-correlation coeffi c is:%e\n",dt[0]);
        printf("\n");
        printf("the f value=%e\n",dt[1]);
        printf("\n");
        printf("estimated values and residuals:\n");
        for (i=0; i<=12; i++)
          printf("ye(%d)=%e  yr(%d)=%e\n",
                  i,ye[i],i,yr[i]);
        printf("\n");
        printf("matrix r:\n");
        for (i=0; i<=4; i++)
          { for (j=0; j<=4; j++)
              printf("%e  ",r[i][j]);
            printf("\n");
          }
        printf("\n");
      }
  }