數據處理與迴歸分析

#include "stdafx.h"
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include "LinearEquation.h"
#include "RegressAnalyze.h"

// 隨機樣本分析
void irhis(double x[], int n, double x0, double h, int m, int l, double dt[], int g[], int q[])
{
 int i,j,k,z;
 double s;
 char a[50];
 dt[0]=0.0;
 for (i=0; i<=n-1; i++)
  dt[0]=dt[0]+x[i]/n;
 dt[1]=0.0;
 for (i=0; i<=n-1; i++)
  dt[1]=dt[1]+(x[i]-dt[0])*(x[i]-dt[0]);
 dt[1]=dt[1]/n;
 dt[2]=sqrt(dt[1]);
 for (i=0; i<=m-1; i++)
 {
  q[i]=0;
  s=x0+(i+0.5)*h-dt[0];
  s=exp(-s*s/(2.0*dt[1]));
  g[i]=n*s*h/(dt[2]*2.5066);
 }
 s=x0+m*h;
 for (i=0; i<=n-1; i++)
  if ((x[i]-x0)>=0.0)
   if ((s-x[i])>=0.0)
   {
    j=(x[i]-x0)/h;
    q[j]=q[j]+1;
   }
 if (l==0)
  return;
 printf("\n");
 printf("n=%d\n",n);
 printf("\n");
 printf("x0=%e  h=%e  m=%d\n",x0,h,m);
 printf("\n");
 printf("xa=%e  s=%e  t=%e\n",dt[0],dt[1],dt[2]);
 printf("\n");
 k=1;
 z=0;
 for (i=0; i<=m-1; i++)
  if (q[i]>z) z=q[i];
 while (z>50)
 {
  z=z/2;
  k=2*k;
 }
 for (i=0; i<=m-1; i++)
 {
  s=x0+(i+0.5)*h;
  for (j=0; j<=49; j++)
   a[j]=' ';
  j=q[i]/k;
  for (z=0; z<=j-1; z++)
   a[z]='X';
  j=g[i]/k;
  if ((j>0)&&(j<=50))
   a[j]='*';
  printf("%e  %7d  ",s,q[i]);
  for (j=0; j<=49; j++)
   printf("%c",a[j]);
  printf("\n");
 }
 printf("\n");
 printf("k=%d\n",k);
 printf("\n");
 return;
}

//一元線性迴歸分析
void isqt1(double x[], double y[], int n, double a[2], double dt[6])
{
 int i;
 double xx,yy,e,f,q,u,p,umax,umin,s;
 xx=0.0;
 yy=0.0;
 for (i=0; i<=n-1; i++)
 {
  xx=xx+x[i]/n;
  yy=yy+y[i]/n;
 }
 e=0.0;
 f=0.0;
 for (i=0; i<=n-1; i++)
 {
  q=x[i]-xx;
  e=e+q*q;
  f=f+q*(y[i]-yy);
 }
 a[1]=f/e;
 a[0]=yy-a[1]*xx;
 q=0.0;
 u=0.0;
 p=0.0;
 umax=0.0;
 umin=1.0e+30;
 for (i=0; i<=n-1; i++)
 {
  s=a[1]*x[i]+a[0];
  q=q+(y[i]-s)*(y[i]-s);
  p=p+(s-yy)*(s-yy);
  e=fabs(y[i]-s);
  if (e>umax) umax=e;
  if (e<umin) umin=e;
  u=u+e/n;
 }
 dt[1]=sqrt(q/n);
 dt[0]=q;
 dt[2]=p;
 dt[3]=umax;
 dt[4]=umin;
 dt[5]=u;
 return;
}

//多元線性迴歸分析
void isqt2(double x[], double y[], int m, int n, double a[], double dt[], double v[])
{
 int i,j,k,mm;
 double q,e,u,p,yy,s,r,pp,*b;
 b=(double *)malloc((m+1)*(m+1)*sizeof(double));
 mm=m+1;
 b[mm*mm-1]=n;
 for (j=0; j<=m-1; j++)
 {
  p=0.0;
  for (i=0; i<=n-1; i++)
   p=p+x[j*n+i];
  b[m*mm+j]=p;
  b[j*mm+m]=p;
 }
 for (i=0; i<=m-1; i++)
  for (j=i; j<=m-1; j++)
  {
   p=0.0;
   for (k=0; k<=n-1; k++)
    p=p+x[i*n+k]*x[j*n+k];
   b[j*mm+i]=p;
   b[i*mm+j]=p;
  }
 a[m]=0.0;
 for (i=0; i<=n-1; i++)
  a[m]=a[m]+y[i];
 for (i=0; i<=m-1; i++)
 {
  a[i]=0.0;
  for (j=0; j<=n-1; j++)
   a[i]=a[i]+x[i*n+j]*y[j];
 }
 achol(b,mm,1,a);
 yy=0.0;
 for (i=0; i<=n-1; i++)
  yy=yy+y[i]/n;
 q=0.0;
 e=0.0;
 u=0.0;
 for (i=0; i<=n-1; i++)
 {
  p=a[m];
  for (j=0; j<=m-1; j++)
   p=p+a[j]*x[j*n+i];
  q=q+(y[i]-p)*(y[i]-p);
  e=e+(y[i]-yy)*(y[i]-yy);
  u=u+(yy-p)*(yy-p);
 }
 s=sqrt(q/n);
 r=sqrt(1.0-q/e);
 for (j=0; j<=m-1; j++)
 {
  p=0.0;
  for (i=0; i<=n-1; i++)
  {
   pp=a[m];
   for (k=0; k<=m-1; k++)
    if (k!=j)
     pp=pp+a[k]*x[k*n+i];
   p=p+(y[i]-pp)*(y[i]-pp);
  }
  v[j]=sqrt(1.0-q/p);
 }
 dt[0]=q;
 dt[1]=s;
 dt[2]=r;
 dt[3]=u;
 free(b);
 return;
}

//逐步迴歸分析
void isqt3(int n, int k, double x[], double f1, double f2, double eps, double xx[], double b[], double v[], double s[], \
     double dt[], double ye[], double yr[], double 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;
}

//半對數數據相關
void ilog1(int n, double x[], double y[], double t, double a[7])
{
 int i;
 double xx,yy,dx,dxy;
 xx=0.0;
 yy=0.0;
 for (i=0; i<=n-1; i++)
 {
  xx=xx+x[i]/n;
  yy=yy+log(y[i])/log(t)/n;
 }
 dx=0.0;
 dxy=0.0;
 for (i=0; i<=n-1; i++)
 {
  a[2]=x[i]-xx;
  dx=dx+a[2]*a[2];
  dxy=dxy+a[2]*(log(y[i])/log(t)-yy);
 }
 a[1]=dxy/dx;
 a[0]=yy-a[1]*xx;
 a[0]=a[0]*log(t);
 a[0]=exp(a[0]);
 a[2]=0.0;
 a[6]=0.0;
 a[4]=0.0;
 a[5]=1.0e+30;
 for (i=0; i<=n-1; i++)
 {
  a[3]=a[1]*x[i]*log(t);
  a[3]=a[0]*exp(a[3]);
  a[2]=a[2]+(y[i]-a[3])*(y[i]-a[3]);
  dx=fabs(y[i]-a[3]);
  if (dx>a[4])
   a[4]=dx;
  if (dx<a[5])
   a[5]=dx;
  a[6]=a[6]+dx/n;
 }
 a[3]=sqrt(a[2]/n);
 return;
}

//對數數據相關
void ilog2(int n, double x[], double y[], double a[7])
{
 int i;
 double xx,yy,dx,dxy;
 xx=0.0;
 yy=0.0;
 for (i=0; i<=n-1; i++)
 {
  xx=xx+log(x[i])/n;
  yy=yy+log(y[i])/n;
 }
 dx=0.0;
 dxy=0.0;
 for (i=0; i<=n-1; i++)
 {
  a[2]=log(x[i])-xx;
  dx=dx+a[2]*a[2];
  dxy=dxy+a[2]*(log(y[i])-yy);
 }
 a[1]=dxy/dx;
 a[0]=yy-a[1]*xx;
 a[0]=exp(a[0]);
 a[2]=0.0;
 a[6]=0.0;
 a[4]=0.0;
 a[5]=1.0e+30;
 for (i=0; i<=n-1; i++)
 {
  a[3]=a[1]*log(x[i]);
  a[3]=a[0]*exp(a[3]);
  a[2]=a[2]+(y[i]-a[3])*(y[i]-a[3]);
  dx=fabs(y[i]-a[3]);
  if (dx>a[4])
   a[4]=dx;
  if (dx<a[5])
   a[5]=dx;
  a[6]=a[6]+dx/n;
 }
 a[3]=sqrt(a[2]/n);
 return;
}
                                      ----根據《C語言經常使用算法程序集》整理

<投票>