X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?p=CreaPhase.git;a=blobdiff_plain;f=octave_packages%2Ftsa-4.2.4%2Farfit2.m;fp=octave_packages%2Ftsa-4.2.4%2Farfit2.m;h=9e9efa39a72c51a24a4db3a54c1ec6f59f47bec1;hp=0000000000000000000000000000000000000000;hb=f5f7a74bd8a4900f0b797da6783be80e11a68d86;hpb=1705066eceaaea976f010f669ce8e972f3734b05

diff --git a/octave_packages/tsa-4.2.4/arfit2.m b/octave_packages/tsa-4.2.4/arfit2.m
new file mode 100644
index 0000000..9e9efa3
--- /dev/null
+++ b/octave_packages/tsa-4.2.4/arfit2.m
@@ -0,0 +1,143 @@
+function [w, MAR, C, sbc, fpe, th] = arfit2(Y, pmin, pmax, selector, no_const)
+% ARFIT2 estimates multivariate autoregressive parameters
+% of the MVAR process Y
+%
+%   Y(t,:)' = w' + A1*Y(t-1,:)' + ... + Ap*Y(t-p,:)' + x(t,:)'
+%
+% ARFIT2 uses the Nutall-Strand method (multivariate Burg algorithm) 
+% which provides better estimates the ARFIT [1], and uses the 
+% same arguments. Moreover, ARFIT2 is faster and can deal with 
+% missing values encoded as NaNs. 
+%
+% [w, A, C, sbc, fpe] = arfit2(v, pmin, pmax, selector, no_const)
+%
+% INPUT: 
+%  v		data - each channel in a column
+%  pmin, pmax 	minimum and maximum model order
+%  selector	'fpe' or 'sbc' [default] 
+%  no_const	'zero' indicates no bias/offset need to be estimated 
+%		in this case is w = [0, 0, ..., 0]'; 
+%
+% OUTPUT: 
+%  w		mean of innovation noise
+%  A		[A1,A2,...,Ap] MVAR estimates	
+%  C		covariance matrix of innovation noise
+%  sbc, fpe	criteria for model order selection 
+%
+% see also: ARFIT, MVAR
+%
+% REFERENCES:
+%  [1] A. Schloegl, 2006, Comparison of Multivariate Autoregressive Estimators.
+%       Signal processing, p. 2426-9.
+%  [2] T. Schneider and A. Neumaier, 2001. 
+%	Algorithm 808: ARFIT-a Matlab package for the estimation of parameters and eigenmodes 
+%	of multivariate autoregressive models. ACM-Transactions on Mathematical Software. 27, (Mar.), 58-65.
+
+%       $Id: arfit2.m 9609 2012-02-10 10:18:00Z schloegl $
+%	Copyright (C) 1996-2005,2008,2012 by Alois Schloegl <alois.schloegl@ist.ac.at>	
+%       This is part of the TSA-toolbox. See also 
+%       http://pub.ist.ac.at/~schloegl/matlab/tsa/
+%       http://octave.sourceforge.net/
+%       http://biosig.sourceforge.net/
+
+%    This program is free software: you can redistribute it and/or modify
+%    it under the terms of the GNU General Public License as published by
+%    the Free Software Foundation, either version 3 of the License, or
+%    (at your option) any later version.
+%
+%    This program 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 General Public License for more details.
+%
+%    You should have received a copy of the GNU General Public License
+%    along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+
+%%%%% checking of the input arguments was done the same way as ARFIT
+if (pmin ~= round(pmin) || pmax ~= round(pmax))
+        error('Order must be integer.');
+end
+if (pmax < pmin)
+        error('PMAX must be greater than or equal to PMIN.')
+end
+
+% set defaults and check for optional arguments
+if (nargin == 3)              	% no optional arguments => set default values
+        mcor       = 1;         % fit intercept vector
+        selector   = 'sbc';	% use SBC as order selection criterion
+elseif (nargin == 4)          	% one optional argument
+        if strcmp(selector, 'zero')
+                mcor     = 0;   % no intercept vector to be fitted
+                selector = 'sbc';	% default order selection 
+        else
+                mcor     = 1;	% fit intercept vector
+        end
+elseif (nargin == 5)		% two optional arguments
+        if strcmp(no_const, 'zero')
+                mcor     = 0;   	% no intercept vector to be fitted
+        else
+                error(['Bad argument. Usage: ', '[w,A,C,SBC,FPE,th]=AR(v,pmin,pmax,SELECTOR,''zero'')'])
+        end
+end
+
+
+
+%%%%% Implementation of the MVAR estimation 
+[N,M]=size(Y);
+    
+if mcor,
+        [m,N] = sumskipnan(Y,1);                    % calculate mean 
+        m = m./N;
+	Y = Y - repmat(m,size(Y)./size(m));    % remove mean    
+end;
+
+[MAR,RCF,PE] = mvar(Y, pmax, 2);   % estimate MVAR(pmax) model
+
+N = min(N);
+
+%if 1;nargout>3;
+ne = N-mcor-(pmin:pmax);
+for p=pmin:pmax, 
+        % Get downdated logarithm of determinant
+        logdp(p-pmin+1) = log(det(PE(:,p*M+(1:M))*(N-p-mcor))); 
+end;
+
+% Schwarz's Bayesian Criterion
+sbc = logdp/M - log(ne) .* (1-(M*(pmin:pmax)+mcor)./ne);
+
+% logarithm of Akaike's Final Prediction Error
+fpe = logdp/M - log(ne.*(ne-M*(pmin:pmax)-mcor)./(ne+M*(pmin:pmax)+mcor));
+
+% Modified Schwarz criterion (MSC):
+% msc(i) = logdp(i)/m - (log(ne) - 2.5) * (1 - 2.5*np(i)/(ne-np(i)));
+
+% get index iopt of order that minimizes the order selection 
+% criterion specified by the variable selector
+if strcmpi(selector,'fpe'); 
+    [val, iopt]  = min(fpe); 
+else %if strcmpi(selector,'sbc'); 
+    [val, iopt]  = min(sbc); 
+end; 
+
+% select order of model
+popt = pmin + iopt-1; % estimated optimum order 
+
+if popt<pmax, 
+        [MAR, RCF, PE] = mvar(Y, popt, 2);
+end;
+%end
+
+C = PE(:,size(PE,2)+(1-M:0));
+
+if mcor,
+        I = eye(M);        
+        for k = 1:popt,
+                I = I - MAR(:,k*M+(1-M:0));
+        end;
+        w = -I*m';
+else
+        w = zeros(M,1);
+end;
+
+if nargout>5,	th=[];	fprintf(2,'Warning ARFIT2: output TH not defined\n'); end;