endless.pas

对于平稳序列,运用ar模型计算其参数

const
datanum =20;           {数据个数}
paranum =2;            {模型参数}
type
a1=array[0..datanum] of real;
a2=array[0..datanum,0..paranum] of real;
a3=array[0..paranum,0..datanum] of real;
a4=array[0..paranum,0..paranum+1] of real;
a5=array[0..paranum] of real;
a7=array[0..paranum+1] of real;
a6=array[0..paranum] of integer;
var
y : a1;                {y为数组}
q : a2;                {q为测量矩阵}
p : a3;                {p为q的转置}
h : a4;                {h为q.T*q}
b,X : a5;              {b为q.T*y；x为待求AR模型参数}
w : a6;
xx :  a7;
filename1,filename2,filename3 ,filename4: text;
aic,aico : real;

Procedure Ein(var m : a1;var q : a2);           {读入数据}
var i,j,k : integer;
begin
  assign(filename1,'f:\shiyan\data.txt');
  reset(filename1);
  for i:=1 to datanum do
    begin
      readln(filename1,y[i]);
      writeln(y[i]);
    end;
  writeln('display q:');
  for j:=1 to datanum-paranum do
  begin
    for k:=1 to paranum do
    begin
      q[j,k]:=y[j+paranum-k];
      write(q[j,k]:12:5);
    end;
    writeln;
  end;
  writeln('display q.t: ' );
end;

Procedure  change(var q : a2);                    {求q的转置}
var i,j,k : integer;
begin
for k:=1 to paranum do
begin
  for i:=1 to datanum-paranum do
    begin
    p[k,i]:=q[i,k];
    write(p[k,i]:10:5,'      ');
    end;
    writeln;
    end;
writeln('display  q.T*q:  ');
end;

procedure mutipl1;                               {求q.T*q}
var m,n,l : integer;
begin
for m:=1 to paranum do
 begin
   for n:=1 to paranum do
    begin
     h[m,n]:=0;
     for l:=1 to datanum-paranum do
     h[m,n]:=h[m,n]+p[m,l]*q[l,n];
     write(h[m,n]:10:5,'      ');
    end;
    writeln;
 end;
 writeln('display q.T*y:');
end;


procedure mutipl2;                          {求q.T*y}
var i,j : integer;
begin
  for i:=1 to paranum do
  begin
    b[i]:=0;
    for j:=1 to datanum-paranum do
    b[i]:=b[i]+p[i,j]*y[j+paranum];
    writeln(b[i]:10:5);
  end;
  writeln('X:');
end;

procedure gause;                              {求解模型估计参数X}
var i,j,k,l,n : integer;
     r,e : real;
begin
for i:=1 to paranum do h[i,paranum+1]:=b[i]; {化为增广矩阵}
for i:=1 to paranum do  {从每行开始做如下循环作业(i循环)}
   begin
      l:=i;   n:=0;  e:=h[i,0];  {p,q,e赋值}
      for j:=i to paranum do    {从节点行以下所有行作如下动作(j循环)}
         for k:=0 to paranum do{k循环}
            if abs(h[j,k])>abs(e) then
            begin e:=h[j,k];
            n:=k;
            l:=j;
            end; {从第j行中求出最大值，配合j循环可得第i行后所有元素中的最大值，并确定其位置p,q}
      if abs(e)>1.0E-10 then{if1}{如此行后最大值小于10-10  ，视此行后所有元素为0，不再计算}{此为if1}
      begin
        if l<>i then {此为if2}{如p=I,即最大值在此行，不调整换行，否则换行}
        begin
        for k:=1 to paranum+1 do begin xx[k]:=H[i,k];H[i,k]:=H[l,k];H[l,k]:=xx[k];end; {换行，i行与此行后最大值行互换}
        end; {if2}
        for j:=1 to paranum do   {j循环}
          begin if (j<>i) and (H[j,n]<>0) then{消去除最大值所在行的其余行最大值所在列元素}
            begin r:=H[j,n]/H[i,n]; for k:=1 to paranum+1 do H[j,k]:=H[j,k]-H[i,k]*r;end;
          end;
        w[i]:=n; {每行非零系数列下标}
      end; {if1}
   end; {(i循环)结束}
   readln;
   for i:=1 to paranum do
     begin
     n:=w[i];
     x[n]:=H[i,paranum+1]/h[i,n];
     end; {求值}

   for i:=1 to paranum do writeln(x[i]:12:5);
   readln;
end;

procedure  jianyan;
var i,j,k : integer;
    temp,demp,sigma2,sigma,aico : real;
    del : a1;
begin
for i:=1 to datanum-paranum do
  begin
    temp:=0;
   for j:=1 to paranum do
    begin
    temp:=temp+y[j+i-1]*x[paranum-j+1];
    end;
  del[i]:=y[paranum+i]-temp;
  writeln(del[i]);
end;
readln;
demp:=0;
for k:=1 to datanum-paranum do
    begin
    demp:=demp+del[k]*del[k];
    end;
  sigma2:=demp/(datanum-paranum);       {sigma2为均方差}
  writeln('sigma2=');
  writeln(sigma2);
  readln;

  sigma:=sqrt(sigma2);
  assign(filename3,'f:\shiyan\cancha.txt');
  append(filename3);
  writeln(filename3,sigma);
  close(filename3);

  assign(filename2,'f:\shiyan\canshu.txt');
  append(filename2);
  writeln(filename2,x[1]/sigma,'       ',x[2]/sigma);
  close(filename2);

  aico:=datanum*LN(sigma2)+2*paranum;
  assign(filename4,'f:\shiyan\aico.txt');
  append(filename4);
  writeln(filename4,aico);
  close(filename4);
  writeln('aico=');
  writeln(aico);
  readln;
end;


begin
                                 {主程序}
  Ein(y,q);
  change(q);
  mutipl1;
  mutipl2;
  gause;
  jianyan;
end.