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);
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
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)
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
28 % AIC > FPE > *MDL* > PHI > SBC > CAT ~ BIC
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)"
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/
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.
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.
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/>.
70 NC=(lc-sum(isnan(e)')')*(NC<lc) + NC.*(NC>=lc); % first part
71 %end;% Pmax=min([100 N/3]); end;
73 %NC=(lc-sum(isnan(e)')')*(NC<lc) + NC.*(NC>=lc); % first part
81 e = e./e(:,ones(1,lc));
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,
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;
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;
113 HAR(2:lc)=-(N-m).*log((N-m).*E(m+1)./(N-m+1)./E(m));
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;
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;
128 %CATcrit= (cumsum(1./E(m+1))/N-1./E(m+1));
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;
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;
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;
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)];
146 C=[FPE;AIC;BIC;SBC;MDL;CATcrit;PHI;JEW;HAR(:)']';