octave_packages/tsa-4.2.4/selmo.m

   1 function [FPE,AIC,BIC,SBC,MDL,CATcrit,PHI,optFPE,optAIC,optBIC,optSBC,optMDL,optCAT,optPHI,p,C]=selmo(e,NC);
   2 % Model order selection of an autoregrssive model
   3 % [FPE,AIC,BIC,SBC,MDL,CAT,PHI,optFPE,optAIC,optBIC,optSBC,optMDL,optCAT,optPHI]=selmo(E,N);
   4 %
   5 % E     Error function E(p)
   6 % N     length of the data set, that was used for calculating E(p)
   7 % show  optional; if given the parameters are shown
   8 %
   9 % FPE   Final Prediction Error (Kay 1987, Wei 1990, Priestley 1981  -> Akaike 1969)
  10 % AIC   Akaike Information Criterion (Marple 1987, Wei 1990, Priestley 1981 -> Akaike 1974)
  11 % BIC   Bayesian Akaike Information Criterion (Wei 1990, Priestley 1981 -> Akaike 1978,1979)
  12 % CAT   Parzen's CAT Criterion (Wei 1994 -> Parzen 1974)
  13 % MDL   Minimal Description length Criterion (Marple 1987 -> Rissanen 1978,83)
  14 % SBC   Schwartz's Bayesian Criterion (Wei 1994; Schwartz 1978)
  15 % PHI   Phi criterion (Pukkila et al. 1988, Hannan 1980 -> Hannan & Quinn, 1979)
  16 % HAR   Haring G. (1975)
  17 % JEW   Jenkins and Watts (1968)
  18 %
  19 % optFPE        order where FPE is minimal
  20 % optAIC        order where AIC is minimal
  21 % optBIC        order where BIC is minimal
  22 % optSBC        order where SBC is minimal
  23 % optMDL        order where MDL is minimal
  24 % optCAT        order where CAT is minimal
  25 % optPHI        order where PHI is minimal
  26 %
  27 % usually is
  28 % AIC > FPE > *MDL* > PHI > SBC > CAT ~ BIC
  29 %
  30 % REFERENCES:
  31 %  P.J. Brockwell and R.A. Davis "Time Series: Theory and Methods", 2nd ed. Springer, 1991.
  32 %  S. Haykin "Adaptive Filter Theory" 3ed. Prentice Hall, 1996.
  33 %  M.B. Priestley "Spectral Analysis and Time Series" Academic Press, 1981.
  34 %  C.E. Shannon and W. Weaver "The mathematical theory of communication" University of Illinois Press, Urbana 1949 (reprint 1963).
  35 %  W.S. Wei "Time Series Analysis" Addison Wesley, 1990.
  36 %  Jenkins G.M. Watts D.G "Spectral Analysis and its applications", Holden-Day, 1968.
  37 %  G. Haring  "Über die Wahl der optimalen Modellordnung bei der Darstellung von stationären Zeitreihen mittels Autoregressivmodell als Basis der Analyse von EEG - Biosignalen mit Hilfe eines Digitalrechners", Habilitationschrift - Technische Universität Graz, Austria, 1975.
  38 %                  (1)"About selecting the optimal model at the representation of stationary time series by means of an autoregressive model as basis of the analysis of EEG - biosignals by means of a digital computer)"
  39 %
  40
  41 %       $Id: selmo.m 9609 2012-02-10 10:18:00Z schloegl $
  42 %       Copyright (C) 1997-2002,2008,2012 by Alois Schloegl <alois.schloegl@ist.ac.at>
  43 %       This is part of the TSA-toolbox. See also
  44 %       http://pub.ist.ac.at/~schloegl/matlab/tsa/
  45 %       http://octave.sourceforge.net/
  46 %       http://biosig.sourceforge.net/
  47 %
  48 %    This program is free software: you can redistribute it and/or modify
  49 %    it under the terms of the GNU General Public License as published by
  50 %    the Free Software Foundation, either version 3 of the License, or
  51 %    (at your option) any later version.
  52 %
  53 %    This program is distributed in the hope that it will be useful,
  54 %    but WITHOUT ANY WARRANTY; without even the implied warranty of
  55 %    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  56 %    GNU General Public License for more details.
  57 %
  58 %    You should have received a copy of the GNU General Public License
  59 %    along with this program.  If not, see <http://www.gnu.org/licenses/>.
  60
  61 [lr,lc]=size(e);
  62 if (lr>1) && (lc>1),
  63         p=zeros(lr+1,9)+NaN;
  64 else
  65         p=zeros(1,9)+NaN;
  66 end;
  67
  68 if nargin<2
  69         NC=lc*ones(lr,1);
  70         NC=(lc-sum(isnan(e)')')*(NC<lc) + NC.*(NC>=lc); % first part
  71 %end;% Pmax=min([100 N/3]); end;
  72         %if NC<lc N=lc; end;
  73         %NC=(lc-sum(isnan(e)')')*(NC<lc) + NC.*(NC>=lc); % first part
  74 else
  75         % NC=NC;
  76 end;
  77
  78 M=lc-1;
  79 m=0:M;
  80
  81 e = e./e(:,ones(1,lc));
  82
  83 for k=0:lr,
  84         if k>0, %
  85                 E=e(k,:);
  86                 N=NC(k);
  87         elseif lr>1
  88                 tmp = e;%(NC>0,:);
  89                 tmp(isnan(tmp)) = 0;
  90                 E = sum(tmp.*(NC*ones(1,lc)))/sum(NC); % weighted average, weigths correspond to number of valid (not missing) values
  91                 N = sum(NC)./sum(NC>0); % corresponding number of values,
  92         else
  93                 E = e;
  94                 N = NC;
  95         end;
  96 FPE = E.*(N+m)./(N-m);  %OK
  97                 optFPE=find(FPE==min(FPE))-1;   %optimal order
  98         if isempty(optFPE), optFPE=NaN; end;
  99 AIC = N*log(E)+2*m;     %OK
 100         optAIC=find(AIC==min(AIC))-1;   %optimal order
 101         if isempty(optAIC), optAIC=NaN; end;
 102 AIC4=N*log(E)+4*m;      %OK
 103         optAIC4=find(AIC4==min(AIC4))-1;        %optimal order
 104         if isempty(optAIC4), optAIC4=NaN; end;
 105
 106 m=1:M;
 107 BIC=[ N*log(E(1)) N*log(E(m+1)) - (N-m).*log(1-m/N) + m*log(N) + m.*log(((E(1)./E(m+1))-1)./m)];
 108 %BIC=[ N*log(E(1)) N*log(E(m+1)) - m + m*log(N) + m.*log(((E(1)./E(m+1))-1)./m)];
 109 %m=0:M; BIC=N*log(E)+m*log(N);          % Hannan, 1980 -> Akaike, 1977 and Rissanen 1978
 110         optBIC=find(BIC==min(BIC))-1;   %optimal order
 111         if isempty(optBIC), optBIC=NaN; end;
 112
 113 HAR(2:lc)=-(N-m).*log((N-m).*E(m+1)./(N-m+1)./E(m));
 114         HAR(1)=HAR(2);
 115         optHAR=min(find(HAR<=(min(HAR)+0.2)))-1;        %optimal order
 116 %       optHAR=find(HAR==min(HAR))-1;   %optimal order
 117         if isempty(optHAR), optHAR=NaN; end;
 118
 119 m=0:M;
 120 SBC = N*log(E)+m*log(N);
 121         optSBC=find(SBC==min(SBC))-1;   %optimal order
 122         if isempty(optSBC), optSBC=NaN; end;
 123 MDL = N*log(E)+log(N)*m;
 124         optMDL=find(MDL==min(MDL))-1;   %optimal order
 125         if isempty(optMDL), optMDL=NaN; end;
 126
 127 m=0:M;
 128 %CATcrit= (cumsum(1./E(m+1))/N-1./E(m+1));
 129 E1=N*E./(N-m);
 130 CATcrit= (cumsum(1./E1(m+1))/N-1./E1(m+1));
 131         optCAT=find(CATcrit==min(CATcrit))-1;   %optimal order
 132         if isempty(optCAT), optCAT=NaN; end;
 133
 134 PHI = N*log(E)+2*log(log(N))*m;
 135         optPHI=find(PHI==min(PHI))-1;   %optimal order
 136         if isempty(optPHI), optPHI=NaN; end;
 137
 138 JEW = E.*(N-m)./(N-2*m-1);      % Jenkins-Watt
 139         optJEW=find(JEW==min(JEW))-1;   %optimal order
 140         if isempty(optJEW), optJEW=NaN; end;
 141
 142 % in case more than 1 minimum is found, the smaller model order is returned;
 143 p(k+1,:) = [optFPE(1), optAIC(1), optBIC(1), optSBC(1), optCAT(1), optMDL(1), optPHI(1), optJEW(1), optHAR(1)];
 144
 145 end;
 146 C=[FPE;AIC;BIC;SBC;MDL;CATcrit;PHI;JEW;HAR(:)']';