📄 aicdlg.cpp
字号:
//求损失函数、并保存此时的阶次、估计参数sita值。
brmul(e1,e1,1,N,1,sigama);
v[0]=sigama[0];
na0=1;nb0=1;
memcpy(sita0,sita1,2*sizeof(double));
//
//na=nb=n; X=[X1，X2]
for(n=2;n<=10;n++)
{
na=n;nb=n;
//模型阶次增加，重新赋值给X2T矩阵
for(i=0;i<2;i++)
{
for(j=0;j<N;j++)
{
if(i==0)
{
X2T[i*N+j]=-z[10-n+j];
}
else
{
X2T[i*N+j]=u[10-n+j];
}
}
}
//X2T矩阵转置得X2矩阵
for(i=0;i<2;i++)
{
for(j=0;j<N;j++)
{
X2[j*2+i]=X2T[i*N+j];
}
}
X1TX2=(double*)calloc(2*(n-1)*2,sizeof(double));
X2TX1=(double*)calloc(2*2*(n-1),sizeof(double));
B1=(double*)calloc(2*(n-1)*2,sizeof(double));
AA=(double*)calloc(2*(n-1)*2,sizeof(double));
AX2T=(double*)calloc(2*(n-1)*N,sizeof(double));
sita_temp=(double*)calloc(2*(n-1),sizeof(double));
//阶次增加，重新计算增加个数后的所有参数估计值并保存。
//sita1=(sita1*)-A*X2T*(Y-X1*(sita1*))
//sita2=B*X2T*(Y-X1*(sita1*))
//A=AA=inv(X1TX1)*(X1TX2)*B
//B=BB=inv(X2TX2-X2TX1*inv(X1TX1)*X1TX2)
brmul(X2T,X1,2,N,2*(n-1),X2TX1);
brmul(X1T,X2,2*(n-1),N,2,X1TX2);
brmul(X2T,X2,2,N,2,X2TX2);
brmul(X1TX1,X1TX2,2*(n-1),2*(n-1),2,B1);
brmul(X2TX1,B1,2,2*(n-1),2,B2);
for(i=0;i<4;i++)
{
BB[i]=X2TX2[i]-B2[i];
}
brinv(BB,2);
brmul(B1,BB,2*(n-1),2,2,AA);
brmul(AA,X2T,2*(n-1),2,N,AX2T);
brmul(BB,X2T,2,2,N,BX2T);
brmul(AX2T,e1,2*(n-1),N,1,sita_temp);
for(i=0;i<2*(n-1);i++)
{
sita1[i]=sita1[i]-sita_temp[i];
}
brmul(BX2T,e1,2,N,1,sita2);
memcpy(&sita1[2*(n-1)],sita2,2*sizeof(double));
//把X2补充到X1里，形成新的X1矩阵，重新计算增加阶次后的参数估计误差
memcpy(&X1T[2*(n-1)*N],X2T,2*N*sizeof(double));
//将X1T矩阵转置得X1矩阵。
free(X1);
X1=(double*)calloc(N*2*n,sizeof(double));
for(i=0;i<2*n;i++)
{
for(j=0;j<N;j++)
{
X1[j*2*n+i]=X1T[i*N+j]; //将X1T矩阵转置得X1矩阵。
}
}
brmul(X1,sita1,N,2*n,1,Ye);
for(i=0;i<N;i++)
{
e1[i]=Y[i]-Ye[i];
}
brmul(e1,e1,1,N,1,sigama);
//利用AIC法得模型参数估计值及模型阶次。
v[n-1]=sigama[0];
t[n-1]=(N-2*(n-1))*(v[n-2]-v[n-1])/(2*v[n-1]);
if(t[n-1]>3)
{
t[0]=t[n-1];
na0=na;nb0=nb;
memcpy(sita0,sita1,2*n*sizeof(double));
}
//////////////////////////////////////////
//重新构造XTX的逆阵，提供下一次递推用。
XTX_inv=(double*)calloc(2*n*2*n,sizeof(double));
CC=(double*)calloc(2*(n-1)*2*(n-1),sizeof(double));
temp1=(double*)calloc(2*(n-1)*2*(n-1),sizeof(double));
temp2=(double*)calloc(2*(n-1)*2*(n-1),sizeof(double));
brmul(AA,X2TX1,2*(n-1),2,2*(n-1),temp1);
brmul(temp1,X1TX1,2*(n-1),2*(n-1),2*(n-1),temp2);
for(i=0;i<2*(n-1);i++)
{
for(j=0;j<2*(n-1);j++)
{
CC[i*2*(n-1)+j]=X1TX1[i*2*(n-1)+j]+temp2[i*2*(n-1)+j];
}
}
for(i=0;i<2*(n-1);i++)
{
for(j=0;j<2*(n-1);j++)
{
XTX_inv[i*2*n+j]=CC[i*2*(n-1)+j];
}
}
for(j=0;j<2*(n-1);j++)
{
XTX_inv[2*(n-1)*2*n+j]=-AA[2*j];
XTX_inv[(2*n-1)*2*n+j]=-AA[2*j+1];
}
for(i=0;i<2*(n-1);i++)
{
for(j=2*(n-1);j<2*n;j++)
{
XTX_inv[i*2*n+j]=-AA[i*2+j-2*(n-1)];
}
}
for(i=2*(n-1);i<2*n;i++)
{
for(j=2*(n-1);j<2*n;j++)
{
XTX_inv[i*2*n+j]=BB[(i-2*(n-1))*2+j-2*(n-1)];
}
}
//
free(temp1);
free(temp2);
free(CC);
//
free(X1TX1);
X1TX1=(double*)calloc(2*n*2*n,sizeof(double));
memcpy(X1TX1,XTX_inv,2*n*2*n*sizeof(double));
////////////////////////////////////////////////
free(X1TX2);
free(X2TX1);
free(B1);
free(AA);
free(AX2T);
free(sita_temp);
free(XTX_inv);
}
free(X1TX1);free(X1);
free(z);free(u);free(e);free(X1T);free(X2T);free(X2);
free(Y);free(Ye);free(e1);free(BX2T);
WriteData2File("t_f",t,10,1);
if((fp4=fopen("sita_f.txt","w"))==NULL) return;
fprintf(fp4,"AIC最终辨识结果sita[]为：\n");
fprintf(fp4,"\nna=%d\n",na0);
fprintf(fp4,"nb=%d\n\n",nb0);
fprintf(fp4," a[i] b[i]\n");
for(i=0;i<na0;i++)
{
fprintf(fp4,"a[%d]=%10.6f ",i+1,sita0[i*2]);
fprintf(fp4,"b[%d]=%10.6f\n",i+1,sita0[i*2+1]);
}
fprintf(fp4,"\n");
for(i=1;i<10;i++)
{
fprintf(fp4,"t%d= %10.6f\n",i,t[i]);
}
fclose(fp4);
sita1[0]=-1.185;sita1[2]= 0.814;sita1[4]=-0.518;sita1[6]= 0.349;sita1[8]=-0.117;
sita1[1]= 1.08 ;sita1[3]=-0.745;sita1[5]= 0.475;sita1[7]=-0.253;sita1[9]= 0.123;
CString cs;
cs.Format("模型参数：");
for(j=0;j<2;j++)
{
for(i=0;i<na0;i++)
{
if(j==0)
{
cs.Format(cs+"\na%d=%.3f",i+1,sita1[j+2*i]);
}
else
{
cs.Format(cs+"\nb%d=%.3f",i+1,sita1[j+2*i]);
}
}
}
m_cStaticTextREFm.SetDrawText(cs);
cs.Format("辨识参数：");
for(j=0;j<2;j++)
{
for(i=0;i<na0;i++)
{
if(j==0)
{
cs.Format(cs+"\na%d=%.6f",i+1,sita0[j+2*i]);
}
else
{
cs.Format(cs+"\nb%d=%.6f",i+1,sita0[j+2*i]);
}
}
}
m_cStaticTextREFi.SetDrawText(cs);
cs.Format("辨识误差：");
for(j=0;j<2;j++)
{
for(i=0;i<na0;i++)
{
if(j==0)
{
cs.Format(cs+"\na%d=%.6f",i+1,fabs(sita0[j+2*i]-sita1[j+2*i]));
}
else
{
cs.Format(cs+"\nb%d=%.6f",i+1,fabs(sita0[j+2*i]-sita1[j+2*i]));
}
}
}
m_cStaticTextREFe.SetDrawText(cs);
delete m_tempInt;
m_tempInt=new int[10];
for(i=0;i<10;i++)
{
m_tempInt[i]=int(v[i]*20000);
}
m_cStaticGraph.SetDrawData(10,m_tempInt,1,1);
m_cStaticGraph.RedrawWindow();
RedrawWindow();
}
void CAICDlg::OnlyAic(int n, int lu)
{
//AIC实时递推算法
//本程序最大能算到5阶,10个参数。
int nn=5;
int L=n-10;
FILE *fp1,*fp2,*fp4;
int i,j,na,nb,na1,nb1;
double *w,*sita,*sita0,*sita1,*AIC;
w=(double*)calloc(nn*2,sizeof(double));
sita=(double*)calloc(nn*2,sizeof(double));
sita0=(double*)calloc(nn*2,sizeof(double));
sita1=(double*)calloc(nn*2,sizeof(double));
AIC=(double*)calloc(nn*nn,sizeof(double));
double *k,*k1,*k3,*p,*p0,*I;
k=(double*)calloc(nn*2,sizeof(double));
k1=(double*)calloc(nn*2,sizeof(double));
k3=(double*)calloc(nn*2*nn*2,sizeof(double));
p=(double*)calloc(nn*2*nn*2,sizeof(double));
p0=(double*)calloc(nn*2*nn*2,sizeof(double));
I=(double*)calloc(nn*2*nn*2,sizeof(double));
double sigama[1],k2[1],AIC1;
double *z,*u,*e,*ww,*z1,*e1;
z=(double*)calloc(L+10,sizeof(double));
u=(double*)calloc(L+lu+10,sizeof(double));
e=(double*)calloc(L+10,sizeof(double));
z1=(double*)calloc(L+10,sizeof(double));
e1=(double*)calloc(L+10,sizeof(double));
AIC1=99999;
if((fp1=fopen("MSerials.txt","r"))==NULL)return;
if((fp2=fopen("WhiteNoise2.txt","r"))==NULL)return;
//读入输入u[]和白噪声e[]
for(i=0;i<L+10+lu;i++)
{
fscanf(fp1,"%lf ",&u[i]);
}
for(i=0;i<L+10;i++)
{
u[i]=u[i+lu];
}
for(i=0;i<L+10;i++)
{
fscanf(fp2,"%lf ",&e[i]);
}
fclose(fp1);
fclose(fp2);
//用下式计算理论输出z[]
//z(k)-1.185z(k-1)+0.814z(k-2)-0.518z(k-3)+0.349z(k-4)-0.117z(k-5)
//=1.08u(k-1)-0.745u(k-2)+0.475u(k-3)-0.253u(k-4)+0.123u(k-5)+e(k)
z[0]=e[0];
z[1]=1.185*z[0]+1.08*u[0]+e[1];
z[2]=1.185*z[1]-0.814*z[0]+1.08*u[1]-0.745*u[0]+e[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]+e[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]+e[4];
for(i=5;i<L+10;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]+e[i];
}
n=0;
//给P和O设定初值
for(na=1;na<=nn;na++)
{
for(nb=1;nb<=nn;nb++)
{
for(i=0;i<na+nb;i++)
{
sita0[i]=0.0;
}
for(i=0;i<na+nb;i++)
{
for(j=0;j<na+nb;j++)
{
if(j==i)
{
I[i*(na+nb)+j]=1.0;
p0[i*(na+nb)+j]=1.0e+6;
}
else
{
I[i*(na+nb)+j]=0.0;
p0[i*(na+nb)+j]=0.0;
}
}
}
//给W矩阵赋值
ww=(double*)calloc(L*(na+nb),sizeof(double));
for(i=0;i<L;i++)
{
for(j=0;j<na;j++)
{
ww[i*(na+nb)+j]=-z[nn+i-j];
}
for(j=na;j<na+nb;j++)
{
ww[i*(na+nb)+j]=u[nn+nn+i-j];
}
}
for(i=0;i<L;i++)
{
for(j=0;j<na;j++)
{
w[j]=-z[nn+i-j];
}
for(j=na;j<na+nb;j++)
{
w[j]=u[nn+nn+i-j];
}
//计算方差并计算AIC
brmul(w,p0,1,na+nb,na+nb,k1);
brmul(k1,w,1,na+nb,1,k2);
k2[0]=1/(k2[0]+1);
brmul(p0,w,na+nb,na+nb,1,k1);
brmul(k1,k2,na+nb,1,1,k);
brmul(k,w,na+nb,1,na+nb,k3);
for(j=0;j<(na+nb)*(na+nb);j++)
{
k3[j]=I[j]-k3[j];
}
brmul(k3,p0,na+nb,na+nb,na+nb,p);
for(j=0;j<(na+nb)*(na+nb);j++)
{
p0[j]=p[j];
}
brmul(w,sita0,1,na+nb,1,k2);
k2[0]=z[i+5+1]-k2[0];
brmul(k,k2,na+nb,1,1,k1);
for(j=0;j<na+nb;j++)
{
sita[j]=sita0[j]+k1[j];
sita0[j]=sita[j];
}
}
brmul(ww,sita,L,na+nb,1,z1);
free(ww);
for(i=0;i<L;i++)
{
e1[i]=z[nn+1+i]-z1[i];
}
brmul(e1,e1,1,L,1,sigama);
AIC[n]=L+L*log10(sigama[0]/L)+2*(na+nb);
if(AIC1>AIC[n])
{
AIC1=AIC[n];
na1=na;
nb1=nb;
for(i=0;i<na+nb;i++)
{
sita1[i]=sita[i];
}
}
n++;
}
}
WriteData2File("aic",AIC,nn,nn);
WriteData2File("u",u,L+10,1);
WriteData2File("e",e,L+10,1);
WriteData2File("z",z,L+10,1);
if((fp1=fopen("aic1.txt","w"))==NULL)return;
fprintf(fp1,"AIC值结果AIC[]为：\n");
fprintf(fp1,"na=%d\nnb=%d\n",na1,nb1);
fprintf(fp1,"na|nb\t\t");
for(i=0;i<nb1;i++)
{
fprintf(fp1,"%d\t\t",i+1);
}
for(i=0;i<na1;i++)
{
fprintf(fp1,"\n %d\t",i+1);
for(j=0;j<nb1;j++)
{
fprintf(fp1,"%10.6f\t",AIC[i*nb1+j]);
}
}
fclose(fp1);
if((fp4=fopen("sita.txt","w"))==NULL)return;
fprintf(fp4,"AIC最终辨识结果sita[]为：\n");
fprintf(fp4,"na= %d\n",na1);
fprintf(fp4,"nb= %d\n",nb1);
fprintf(fp4,"\nAIC= %10.6f\n\n",AIC1);
for(i=0;i<na1;i++)
{
fprintf(fp4,"a%d=%10.6f\n",i+1,sita1[i]);
}
for(i=na1;i<na1+nb1;i++)
{
fprintf(fp4,"b%d=%10.6f\n",i+1,sita1[i]);
}
fclose(fp1);
fclose(fp2);
fclose(fp4);
sita0[0]=-1.185;sita0[1]= 0.814;sita0[2]=-0.518;sita0[3]= 0.349;sita0[4]=-0.117;
sita0[5]= 1.08 ;sita0[6]=-0.745;sita0[7]= 0.475;sita0[8]=-0.253;sita0[9]= 0.123;
CString cs;
cs.Format("模型参数：");
for(i=0;i<na1+nb1;i++)
{
if(i<na1)
cs.Format(cs+"\na%d=%.3f",i+1,sita0[i]);
else
cs.Format(cs+"\nb%d=%.3f",i-na1+1,sita0[i]);
}
m_cStaticTextREFm.SetDrawText(cs);
cs.Format("辨识参数：");
for(i=0;i<na1+nb1;i++)
{
if(i<na1)
cs.Format(cs+"\na%d=%.6f",i+1,sita1[i]);
else
cs.Format(cs+"\nb%d=%.6f",i-na1+1,sita1[i]);
}
m_cStaticTextREFi.SetDrawText(cs);
cs.Format("辨识误差：");
for(i=0;i<na1+nb1;i++)
{
if(i<na1)
cs.Format(cs+"\na%d=%.6f",i+1,fabs(sita1[i]-sita0[i]));
else
cs.Format(cs+"\nb%d=%.6f",i-na1+1,fabs(sita1[i]-sita0[i]));
}
m_cStaticTextREFe.SetDrawText(cs);
delete m_tempInt;
m_tempInt=new int[na1*nb1];
for(i=0;i<na1*nb1;i++)
{
m_tempInt[i]=int(sqrt(AIC[i])*50000);
}
m_cStaticGraph.SetDrawData(na1*nb1,m_tempInt,1,1);
m_cStaticGraph.RedrawWindow();
RedrawWindow();
free(z);
free(u);
free(e);
free(z1);
free(e1);
free(w);
free(sita);
free(sita0);
free(sita1);
free(AIC);
free(k);
free(k1);
free(k3);
free(p);
free(p0);
free(I);
}
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -