msls6.c

在系统辨识过程中的多步最小二乘方法,比一步最小二乘精度更高

#include "stdio.h"
#include "stdlib.h"
#include "math.h"
#include "brmul.c"
#include "yrinv.c"
int main()
{ 
	FILE *fp1,*fp2,*fp3,*fp4;
	static double h[51][1],u[651],e[651],z[651],z1[601][1],y[651],y1[600][1],v[651],v1[651],pp[5][5],ss[5][1];
	static double u1[601][51],u2[51][601],w[5][1],w1[1][5],s[5][1],s1[5][1],c[2][1],o[1][2],o1[2][1],p[5][5];
	static double q[51][51],qu[51][601],w1p[1][5],pw[5][1],k[5][1],g[2][2],c1[2][1],gg[2][2];
	static double a,b,wpw[1],w1s[1],k1,err,ogo[1],o1c[1],o1g[1][2],go[2][1],k2[2][2],b1;
	int i,j,n,m;
	fp1=fopen("h1.txt","w");
	fp2=fopen("M.txt","r");
	fp3=fopen("wnoise1.txt","r");
	fp4=fopen("msls6.txt","w");

	for(i=0;i<651;i++)
	fscanf(fp2,"%lf",&u[i]);
	for(i=0;i<651;i++)
	fscanf(fp3,"%lf",&e[i]);


	v[0]=e[0];
	v[1]=-1.0*v[0]+e[1];
	for(i=2;i<651;i++)
	v[i]=-1.0*v[i-1]-0.41*v[i-2]+e[i];

	z[0]=v[0];
	z[1]=-0.9*z[0]+0.7*u[0]+v[1];
	z[2]=-0.9*z[1]-0.15*z[0]+0.7*u[1]-1.5*u[0]+v[2];
	for(i=3;i<651;i++)
		z[i]=-0.9*z[i-1]-0.15*z[i-2]-0.02*z[i-3]+0.7*u[i-1]-1.5*u[i-2]+v[i];
	for(i=0;i<601;i++)
		z1[i][0]=z[i+50-1];

	w1s[0]=0.0;
	wpw[0]=0.0;
	for(i=0;i<5;i++)
		p[i][i]=1.0e+6;

	for(i=0;i<=600;i++)
		for(j=0;j<=50;j++)
			u1[i][j]=u[50-j+i];
	for(i=0;i<=50;i++)
		for(j=0;j<=600;j++)
			u2[i][j]=u1[j][i];

	brmul(u2,u1,51,601,51,q);
	yrinv(q,51);
	brmul(q,u2,51,51,601,qu);
	brmul(qu,z1,51,601,1,h);
	for(i=0;i<51;i++)
		fprintf(fp1,"%lf\n",h[i][0]);
	fclose(fp1);
	fclose(fp2);
	fclose(fp3);

	for(i=0;i<651;i++)
	{
		a=0.0;
		b=0.0;
		if(i<50)
		{
			for(j=0;j<=i;j++)
				a=a+h[j][0]*u[i-j];
			y[i]=a;
		}
		else
		{	for(j=0;j<=50;j++)
				b=b+h[j][0]*u[i-j];
			y[i]=b;
		}
	}

	w[0][0]=-z[0];
	w[3][0]=u[0];
	n=0;

	for(m=0;m<600;m++)
	{ 
		for(i=0;i<5;i++)
			s[i][0]=s1[i][0];
		for(i=0;i<5;i++)
		w1[0][i]=w[i][0];

		brmul(w1,p,1,5,5,w1p);
		brmul(w1p,w,1,5,1,wpw);
		k1=1.0/(wpw[0]+1.0);
		brmul(p,w,5,5,1,pw);
		for(i=0;i<5;i++)
			k[i][0]=pw[i][0]*k1;
		brmul(w1,s,1,5,1,w1s);
		b=z[n+1]-w1s[0];
		for(i=0;i<5;i++)
			s1[i][0]=s[i][0]+k[i][0]*b;
		brmul(pw,w1p,5,1,5,pp);
		for(i=0;i<5;i++)
			for(j=0;j<5;j++)
			p[i][j]=p[i][j]-pp[i][j]/(1.0+wpw[0]);
		n=n+1;
		w[0][0]=-z[n];
		w[1][0]=-z[n-1];
		w[2][0]=-z[n-2];
		w[3][0]=u[n];
		w[4][0]=u[n-1];
	}
	for(i=0;i<5;i++)
	{
		printf("%lf\n",s1[i][0]);
		fprintf(fp4,"%lf  ",s1[i][0]);
	}
	fprintf(fp4,"\n");
	ss[0][0]=0.9;
	ss[1][0]=0.15;
	ss[2][0]=0.02;
	ss[3][0]=0.7;
	ss[4][0]=-1.5;
	err=0.0;
	for(i=0;i<5;i++)
		err=err+(s1[i][0]-ss[i][0])*(s1[i][0]-ss[i][0]);
	printf("误差平方和:%lf\n",err);
	o1c[0]=0.0;
	ogo[0]=0.0;
	n=0;
	for(i=0;i<2;i++)
		g[i][i]=1.0e+6;
	v1[0]=z[0]+0.0*u[0];
	v1[1]=z[1]+s1[0][0]*z[0]-0.0*u[1]-s1[3][0]*u[0];
	v1[2]=z[2]+s1[0][0]*z[1]+s1[1][0]*z[0]-0.0*u[2]-s1[3][0]*u[1]-s1[4][0]*u[0];
	for(i=3;i<651;i++)
		v1[i]=z[i]+s1[0][0]*z[i-1]+s1[1][0]*z[i-2]+s1[2][0]*z[i-3]-0.0*u[i]-s1[3][0]*u[i-1]-s1[4][0]*u[0];
	for(i=0;i<651;i++)
		v1[i]=v1[i];
	o[0][0]=v1[0];
	o[0][1]=0;
	for(m=0;m<600;m++)
	{ 
		for(i=0;i<2;i++)
			c[i][0]=c1[i][0];
		for(i=0;i<2;i++)
			o1[i][0]=o[0][i];
		brmul(o1,g,1,2,2,o1g);
		brmul(o1g,o,1,2,1,ogo);
		k1=1.0/(ogo[0]+1.0);
		brmul(g,o,2,2,1,go);
		for(i=0;i<2;i++)
			k2[i][0]=go[i][0]*k1;
		brmul(o1,c,1,2,1,o1c);
		b1=v[n+1]-o1c[0];
		for(i=0;i<2;i++)
			c1[i][0]=c[i][0]+k2[i][0]*b1;
		brmul(go,o1g,2,1,2,gg);
		for(i=0;i<2;i++)
			for(j=0;j<2;j++)
				g[i][j]=g[i][j]-gg[i][j]/(1.0+ogo[0]);
		n=n+1;
		o[0][0]=-v[n];
		o[0][1]=-v[n-1];
	}
	for(i=0;i<2;i++)
	{
		printf("%lf\n",c1[i][0]);
		fprintf(fp4,"%lf  ",c1[i][0]);
	}
	fclose(fp4);
	err=0.0;
	err=(c1[0][0]-1.0)*(c1[0][0]-1.0)+(c1[1][0]-0.41)*(c1[1][0]-0.41);
	printf("误差平方和为：%lf\n",err);
	return 0;
}