当前位置：文档之家 › AIC法定阶的依阶次递推算法程序.

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

  • 格式：doc
  • 大小：39.50 KB
  • 文档页数：6

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 <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];doubleq24[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();}。

合集下载

应用最小二乘一次完成法和递推最小二乘法算法的系统辨识讲解

应用最小二乘一次完成法和递推最小二乘法算法的系统辨识讲解

1最小二乘法的理论基础1.1最小二乘法设单输入单输出线性定长系统的差分方程表示为：其中δ(k)为服从N(0,1)的随机噪声,现分别测出n+N 个输出输入值y(1),y(2),…,y(n+N),u(1),u(2),…,u(n+N),则可写出N 个方程，写成向量－矩阵形式(4.1.1)()()()()()()()()1201121n n y k a y k a y k a y k n b u k b u k b u k n k ξ=-------++-++-+()()()()()()101122,,n n a y n n y n a n y b y n N n N b ξξθξξ⎡⎤⎢⎥++⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥++⎢⎥⎢⎥===⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥++⎢⎥⎢⎥⎣⎦⎣⎦⎢⎥⎢⎥⎣⎦()()()()()()()()()()()()()()()()()()10111121222112n n y n y n y u n u y n y n y u n u y n N y n N y N u n N u N a n a n b n N b ξξξ+--+⎡⎤⎡⎤⎢⎥⎢⎥+-+-+⎢⎥⎢⎥=⨯⎢⎥⎢⎥⎢⎥⎢⎥+-+--+⎢⎥⎢⎥⎣⎦⎣⎦⎡⎤⎢⎥+⎡⎤⎢⎥⎢⎥⎢⎥+⎢⎥+⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥+⎢⎥⎣⎦⎢⎥⎢⎥⎣⎦则式(1.1.1)可写为 （4.1.2） 式中：y 为N 维输出向量；ξ为N 为维噪声向量；θ为（2n+1）维参数向量；Φ为N ×（2n+1）测量矩阵。

因此，式（4.1.1）是一个含有（2n+1）个未知参数，由N 个方程组成的联立方程组。

11y θφφξ--=-在给定输出向量y 和测量矩阵Φ的条件下求参数θ的估计，这就是系统辨识问题。

设 表示 θ 的估计值，ŷ表示y 的最优估计，则有 （4.1.3） 式中：()()()10ˆˆ1ˆˆ2ˆˆ,ˆˆˆn n ay n a y n y b y n N b θ⎡⎤⎢⎥+⎡⎤⎢⎥⎢⎥⎢⎥+⎢⎥⎢⎥==⎢⎥⎢⎥⎢⎥⎢⎥+⎢⎥⎣⎦⎢⎥⎢⎥⎣⎦设e(k)=y(k)- ŷ(k), e(k)称为残差，则有e=y- ŷ=y-Φθ 最小二乘估计要求残差的平方和最小，即按照指数函数（4.1.4）求Ｊ对 的偏导数并令其等于0可得：（4.1.5）由式（4.1.5）可得的 θ 最小二乘估计：（4.1.6）J 为极小值的充分条件是：即矩阵ΦT Φ为正定矩阵，或者说是非奇异的。

系统辨识与参数估计_Fixed

系统辨识与参数估计_Fixed

第一题 递推最小二乘估计参数考虑如下图所示的仿真对象，图中)(k v 是服从)1,0(N 正态分布的不相关随机噪声；输入信号)(k u 采用4阶逆M 序列，特征多项式取41)(s s s F ⊕⊕=，幅度为1，循环周期为bit N p 62=；控制λ值，使数据的噪信比为73％。

加权因子1)(=Λk ；数据长度L=500;初始条件取001.0)0(ˆ,10)0(6==θI P ， （1） 利用递推最小二乘算法在线估计参数，（2） 利用模型阶次辨识方法（AIC 准则），确定模型的阶次。

（3） 估计噪声)(k v 的方差和模型静态增益K （4） 作出参数估计值随时间的变化图 答：设过程的输入输出关系可以描述成()()()T z k h k n k θ=+()z k 是输出量，()h k 是可观测的数据向量，n (k)是均值为0的随机噪声 []()(1),(2),(1),(2)Th k z k z k u k u k =------[]1212,,,Ta ab b θ=选取的模型为结构是1212()(1)(2)(1)(2)z k a z k a z k bu k bu k =----+-+- 12121.5,0.7, 1.0,0.5a a b b =-===加权最小二乘参数估计递推算法RWLS 的公式如下，11()(1)()()(1)()()()(1)()()()(1)()()()(1)T T TK k p k h k h k p k h k k k k K k z k h k k p k I K k h k p k θθθ-⎡⎤=--+⎢⎥Λ⎣⎦⎡⎤=-+--⎣⎦⎡⎤=--⎣⎦为了把p(k)的对称性，可以把p(k)写成1()(1)()()()(1)()()T T p k p k K k K k h k p k h k k ⎡⎤=---+⎢⎥Λ⎣⎦ 如果把()k Λ设成1的时候，加权最小二乘法就退化成最小二乘法。

组合数学(4)递推递归母函数

组合数学(4)递推递归母函数

ACM 暑期集训 组合数学（4） 递推 递归 母函数1 递推关系序列{a n }=a 0，a 1，…，a n ，…，把 a n 与某些a i (i <n )联系起来的等式叫做关于序列{a n }的递推方程。

当给定递推方程和适当的初值就唯一确定了序列。

递推关系分类： （1）按常量部分：齐次递推关系：指常量＝0，如F(n)=F(n-1)+F(n-2) 非齐次递推关系：指常量≠0，如F(n)=2*F(n-1)+1 （2）按运算关系：线性关系，如上面的两个；非线性关系，如F(n)=F(n-1)*F(n-2)。

（3）按系数：常系数递推关系，如（1）中的两个；变系数递推关系，如D(n)=(n-1)(D(n-1)+D(n-2)。

（4）按数列的多少一元递推关系，只涉及一个数列，上面的均为一元； 多元递推关系，涉及多个数列，如⎩⎨⎧+=+=----111177n n nn n n a b b b a a Fibonacci 数列为1,1,2,3,5,8,13,.....long long data[100]; data[1]=1; data[2]=1;for(int i=3;i<=50;i++) data[i]=data[i-1]+data[i-2]; while(cin>>n) cout<<data[n]<<endl;例1：直线割平面问题。

在一个无限的平面上有N 条直线，试问这些直线最多能将平面分割成多少区域？F(1) = 2; F(2) = 4; F(3) = 7; F(n)=F(n-1)+n; (n>1)int recurrence(int n) //递推 {f[1]=2;for(i=2;i<=n;i++) f[n]=f[n-1]+n; return f[n]; }int recursion(int n) 递归： {if(n==1) return 2;//递归终止条件 else return recursion(n-1)+n; }更快的方法是求出通项：F(n)=n^(n+1)/2+1例2：HDOJ2050 折线割平面问题在一个无限的平面上有N 条折线，试问这些折线最多能将平面分割成多少区域？F(n)=F(n-1)+4n-3; F(n)=2*n^2-n+1;例3：椭圆割平面问题。

递推算法的一般步骤

递推算法的一般步骤

递推算法的一般步骤递推算法是一种常用的计算方法，通过已知条件推导出未知结果的一种数学推理过程。

它可以用于解决一些复杂的问题，特别是涉及到递归关系的情况。

递推算法的一般步骤如下：1. 首先，明确问题的递推关系。

递推算法的核心就是通过已知条件推导出未知结果，因此需要通过观察问题的特点，找到问题的递推关系。

这个递推关系可以是一个公式、一个递归函数或者一系列等差、等比数列等。

2. 然后，确定初始条件。

递推算法是从已知条件推导出未知结果的过程，因此需要确定问题的初始条件。

这个初始条件可以是问题中给出的初始值，也可以是手动给定的一些值。

3. 接下来，进行递推计算。

根据问题的递推关系和初始条件，可以开始逐步计算得到未知结果。

递推计算的过程就是反复应用递推关系，将已知的结果代入递推关系中得到新的结果，然后再将新的结果代入递推关系中得到更新的结果，依此类推，直到得到想要的未知结果。

4. 最后，进行结果验证。

得到未知结果后，需要对结果进行验证，确保结果的正确性。

这可以通过与已知条件进行对比、与问题的实际情况进行比较等方法来进行。

如果结果符合要求，即可确认递推算法的正确性。

除了以上步骤，还有一些在使用递推算法时需要注意的问题：1. 递归深度的控制。

递推算法中经常使用递归函数来表达递推关系，但要注意递归的深度，避免出现无限递归的情况。

可以通过设置递归深度的限制条件或者使用迭代方法来避免这种情况的发生。

2. 边界条件的处理。

在确定初始条件时，要考虑问题的边界条件。

边界条件是递推算法中的特殊情况，需要单独处理，以确保算法的正确性。

边界条件可以是问题中给出的某种特殊情况，或者是根据问题的特点手动给定的一些条件。

3. 问题的优化。

在实际应用中，递推算法可能会面临一些效率问题，特别是在处理大规模数据时。

因此，需要对算法进行优化，提高算法的执行效率。

可以通过改进递推关系、使用动态规划等方法来实现问题的优化。

综上所述，递推算法是一种通过已知条件推导出未知结果的计算方法。