aic.cpp

估计模型出阶次的AIC法的VC程序

#include <iostream.h>
#include <math.h>
#include <stdlib.h>
#include<fstream.h>  

//矩阵求逆函数

  int brinv(double f[],int n)
  { int *is,*js,i,j,k,l,u,v;
    double d,p;
    is=(int *)malloc(n*sizeof(int));
    js=(int *)malloc(n*sizeof(int));
    for (k=0; k<=n-1; k++)
      { d=0.0;
        for (i=k; i<=n-1; i++)
        for (j=k; j<=n-1; j++)
          { l=i*n+j; p=fabs(f[l]);
            if (p>d) { d=p; is[k]=i; js[k]=j;}
          }
        if (d+1.0==1.0)
          { free(is); free(js); cout<<"err**not inv\n";
            return(0);
          }
        if (is[k]!=k)
          for (j=0; j<=n-1; j++)
            { u=k*n+j; v=is[k]*n+j;
              p=f[u]; f[u]=f[v]; f[v]=p;
            }
        if (js[k]!=k)
          for (i=0; i<=n-1; i++)
            { u=i*n+k; v=i*n+js[k];
              p=f[u]; f[u]=f[v]; f[v]=p;
            }
        l=k*n+k;
        f[l]=1.0/f[l];
        for (j=0; j<=n-1; j++)
          if (j!=k)
            { u=k*n+j; f[u]=f[u]*f[l];}
        for (i=0; i<=n-1; i++)
          if (i!=k)
            for (j=0; j<=n-1; j++)
              if (j!=k)
                { u=i*n+j;
                  f[u]=f[u]-f[i*n+k]*f[k*n+j];
                }
        for (i=0; i<=n-1; i++)
          if (i!=k)
            { u=i*n+k; f[u]=-f[u]*f[l];}
      }
    for (k=n-1; k>=0; k--)
      { if (js[k]!=k)
          for (j=0; j<=n-1; j++)
            { u=k*n+j; v=js[k]*n+j;
              p=f[u]; f[u]=f[v]; f[v]=p;
            }
        if (is[k]!=k)
          for (i=0; i<=n-1; i++)
            { u=i*n+k; v=i*n+is[k];
              p=f[u]; f[u]=f[v]; f[v]=p;
            }
      }
    free(is); free(js);
    return(1);
  }
  //矩阵相乘函数
 int brmul(double a[],double b[],int m,int n,int k,double c[])
    { int i,j,l,u;
    for (i=0; i<=m-1; i++)
    for (j=0; j<=k-1; j++)
      { u=i*k+j; c[u]=0.0;
        for (l=0; l<=n-1; l++)
          c[u]=c[u]+a[i*n+l]*b[l*k+j];
      }
    return 0;
  }

//主函数
void main(){
     int i,k,u[1000],j,k1;
     double z[1000],v[1000],z1[600],q11[7200],q1[7200],qn[144],q3[7200],r[12],B[4];
	 double q24[24],r12[2],q29[2],A[24],q21[1200],q2[1200],q22[24],q25[7200],q221[4],q23[1200],q231[4],q26[600],q27[600],q28[12],r11[12],r0[14],q[8400],r01[600],w[601],aic,b;
     //double para[10]={-1.185,0.814,-0.518,0.349,-0.117,1.08,-0.745,0.475,-0.253,0.123}
     
	 ofstream fop("data.txt");
     ifstream fip1("m.txt");
     ifstream fip2("Gauss.txt");
	 
     for(i=0;i<1000;i++){
	    fip1>>u[i];
	    fip2>>v[i];	
	 }
	 
     z[0]=v[0];
     z[1]=1.185*z[0]+1.08*u[0]+v[1];
     z[2]=1.185*z[1]-0.814*z[0]+1.08*u[1]-0.745*u[0]+v[2];
     z[3]=1.185*z[2]-0.814*z[1]+0.518*z[0]+1.08*u[2]-0.745*u[1]+0.475*u[0]+v[3];
     z[4]=1.185*z[3]-0.814*z[2]+0.518*z[1]-0.349*z[0]+1.08*u[3]-0.745*u[2]+0.475*u[1]-0.253*u[0]+v[4];
     for(i=5;i<1000;i++)
	 {
	    z[i]=1.185*z[i-1]-0.814*z[i-2]+0.518*z[i-3]-0.349*z[i-4]+0.117*z[i-5]+1.08*u[i-1]-0.745*u[i-2]+0.475*u[i-3]-0.253*u[i-4]+0.123*u[i-5]+v[i];
		
	 }
     i=0;
     for(k=0;k<600;k++)
	 {
	     q2[2*i]=-z[k];
		 q2[2*i+1]=u[k];
         i=i+1;
		 
	 }
	 
	 i=0;
     for(j=0;j<600;j++)
	 {
		q21[i]=-z[j];
		i=i+1;
	 }
     for(j=0;j<600;j++)
	 {
		q21[i]=u[j];
		i=i+1;
        
	 }
	 //for(i=0;i<200;i++)
		 //cout<<q21[i]<<endl;
     for(k1=1;k1<7;k1++)
	 {
        i=0;
	    for(k=0;k<600;k++){
	        for(j=0;j<k1;j++){
	           q[i]=-z[k+k1-1-j];
               i=i+1;
			}
	       for(j=0;j<k1;j++){
		      q[i]=u[k+k1-1-j];
		      i=i+1;
		   }
		}
        i=0;
	    for(k=0;k<600;k++){
           z1[k]=z[k+k1];
	       for(j=0;j<k1-1;j++){
		      q1[2*i]=-z[k1-1+k-j];
		      q1[2*i+1]=u[k1-1+k-j];              
		      i=i+1;
		   }
		}
        
		 
	 i=0;
     for(k=0;k<k1-1;k++)
	 {
		 for(j=0;j<600;j++)
		 {
			 q11[i]=-z[j+k1-1-k];
			 i=i+1;
		 }
		 for(j=0;j<600;j++)
		 {
             q11[i]=u[j+k1-1-k];
			 i=i+1;
		 }
	 }
	 
      brmul(q11,q1,2*k1-2,600,2*k1-2,qn);
      brinv(qn,2*k1-2);
      brmul(qn,q11,2*k1-2,2*k1-2,600,q3);
      brmul(q3,z1,2*k1-2,600,1,r);
      brmul(q21,q2,2,600,2,q221);
	  brmul(q21,q1,2,600,2*k1-2,q22);
      brmul(q22,q3,2,2*k1-2,600,q23);
      brmul(q23,q2,2,600,2,q231);
	  for(i=0;i<4;i++)
		  B[i]=q221[i]-q231[i];
	  brinv(B,2);
      brmul(q3,q2,2*k1-2,600,2,q24);
	  brmul(q24,B,2*k1-2,2,2,A);
      brmul(A,q21,2*k1-2,2,600,q25);
      brmul(q1,r,600,2*k1-2,1,q26);
	  for(i=0;i<600;i++){
		  q27[i]=z1[i]-q26[i];
	  }
      brmul(q25,q27,2*k1-2,600,1,q28);
	  for(i=0;i<2*k1-2;i++)
	  {
		  r11[i]=r[i]-q28[i];
	  }
      brmul(q21,q27,2,600,1,q29);
      brmul(B,q29,2,2,1,r12);
	  j=0;
	  for(i=0;i<k1-1;i++)
	  {
		  r0[j]=r11[2*i];
		  j=j+1;
	  }
	  j=j+1;
	  for(i=0;i<k1-1;i++)
	  {
		  r0[j]=r11[2*i+1];
		  j=j+1;
	  }
	  r0[k1-1]=r12[0];
	  r0[2*k1-1]=r12[1];
	  
	  
     brmul(q,r0,600,2*k1,1,r01);
	 w[0]=0;
	 for(j=0;j<600;j++)
	 {
		w[j+1]=w[j]+pow(z1[j]-r01[j],2);
	 }
	  b=w[600]/600;
	  aic=600*(log(b))+4*k1;

      for(i=0;i<2*k1;i++)
	  {cout<<r0[i]<<endl;fop<<r0[i]<<endl;}
         cout<<k1<<endl;
		 cout<<"判断阶次的值"<<"  "<<aic<<endl;
		 fop<<k1<<endl;
		 fop<<"判断的阶次的值"<<"  "<<aic<<endl;	  
	  
	 }
	 fop.close();
}