Global Index (short | long) | Local contents | Local Index (short | long)
[th,it_inf]=pem(z,nn,index,maxiter,tol,lim,maxsize,Tsamp,trc)
PEM Computes the prediction error estimate of a general linear model.
TH = PEM(Z,THSTRUC) or TH = PEM(Z,THSTRUC,'trace')
TH : returned as the estimated parameters of the model
along with estimated covariances and structure information.
For the exact format of TH see also THETA.
Z : the output input data [y u], with y and u as column vectors
For multi-variable systems, Z=[y1 y2 ... yp u1 u2 ... un]
THSTRUC: A matrix that defines the model structure. Typically created
by (see also) POLY2TH, MS2TH or MF2TH or by some estimation routine.
The minimization is then initialized at the parameters given in THSTRUC.
A general, black-box multi-input structure
A(q) y(t) = [B(q)/F(q)] u(t-nk) + [C(q)/D(q)] e(t)
is obtained for THSTRUC = [na nb nc nd nf nk], indicating the orders
in the above model.
By TH=PEM(Z,THSTRUC,INDEX) only parameters corresponding to the indices
in the row vector INDEX are estimated.
If a *last* argument to pem is entered as 'trace', - like
TH=PEM(Z,THSTRUC,INDEX,'trace') - , information about the criterion
minimization is given in the MATLAB Command Window.
Some parameters associated with the algorithm are accessed by
[TH,IT_INF] = PEM(Z,THSTRUC,index,maxiter,tol,lim,maxsize,T)
See also AUXVAR for an explanation of these and their default values.
See also ARX, ARMAX, BJ, IV4, N4SID, OE.
| This function calls | This function is called by |
|---|---|
function [th,it_inf]=pem(z,nn,index,maxiter,tol,lim,maxsize,Tsamp,trc)
% L. Ljung 10-1-86, 7-25-94
% Copyright (c) 1986-98 by The MathWorks, Inc.
% $Revision: 2.6 $ $Date: 1997/12/03 23:16:04 $
if nargin<2
disp('Usage: TH = PEM(Z,NN)')
disp(' TH = PEM(Z,NN,''trace'')')
disp(' TH = PEM(Z,NN,index,maxiter,tol,lim,maxsize,T,''trace'')')
return
end
% *** Set up default values ***
global XIDiter
gui=0;display=0;stopp=0;flag=[];
eval('flag=findobj(get(0,''children''),''flat'',''tag'',''sitb22'');','');
if ~isempty(flag),if strcmp(get(flag,'vis'),'on')
gui=1;set(XIDiter(7),'userdata',0);
end,end
[Ncap,nz]=size(z);
maxsdef=idmsize(Ncap);
if nargin==9,if strcmp(lower(trc),'trace'),display=1;end,end
if nargin==8,if strcmp(lower(Tsamp),'trace'),display=1;Tsamp=[];end,end
if nargin==7,if strcmp(lower(maxsize),'trace'),display=1;maxsize=[];end,end
if nargin==6,if strcmp(lower(lim),'trace'),display=1;lim=[];end,end
if nargin==5,if strcmp(lower(tol),'trace'),display=1;tol=[];end,end
if nargin==4,if strcmp(lower(maxiter),'trace'),display=1;maxiter=[];end,end
if nargin==3,if strcmp(lower(index),'trace'),display=1;index=[];end,end
if nargin<8, Tsamp=[];end
if nargin<7, maxsize=[];end
if nargin<6, lim=[];end
if nargin<5, tol=[];end
if nargin<4, maxiter=[];end
if nargin<3, index=[];end
[nr,cc]=size(nn);
if isempty(Tsamp),if nr>1,Tsamp=nn(1,2);else Tsamp=1;end,end
if isempty(maxsize),maxsize=maxsdef;end,if maxsize<0,maxsize=maxsdef;end
if isempty(lim),lim=1.6;end,if lim<0,lim=1.6;end
if isempty(tol),tol=0.01;end,if tol<0,tol=0.01;end
if isempty(maxiter),maxiter=10;end,if maxiter<0,maxiter=10;end
if nr>1,if isthss(nn)
eval(['[th,it_inf]=pemss(z,nn,index,maxiter,tol,',...
'lim,maxsize,Tsamp,display,gui);'])
return,end,end
if Tsamp<0, Tsamp=1;end
% *** Do some consistency tests ***
nu=nz-1;
if nz>Ncap, error('The data should be organized in column vectors.')
return,end
nb=nn(2:nz);
if nz==1
error('For a scalar time series, use ARMAX instead.')
end
if nr==1
if cc~=3+3*nu,
disp('Incorrect number of orders specified.')
disp('You should have nn=[na nb nc nd nf nk],')
disp('where nb, nf and nk are row vectors of length equal to')
disp(' the number of inputs'),error(' ')
end
if any(nn<0)
error('All orders must be non-negative.')
end
kz=find(nb(1:nu)==0);
if ~isempty(kz),nn(nu+kz+3)=zeros(1,length(kz));end % Let the corre-
% sponding nf also be zero
if sum(nb)==0&nu>0
if sum([nn(1:3+nu)])==0 % This would be an OE structure without B
error('For this structure it does not make sense to have nb all zero.')
else
[th,it_inf]=armax(z(:,1),[nn(1)+nn(3+nu) nn(2+nu)],maxiter,tol,...
lim,maxsize,Tsamp,display,gui);
return
end
end
end
%*** if nn has more than 1 row, then it is a theta matrix
% and we jump directly to mincrit ***
%
if nr>1
th=nn;
if th(1,3)~=nu
error(['The number of inputs in the data set is not',...
'consistent with the chosen model structure.']);
end
end
if nr==1
th=inival(z,nn,maxsize,Tsamp);
if isempty(th)
disp('WARNING: Your input signal is apparently not persistently exciting.')
disp(' Use other input or lower model orders.')
disp(' No model returned.')
th=[];return
end
end
% *** Display the initial estimate ***
if gui,drawnow,display=get(XIDiter(2),'value');end
if display
%clc
disp([' INITIAL ESTIMATE'])
disp(['Current loss:' num2str(th(1,1))])
disp(['th-vector'])
nu=th(1,3); n=sum(th(1,4:6+2*nu));
disp(th(3,1:n)')
end
if maxiter==0,return,end
% *** Minimize the prediction error criterion ***
[th,it_inf]=mincrit(z,th,index,maxiter,tol,lim,maxsize,Tsamp,display,gui);