AIC法定阶的依阶次递推算法程序.

AIC 法定阶的依阶次递推算法程序

依阶次递推算法所得到估计θ

ˆ，再按下式计算残差方差的估计值： ∑==N j j e N 122

)ˆ,(1ˆθσ 由上式的结果计算AIC ：

)(2ˆlg AIC 2b a n n N ++=σ

在结果中找到AIC 最小的模型（阶次和参数）就是估计的模型。

由输出数据可知当k1=5时aic 的值最小。所以最后的辨识结果取阶次为5，参数为：

–1.18394，0.813938，–0.518174，0.348744，–0.116818，

1.07998，–0.74386，0.475444，–0.253022，0.122781

判断的阶次的最小aic 值： aic= – 8981.58发生在阶次为5时。

源程序：

#include

#include

#include

#include

//矩阵求逆函数

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;

}