Mercurial > forge
changeset 984:3451c5f4a428 octave-forge
aincluded aarmam and update INDEX
| author | schloegl |
|---|---|
| date | Tue, 24 Jun 2003 19:27:13 +0000 |
| parents | 9c2a30bbaa0b |
| children | e347f09afbbb |
| files | extra/tsa/INDEX extra/tsa/aarmam.m |
| diffstat | 2 files changed, 246 insertions(+), 0 deletions(-) [+] |
line wrap: on
line diff
--- a/extra/tsa/INDEX Tue Jun 24 19:22:37 2003 +0000 +++ b/extra/tsa/INDEX Tue Jun 24 19:27:13 2003 +0000 @@ -20,14 +20,19 @@ ar_spa detrend flix + Multivariate stationary analysis mvar mvfilter mvfreqz arfit2 histo4 + Adaptive (time-varying) analysis aar + aarmam + adim + Conversions between forms ac2poly ac2rc
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/tsa/aarmam.m Tue Jun 24 19:27:13 2003 +0000 @@ -0,0 +1,241 @@ +function [z,e,REV,ESU,V,Z,SPUR] = aarmam(y, Mode, MOP, UC, z0, Z0, V0, W); +% Estimating Adaptive AutoRegressive-Moving-Average-and-mean model (includes mean term) +% [z,E,REV,ESU,V,Z,SPUR] = amarma(y, mode, MOP, UC, z0, Z0, V0, W); +% Estimates AAR parameters with Kalman filter algorithm +% y(t) = sum_i(a_i(t)*y(t-i)) + m(t) + e(t) + sum_i(b_i(t)*e(t-i)) +% +% State space model +% z(t) = G*z(t-1) + w(t) w(t)=N(0,W) +% y(t) = H*z(t) + v(t) v(t)=N(0,V) +% +% G = I, +% z = [m(t),a_1(t-1),..,a_p(t-p),b_1(t-1),...,b_q(t-q)]; +% H = [1,y(t-1),..,y(t-p),e(t-1),...,e(t-q)]; +% W = E{(z(t)-G*z(t-1))*(z(t)-G*z(t-1))'} +% V = E{(y(t)-H*z(t-1))*(y(t)-H*z(t-1))'} +% +% +% Input: +% y Signal (AR-Process) +% Mode determines the type of algorithm +% +% MOP Model order [m,p,q], default [0,10,0] +% m=1 includes the mean term, m=0 does not. +% p and q must be positive integers +% it is recommended to set q=0. +% UC Update Coefficient, default 0 +% z0 Initial state vector +% Z0 Initial Covariance matrix +% +% Output: +% z AR-Parameter +% E error process (Adaptively filtered process) +% REV relative error variance MSE/MSY +% +% REFERENCE(S): +% [1] A. Schloegl (2000), The electroencephalogram and the adaptive autoregressive model: theory and applications. +% ISBN 3-8265-7640-3 Shaker Verlag, Aachen, Germany. +% +% More references can be found at +% http://www.dpmi.tu-graz.ac.at/~schloegl/publications/ + +% Version 2.99 07 Aug 2002 +% Copyright (c) 1998-2002 by Alois Schloegl <a.schloegl@ieee.org> +% +% This library is free software; you can redistribute it and/or +% modify it under the terms of the GNU Library General Public +% License as published by the Free Software Foundation; either +% version 2 of the License, or (at your option) any later version. +% +% This library is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +% Library General Public License for more details. +% +% You should have received a copy of the GNU Library General Public +% License along with this library; if not, write to the +% Free Software Foundation, Inc., 59 Temple Place - Suite 330, +% Boston, MA 02111-1307, USA. + +% 12.04.1999 ESU included +% 19.10.2000 aMode 13,14 included +% NaN handling +% 06.07.2002 Docu improved, included into TSA +% 11.07.2002 nanmean replaced by mean +% 07.08.2002 cosmetic changes + +%#realonly +%#inbounds + + +[nc,nr]=size(y); + +if nargin<2 Mode=0; +elseif ischar(Mode) Mode=bin2dec(Mode); +elseif isnan(Mode) return; end; +if nargin<3 MOP=[0,10,0]; end; +if length(MOP)==0, m=0;p=10; q=0; MOP=p; +elseif length(MOP)==1, m=0;p=MOP(1); q=0; MOP=p; +elseif length(MOP)==2, fprintf(1,'Error AMARMA: MOP is ambiguos\n'); +elseif length(MOP)>2, m=MOP(1); p=MOP(2); q=MOP(3);MOP=m+p+q; +end; + +if prod(size(Mode))>1 + aMode=Mode(1); + eMode=Mode(2); +end; +%fprintf(1,['a' int2str(aMode) 'e' int2str(eMode) ' ']); + + +e=zeros(nc,1); +V=zeros(nc,1);V(1)=V0; +T=zeros(nc,1); +ESU=zeros(nc,1)+nan; +SPUR=zeros(nc,1)+nan; +z=z0(ones(nc,1),:); + +arc=poly((1-UC*2)*[1;1]);b0=sum(arc); % Whale forgetting factor for Mode=258,(Bianci et al. 1997) + +dW=UC/MOP*eye(MOP); % Schloegl + + +%------------------------------------------------ +% First Iteration +%------------------------------------------------ + +H = zeros(MOP,1); +if m, + H(1) = 1;%M0; + if m~=1, + fprintf(2,'Warning AARMAM: m must be 0 or 1\n'); + return; + end; +end; +if (p<0) | (q<0) | (round(p)~=p) | (round(q)~=q), + fprintf(2,'Error AARMAM: p and q must be positive integers\n'); + return; +end; + +E = 0; +Z = Z0; +zt= z0; + +A1 = zeros(MOP); A2 = A1; + +y_1=0; + +%------------------------------------------------ +% Update Equations +%------------------------------------------------ + +for t=1:nc, + + % make measurement matrix + if 0, + if t>1, + y_1 = y(t-1); + end; + H=[1; y_1; H(m+(1:p-1)'); E(1:min(1,q-1)) ; H(p+m+(1:q-1)')]; % shift y and e + + else % this seem to be slightly faster + if t<=p, H(m+(1:t-1)) = y(t-1:-1:1); % Autoregressive + else H(m+(1:p)) = y(t-1:-1:t-p); + end; + + if t<=q, H(m+p+(1:t-1)) = e(t-1:-1:1); % Moving Average + else H(m+p+(1:q)) = e(t-1:-1:t-q); + end; + end; + + % Prediction Error + E = y(t) - zt*H; + e(t) = E; + + if ~isnan(E), + E2 = E*E; + AY = Z*H; + + ESU(t) = H'*AY; + + if eMode==1 + V0 = V(t-1); + V(t) = V0*(1-UC)+UC*E2; + elseif eMode==2 + V0 = 1; + V(t) = V0; %V(t-1)*(1-UC)+UC*E2; + elseif eMode==3 + V0 = 1-UC; + V(t) = V0; %(t-1)*(1-UC)+UC*E2; + elseif eMode==4 + V0 = V0*(1-UC)+UC*E2; + V(t) = V0; + elseif eMode==5 + V(t)=V0; + %V0 = V0; + elseif eMode==6 + if E2>ESU(t), + V0=(1-UC)*V0+UC*(E2-ESU(t)); + end; + V(t)=V0; + elseif eMode==7 + V0=V(t); + if E2>ESU(t) + V(t) = (1-UC)*V0+UC*(E2-ESU(t)); + else + V(t) = V0; + end; + elseif eMode==8 + V0=0; + V(t) = V0; % (t-1)*(1-UC)+UC*E2; + end; + + k = AY / (ESU(t) + V0); % Kalman Gain + zt = zt + k'*E; + %z(t,:) = zt; + + if aMode==2 + T(t)=(1-UC)*T(t-1)+UC*(E2-Q(t))/(H'*H); % Roberts I 1998 + Z=Z*V(t-1)/Q(t); + if T(t)>0 W=T(t)*eye(MOP); else W=zeros(MOP);end; + elseif aMode==5 + Q_wo = (H'*C*H + V(t-1)); % Roberts II 1998 + T(t)=(1-UC)*T(t-1)+UC*(E2-Q_wo)/(H'*H); + if T(t)>0 W=T(t)*eye(MOP); else W=zeros(MOP); end; + elseif aMode==6 + T(t)=(1-UC)*T(t-1)+UC*(E2-Q(t))/(H'*H); + Z=Z*V(t)/Q(t); + if T(t)>0 W=T(t)*eye(MOP); else W=zeros(MOP); end; + elseif aMode==11 + %Z = Z - k*AY'; + W = sum(diag(Z))*dW; + elseif aMode==12 + W = UC*UC*eye(MOP); + elseif aMode==13 + W = UC*diag(diag(Z)); + elseif aMode==14 + W = (UC*UC)*diag(diag(Z)); + elseif aMode==15 + W = sum(diag(Z))*dW; + elseif aMode==16 + W = UC*eye(MOP); % Schloegl 1998 + %elseif aMode==17 + %W=W; + end; + + Z = Z - k*AY'; % Schloegl 1998 + else + + V(t) = V0; + + end; + + if any(any(isnan(W))), W=UC*Z;end; + + z(t,:) = zt; + Z = Z + W; % Schloegl 1998 + SPUR(t) = trace(Z); +end; + +REV = mean(e.*e)/mean(y.*y); +if any(~isfinite(Z(:))), REV=inf; end; +