时间序列达人们来讨论一下matlab中predict函数
打开程序也看不懂，600多行代码，请达人们帮助解读一下，先谢谢各位仁兄了～～ 贴代码上来哇 你想自己编一个predict函数？ 恩，贴代码上来大家看看 回复 4# climber谢谢大家的关注～～
function =predict(varargin)
%PREDICT Computes the k-step ahead prediction.
% YP = PREDICT(MODEL,DATA,K)
%
% DATA: The output - input data as an IDDATA object, for which the
% prediction is computed.
%
% MODEL: The model as any IDMODEL object, IDPOLY, IDSS, IDARX or IDGREY.
%
% K: The prediction horizon. Old outputs up to time t-K are used to
% predict the output at time t. All relevant inputs are used.
% K = Inf gives a pure simulation of the system.(Default K=1).
% YP: The resulting predicted output as an IDDATA object. If DATA
% contains multiple experiments, so will YP.
%
% YP = PREDICT(MODEL,DATA,K,INIT)or
% YP = PREDICT(MODEL,DATA,K,'InitialState',INIT) allows the choice of
% initial state vector:
% INIT: The initialization strategy: one of
% 'e': Estimate initial state so that the norm of the
% prediction errors is minimized.
% This state is returned as X0e, see below. For multiexperiment
% DATA, X0e is a matrix whose columns contain the initial states
% for each experiment.
% 'z': Take the initial state as zero
% 'm': Use the model's internal initial state.
% 'd': Same as 'e', but if Model.InputDelay is non-zero, these delays
% are first converted to explicit model delays, so that the are
% contained in X0 for the calculation of YP.
% X0: a column vector of appropriate length to be used as initial value.
% For multiexperiment DATA, X0 may be a matrix whose columns give
% different initial states for each experiment.
% If INIT is not specified, Model.InitialState is used, so that
% 'Estimate', 'Backcast' and 'Auto' gives an estimated initial state,
% while 'Zero' gives 'z' and 'Fixed' gives 'm'. If Model is an IDARX
% model the default initial state is 'z'.
%
% With = PREDICT(MODEL,DATA,K) the initial state(s) and the
% predictor MPRED are returned.Note that if MODEL is continuous time, X0e
% is returned as the states of this model. These may differ from the initial
% states of the discrete time model that has the sampling interval(s) of
% DATA. Also when INIT = 'd' only the states of MODEL are returned in
% X0e. (To obtain the full X0, apply first INPD2NK to a discrete time
% version of the model.)
% MPRED is a cell array of IDPOLY objects,
% such that MPRED{ky} is the predictor for the ky:th output. The matching
% the channels of MPRED with data follows from its InputNames.
% See also COMPARE and IDMODEL/SIM.
% L. Ljung 10-1-89,9-9-94
% Copyright 1986-2008 The MathWorks, Inc.
% $Revision: 1.1.8.9.2.1 $$Date: 2008/08/01 15:32:39 $
% First find if there is any pair 'InitialState', init in the argument
% list:
nr = find(strncmpi(varargin,'in',2));
init =[]; xi = [];
if ~isempty(nr)
if length(varargin)<nr+1
error('Ident:general:optionsValuePair',...
'The "%s" option of the "%s" command must be specified using a name/value pair. Type "help %s" for more information.',...
'InitialState','predict','predict')
end
init = varargin{nr+1};
if ~isa(init,'char') && ~isa(init,'double')
error('Ident:analysis:IniSpec2',...
'The "InitialState" option of the "%s" command must be set to a valid character or a vector/matrix. Type "help %s" for more information.','predict','predict')
end
varargin(nr:nr+1)=[];%
end
if length(varargin)<2
disp('Usage: YP = PREDICT(MODEL,DATA)')
disp(' YP = PREDICT(MODEL,DATA,M)')
yhat = []; xi=[];
return
end
data = varargin{1};
theta = varargin{2};
if length(varargin)>2
m = varargin{3};
if ~isnumeric(m) || ~isscalar(m) || ~isreal(m) || round(m)~=m || m<1
error('Ident:analysis:predictInvalidHorizon',...
'In the "predict(MODEL,DATA,K)" command, the prediction horizon K must be a positive integer.')
end
else
m = 1;
end
if length(varargin)>3 && isempty(init)
init = varargin{4};
end
thpred = [];
yhat = [];
if isa(data,'idmodel') % Forgive order
data1 = theta;
theta = data;
data = data1;
end
nxorig = length(ssdata(theta));
= size(theta);
if isa(data,'frd') || isa(data,'idfrd')
error('Ident:analysis:predictIDFRD',...
'The "predict" command is not applicable to frequency response data (IDFRD object).')
end
if isa(data,'iddata')
= size(data);
if ny~=nym || nu~=num
error('Ident:general:modelDataDimMismatch',...
'The number of inputs and outputs of the model must match that of the data.')
end
end
if isempty(m), m=1; end
if ~isinf(m),
if m<1 || m~=floor(m)
error('Ident:analysis:predictInvalidHorizon',...
'In the "predict(MODEL,DATA,K)" command, the prediction horizon K must be a positive integer.')
end
end
% First deal with the special case of continuous time data and continuous
% model:
if isa(data,'iddata') && any(cell2mat(pvget(data,'Ts'))==0)
if ~isinf(m)
error('Ident:analysis:predictCTData',...
'The "predict" command cannot be used with continuous time data.')
end
yhat = sim(data,theta);%没懂
if nargout == 0
utidplot(theta,yhat,'Simulated')
clear yhat
end
return
end
if isempty(init)
if isa(theta,'idarx')
init='z';
else
init=pvget(theta,'InitialState');
end
init=lower(init(1));
if init=='d'
inpd= pvget(theta,'InputDelay');
if norm(inpd)==0
init = 'e';
end
end
if init=='a' || init=='b'
init = 'e';
elseif init=='f'
init = 'm';
end
elseif ischar(init)
init=lower(init(1));
if ~any(strcmp(init,{'m','z','e','d'}))
error('Ident:analysis:predictIni3',...
'The "InitialState" option of the "idmodel/predict" command must be set to ''E(stimate)'', ''Z(ero)'', ''M(odel)'', ''D(elayconvert)'', or a vector.')
end
if init=='m' && isa(theta,'idpoly')
warning('Ident:analysis:predictIni4',...
'Setting the "InitialState" option to ''m'' in the "predict" command is same as setting it to ''zero'' for IDPOLY models.')
end
end
%% Now init is either a vector or the values 'e', 'm', 'z' 'd'
if isa(data,'iddata');
if strcmp(pvget(data,'Domain'),'Frequency')
if m<inf
error('Ident:analysis:predictFreqDataFiniteK',...
'Prediction with finite horizon is not applicable to frequency domain data.')
end
end
end
if isa(data,'iddata')
= dunique(data);
iddflag = 1;
else
iddflag = 0;
uni = 1;
Tsd = 1;
inters = 'z';
end
%% Do away with continuous time right away:
Ts = pvget(theta,'Ts');
if Ts == 0
theta = pvset(theta,'CovarianceMatrix',[]); % Less calculations
ms = idss(theta);
nc = size(ms,'nx');
if ischar(init)
if init=='m';
init=pvget(ms,'X0');
end
if init=='z'
init=zeros(nc,1);
end
end
if ~isa(data,'iddata')
error('Ident:analysis:predictDataFormat',...
['For a continuous time model, data must be specified using an IDDATA object.',...
'Type "help iddata" for information on how to create the IDDATA object.'])
end
Nex = size(data,'Ne');
u = pvget(data,'InputData');
u1 = [];
for kexp = 1:length(u);
u1 = ;
end
nu = size(data,'nu');
if uni
nc = size(ms,'Nx');
= c2d(ms,Tsd,inters);
if init=='d', md = inpd2nk(md); init='e';end
% First work on init
if isa(init,'double')
=size(init);
if nxi~=nc
error('Ident:analysis:IniRows',...
'The value of the "InitialState" option should have %d row(s).',nc)
end
if mexp==1 && Nex>1
init = init*ones(1,Nex);
mexp = Nex;
end
if mexp~=Nex
error('Ident:analysis:IniSize',...
'The "InitialState" option value must have either 1 or Ne columns, where Ne = no. of data experiments.');
end
if ~isempty(Gll)
init = Gll*;
end
end
if init == 'e'
inite = x0iniest(md,data);
else
inite = init;
end
yhat = predict(md,data,m,inite); %Check XIC!
if ~isempty(Gll)
xi = inite(1:nc,:)-Gll(1:nc,nc+1:nc+nu)*u1;
else
xi = inite(1:nc,:);
end
else % Different sampling intervals
Tsd = pvget(data,'Ts');
ints = pvget(data,'InterSample');
%xic=zeros(nc,Nex);
for kexp = 1:Nex
= c2d(theta,Tsd{kexp},ints{1,kexp});
if init == 'd', md = inpd2nk(md); init = 'e';end
if isa(init,'double')
if ~isempty(Gll)
if size(init,2)==1
initk = Gll*;
else
if size(init,2)~=Nex
error('Ident:analysis:IniSize',...
'The "InitialState" option value must have either 1 or Ne columns, where Ne = no. of data experiments.');
end
initk = Gll*;
end
end
else
initk = init;
end
if init == 'e'
initk = x0iniest(md,getexp(data,kexp));
end
if ~isempty(Gll)
xi(:,kexp) = initk(1:nc)-Gll(1:nc,nc+1:nc+nu)*u1(:,kexp);
else
xi(:,kexp) = initk(1:nc);
end
yhatk = predict(md,getexp(data,kexp),m,initk);
if isempty(yhat)
yhat = yhatk;
else
yhat=merge(yhat,yhatk);
end
end
end
if nargout == 0
utidplot(md,yhat,'Predicted')
clear yhat
end
return
else %Ts>0
if init=='d'
if any(pvget(theta,'InputDelay')>0)
theta = inpd2nk(theta);
end
init = 'e';
end
end
if isinf(m)
if ischar(init) && init=='e'
X0 = x0iniest(theta,data);
init = X0;% Note the difference with COMPARE, which fits
% X0 with K = 0;
end
if ~isa(data,'iddata')
data = data(:,nym+1:end);
end
= lastwarn; was = warning('off','Ident:analysis:unstableSim');
try
yhat = sim(theta,data,init);
catch E
warning(was)
lastwarn(lw0,id0);
throw(E)
end
warning(was)
lastwarn(lw0,id0);
thpred = theta;
if nargout == 0
utidplot(theta,yhat,'Predicted')
clear yhat
elseif nargout>1
xi = init;
end
return
end
%Inpd = pvget(theta,'InputDelay');
nu = size(theta,'nu');
if ~uni
Tsd = pvget(data,'Ts'); Tsd = Tsd{1};
warning('Ident:estimation:nonUniqueDataTs',...
['The data set contains experiments with different sampling intervals.\n',...
'Sampling interval from the first experiment (=%g) will be used.'],Tsd)
end
if iddflag && any(abs(Ts-Tsd)>1e4*eps)
warning('Ident:analysis:dataModelTsMismatch',...
'The data and model sampling intervals are different. The data sampling interval will be used.')
end
Inpd = pvget(theta,'InputDelay');
if init=='d'
init = 'e';
theta = inpd2nk(theta);
Inpd = zeros(size(Inpd));
end
if isa(data,'iddata')
if isnan(data)
error('Ident:general:idprepMissingData',...
'Data contains NaNs which represent missing values. Use the "misdata" command to fill in the missing values before using.');
end
data = nkshift(data,Inpd,'append');
theta = pvset(theta,'InputDelay',zeros(nu,1));
= idprep(data,0,'dummy');
%if ~isempty(errflag.message), error(errflag), end
if ~isempty(Name), data.Name = Name; end
iddatflag = 1;
else
if norm(Inpd)>eps
if iscell(data)
data = data{1};
end
= size(data);
nk1 = Inpd;
ny = nudum - length(nk1);
Ncc = min();
for ku = 1:length(nk1)
u1 = data(max()-nk1(ku)+1:Ncc-nk1(ku),ny+ku);
newsamp = Ncap-length(u1);
if nk1(ku)>0
u1= ;
else
u1 = ;
end
data(:,ny+ku) = u1;
end
end
iddatflag= 0;
if ~iscell(data), ze ={data};else ze = data;end
Ne = length(ze);
nz = size(ze{1},2);
for kexp = 1:Ne
Ncaps(kexp) = size(ze{kexp},1);
end
%Ts =[];
end
if m>=min(Ncaps)
error('Ident:analysis:predictLargeK',...
'Prediction horizon must be smaller than the number of data samples.')
end
=ssdata(theta);
=size(C);nut=size(B,2);
if iddatflag
if ny~=nyt || nu~=nut
error('Ident:general:modelDataDimMismatch',...
'The number of inputs and outputs of the model must match that of the data.')
end
else
if nz~=nyt+nut
error('Ident:general:modelDataDimMismatch',...
'The number of inputs and outputs of the model must match that of the data.')
end
ny = nyt; nu = nut;
end
if strcmp(init,'e')
X0 = x0iniest(theta,ze);
elseif strcmp(init,'z')
X0 = zeros(size(X0));
end
if ~ischar(init)
X0 = init;
end
if nargout>1
xi = X0;
end
= size(X0);
if xnc~=1 && xnc~=Ne
error('Ident:analysis:IniSize',...
'The "InitialState" option value must have either 1 or Ne columns, where Ne = no. of data experiments.');
end
if xnr~=nx
error('Ident:analysis:IniRows',...
'The value of the "InitialState" option should have %d row(s).',nx)
end
if xnc==1 && Ne>1
X0= X0*ones(1,Ne);
end
for kexp = 1:Ne
z = ze{kexp};
u = z(:,1+ny:end);
Ncap = Ncaps(kexp);
if m==inf,
yhat=sim(theta,u,X0(:,kexp));
else
x=ltitr(A-K*C,, z, X0(:,kexp));
if m==1,
yhat=(C*x.').';
if ~isempty(D),yhat=yhat + (D*u.').';end
else
F=D;Mm=eye(length(A));
for km=1:m-1
F=;
Mm=A*Mm;
end
yhat=zeros(Ncap,ny);%corr 911111
for ky=1:ny
for ku=1:nu
yhat(:,ky)=yhat(:,ky)+filter(F(ky,ku:nu:m*nu),1,u(:,ku));
end
end
if isempty(yhat),yhat=zeros(Ncap,ny);end
yhat(m:Ncap,:)=yhat(m:Ncap,:)+(C*Mm*x(1:Ncap-m+1,:).').';
if nu>0
x=ltitr(A,B,u(1:m,:),X0(:,kexp));
yhat(1:m,:)=(C*x.').';
end
if ~isempty(D),yhat(1:m,:)=yhat(1:m,:) + (D*u(1:m,:).').';end
end
end
yhatc{kexp} = yhat;
end
if iddatflag
yhat = data;
yhat = pvset(yhat,'OutputData',yhatc,'InputData',[]);
else
yhat = yhatc;
end
if nargout >2
thpred = polypred(theta,m);
end
if nargout > 1
try
xi = xi(1:nxorig,:);
end
end
if (nargout == 0)
% Plot y and yhat.
utidplot(theta,yhat,'Predicted');
clear yhat x0;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% LOCAL FUNCTIONS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function xi = x0iniest(m,data)
if isa(data,'iddata') && strcmpi(pvget(data,'Domain'),'frequency')
= pe_f(m,data);
return
end
Inpd = pvget(m,'InputDelay');
if isa(data,'iddata')
data = nkshift(data,Inpd,'append');
= idprep(data,0,'dummy');
%if ~isempty(errflag.message), error(errflag), end
if ~isempty(Name)
data.Name = Name;
end
else
if iscell(data)
ze = data;
else
ze = {data};
end
end
alg=pvget(m,'Algorithm');
maxsize=alg.MaxSize;
=ssdata(m);
nx=length(A);
if ischar(maxsize)
maxsize = idmsize(length(ze{1}),nx);
end
AKC=A-K*C;
=size(C);
nu=size(B,2);
el=zeros(0,ny);
xic = zeros(nx,0);
Ne = length(ze);
for kexp = 1:Ne
z = ze{kexp};
=size(z);
if nu+ny~=nz
error('Ident:general:modelDataDimMismatch',...
'The number of inputs and outputs of the model must match that of the data.')
end
nz=ny+nu; = size(z); n = nx;
rowmax = nx+nz; X0 = zeros(nx,1);
M=floor(maxsize/rowmax);
if ny>1 || M<Ncap
R=zeros(n,n);Fcap=zeros(n,1);R1=[];
for kc=1:M:Ncap
jj=(kc:min(Ncap,kc-1+M));
if jj(length(jj))<Ncap
jjz = ;
else
jjz=jj;
end
psitemp=zeros(length(jj),ny);
psi=zeros(ny*length(jj),n);
x=ltitr(AKC,,z(jjz,:),X0);
yh=(C*x(1:length(jj),:).').';
nm=pvget(m,'NoiseVariance');
if isempty(nm) || norm(nm)==0 % To handle models without noise
nm = eye(ny);
end
sqrlam=pinv(sqrtm(nm));
if ~isempty(D),yh=yh+(D*z(jj,ny+1:ny+nu).').';end
e=(z(jj,1:ny)-yh)*sqrlam;
=size(x);X0=x(nxr,:).';
evec=e(:);
kl=1;
for kx=1:nx
if kc==1
x0dum=zeros(nx,1);x0dum(kx,1)=1;
else
x0dum=X00(:,kl);
end
psix=ltitr(AKC,zeros(nx,1),zeros(length(jjz),1),x0dum);
=size(psix);
X00(:,kl)=psix(rp,:).';
psitemp=(C*psix(1:length(jj),:).').'*sqrlam;
psi(:,kl)=psitemp(:);kl=kl+1;
end
if ~isempty(R1)
if size(R1,1)<n+1
error('Ident:estimation:predictSmallMaxSize',...
'The value of the algorithm property "MaxSize" is too small to estimate the initial states. Either increase MaxSize value in the model or set the "InitialState" option of the "predict" command to a value other than ''Estimate''.')
end
R1=R1(1:n+1,:);
end
H1 = ];R1 = triu(qr(H1));
end
try
xi(:,kexp) = pinv(R1(1:n,1:n))*R1(1:n,n+1);
catch
warning('Ident:estimation:X0EstFailed',...
'Failed to estimate initial conditions. Check model stability. \nThe "InitialState" option has been set to zero.')
xi(:,kexp) = zeros(n,1);
end
else
%% First estimate new value of xi
x=ltitr(AKC,,z);
y0=x*C';
if ~isempty(D),
y0=y0+(D*z(:,ny+1:ny+nu).').';
end
psix0=ltitr(AKC.',C.',);
psix0=psix0(2:end,:);
try
xi(:,kexp) = pinv(psix0)*(z(:,1)-y0);
catch
warning('Ident:estimation:X0EstFailed',...
'Failed to estimate initial conditions. Check model stability. \nThe "InitialState" option has been set to zero.')
xi(:,kexp) = zeros(n,1);
end
end
end
%%*************************************************************************
function thpred = polypred(theta,m)
% Note that cross couplings between output
% channels are ignored in these calculations.
= size(theta);
yna = pvget(theta,'OutputName');
una = pvget(theta,'InputName');
yu = pvget(theta,'OutputUnit');
uu = pvget(theta,'InputUnit');
for ky = 1:ny
= polydata(theta(ky,:));
if nu>0
ff = 1;
for ku = 1:nu,
bt = b(ku,:);
for kku = 1:nu, if kku~=ku, bt = conv(bt,f(kku,:)); end, end
bb(ku,:) = bt;
ff = conv(ff,f(ku,:));
end
a = conv(conv(a,ff),d);c=conv(c,ff);
else
a = conv(a,d);
end
na = length(a); nc = length(c); nn = max(na,nc);
a = ; c = ;
= deconv(conv(,c),a);
ng=length(g);
if nu>0,
df=conv(d,f);
for ku=1:nu
bf(ku,:)=conv(bb(ku,:),df);
end
nbf=length(bf(1,:));nn=max(ng,nbf);
gg=[;];
else
gg=g;
end
th1 = idpoly(c,gg);
th1 = pvset(th1,'InputName',,'OutputName',,...
'InputUnit',,'OutputUnit',yu(ky),'Ts',pvget(theta,'Ts'),...
'TimeUnit',pvget(theta,'TimeUnit'),'EstimationInfo',pvget(theta,'EstimationInfo'));
thpred{ky} = th1;
end
%%%************************************************************************
页:
[1]
