ch9.h

将C语言的常用程序集移植到VC开发环境的源代码

/************************************************
 Expect bugs!
 Please use and enjoy, and let me know of any bugs/mods/improvements 
 that you have found/implemented and I will fix/incorporate them into 
 this file. Thank Mr. Xushiliang once again.

					hujinshan@2002.11.3
				Airforce Engineering University
************************************************/

/*****  #include "CH9.h"  数据处理与回归分析*****/
#ifndef CH9_H_
#define CH9_H_

#include "stdlib.h"
#include "math.h"
#include "stdio.h"

#include "ch1.h"
//*******************************************************************

void irhis(double x[],int n,double x0,int h,int m,int l,double dt[],int g[],int q[],char str[])//随机样本分析
void isqt1(double x[],double y[],int n,double a[2],double dt[6]);//一元回归分析
void isqt2(double x[],double y[],int m,int n,double a[],double dt[],double v[]);//多元回归分析
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[]);//逐步回归分析
void ilog1(int n,double x[],double y[],double t,double a[7]);//半对数数据相关
void ilog2(int n,double x[],double y[],double a[]);//对数数据相关

//*******************************************************************
void irhis(double x[],int n,double x0,int h,int m,int l,double dt[],int g[],int q[],char str[])
{ 
	int i,j,k,z;
    double s;
    char a[100];
    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;
	//*****************
	char buffer[100];
	
	strcat( str, "\r\n" );
	sprintf( buffer,"n=%d\r\n",n);
	strcat( str, buffer );

	sprintf( buffer,"x0=%e  h=%d  m=%d \r\n",x0,h,m);
    strcat( str, buffer );

    sprintf( buffer,"xa=%e  s=%e  t=%e\r\n",dt[0],dt[1],dt[2]);
    strcat( str, buffer );

    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]='*';

		sprintf( buffer,"%e  %7d  ",s,q[i]);
		strcat( str, buffer );

        for (j=0; j<=49; j++)
		{
			sprintf( buffer,"%c",a[j]);
			strcat( str, buffer );
		}
		strcat( str, "\r\n" );
      }
	sprintf( buffer,"k=%d\r\n",k);
	strcat( str, buffer );
    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[])
{ 
	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;
}
#endif
/////////////////////////////////////////////////////////////
/*
void irhis(double x[],int n,double x0,int h,int m,int l,double dt[],int g[],int q[],CString& strRet)
{
	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;
}
/**/