0001 function varargout = ltpda_loglike(varargin)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024 ALGONAME = mfilename;
0025 VERSION = '$Id: ltpda_loglike.m,v 1.1 2008/03/11 16:52:56 adrien Exp $';
0026 CATEGORY = 'STATESPACE';
0027 display(['starting ' ALGONAME]);
0028
0029 if not(isempty(varargin))
0030 if isequal( varargin{1}, 'Version')
0031 varargout = VERSION;
0032 return;
0033 elseif isequal(varargin{1}, 'Params')
0034 varargout = plist('lamba', 1.025 );
0035 return;
0036 elseif isequal(varargin{1}, 'Category')
0037 varargout = CATEGORY;
0038 return;
0039 end
0040 end
0041
0042 crit_out = plist();
0043 ao_exp = varargin{1};
0044 ao_est = varargin{2};
0045 options = varargin{3};
0046 for i_ao_exp = 1:length(ao_exp)
0047 for i_ao_est = 1:length(ao_est)
0048 for i_options = 1:length(options)
0049 options(i_options) = combine( options(i_options), plist('lamba', 1.025 ));
0050
0051 lambda = find(options(i_options),'lambda');
0052 Y = transpose(ao_exp(i_ao_exp).data.y);
0053 Y_est = transpose(ao_est(i_ao_est).data.y);
0054 t = size(Y,2);
0055 dout = size(Y,1);
0056
0057 if not(isequal(size(Y),size(Y_est)))
0058 msg = 'error because input aos should be the same size';
0059 error(msg);
0060 end
0061
0062 NT = round(log(2)/log(lambda))+1;
0063
0064 NNT = floor(log(t)/log(2));
0065
0066 Rcrush = zeros(dout,dout,NT,NNT);
0067 X1 =Y;
0068 longueur = zeros(1,NNT+1);
0069 longueur(1) = t;
0070 Tcrush = zeros(NNT,NT);
0071 for i = 1:NNT
0072 for j= 1:NT
0073 Tcrush(i, j)= 2^i*j;
0074 if (longueur(i)-Tcrush(i,j)>0)&&( j==1||2*j>NT )
0075 Ricrush = norm2lag0(X1, Tcrush(i,j));
0076 Rcrush(:,:,i,j)=(Ricrush+1e-13*eye(dout))^(-1);
0077 end
0078 end
0079 X1 =crush(X1);
0080 longueur(i+1) = size(X1,2);
0081 end
0082 Y_err_crush = Y-Y_est;
0083
0084 criterioncrush = 0;
0085 for i = 1:NNT
0086 for j= 1:NT
0087 if (longueur(i)-Tcrush(i,j)>0)&&(j==1||2*j>NT)
0088 criterioncrush = criterioncrush + norm(norm2lag(Y_err_crush, Rcrush(:,:,i,j), Tcrush(i,j)));
0089 end
0090 end
0091 Y_err_crush =crush(Y_err_crush);
0092 end
0093
0094 crit_out(i_ao_exp, i_ao_est, i_options )= plist('criterion', criterioncrush, 'lambda', lambda);
0095 end
0096 end
0097 end
0098 varargout = {crit_out};
0099 end
0100
0101
0102 function a = crush(b)
0103
0104 n = size(b,1);
0105 N = floor(size(b,2)/2);
0106 a = zeros(n,N);
0107 for i=1:N
0108 for j=1:n
0109 a(j,i) = b(j, 2*i) + b(j, 2*i-1);
0110 end
0111 end
0112 end
0113
0114 function p = norm2lag(y,S,lag)
0115
0116 [h,l] = size(y);
0117 p = zeros(h,h);
0118 for j = 1:h
0119 for k = 1:h
0120 x1 = y(k,1+lag:l);
0121 x2 = transpose(y(j,1:l-lag));
0122 p(j,k) = p(j,k) + S(j,k)*(x1*x2);
0123 end
0124 end
0125 end
0126
0127
0128
0129 function p = norm2lag0(y,lag)
0130
0131 [h,l] = size(y);
0132 p = zeros(h,h);
0133 for i=1:h
0134 for j =1:h
0135 if lag<l
0136 p(i,j) = y(i,1+lag:l)*transpose(y(j,1:l-lag));
0137 end
0138 end
0139 end
0140 end