Home > m > timetools > statespacefunctions > ltpda_loglike.m

ltpda_loglike

PURPOSE ^

Log likelihood criterion

SYNOPSIS ^

function varargout = ltpda_loglike(varargin)

DESCRIPTION ^

Log likelihood criterion
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 DESCRIPTION: ltpda_loglike computes the "crunched" log likelihood
 criterion of two time series.

 CALL: crit_out = simulate_discrete(ao_exp, ao_est, options)

 INPUTS: ao_exp - (array of) ao containing experimental data stream
 ao_est - (array of) ao containing estimated (=filered) data stream
 options - computation options

 OUTPUTS: crit_out - (array of) plist containing output criterion value

 &&VERSION: $Id: ltpda_loglike.m,v 1.1 2008/03/11 16:52:56 adrien Exp $

 HISTORY: 07-03-2008 A Grynagier

 TO DO: tests on frequency of data, evaluate processing time... 
 what if time sample are not evenly scaled
 cut if time series are not the same length
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function varargout = ltpda_loglike(varargin)
0002 %Log likelihood criterion
0003 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0004 %
0005 % DESCRIPTION: ltpda_loglike computes the "crunched" log likelihood
0006 % criterion of two time series.
0007 %
0008 % CALL: crit_out = simulate_discrete(ao_exp, ao_est, options)
0009 %
0010 % INPUTS: ao_exp - (array of) ao containing experimental data stream
0011 % ao_est - (array of) ao containing estimated (=filered) data stream
0012 % options - computation options
0013 %
0014 % OUTPUTS: crit_out - (array of) plist containing output criterion value
0015 %
0016 % &&VERSION: $Id: ltpda_loglike.m,v 1.1 2008/03/11 16:52:56 adrien Exp $
0017 %
0018 % HISTORY: 07-03-2008 A Grynagier
0019 %
0020 % TO DO: tests on frequency of data, evaluate processing time...
0021 % what if time sample are not evenly scaled
0022 % cut if time series are not the same length
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             % retriveing data
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             % checking data size is correct
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             % number of time lags per crusched series
0062             NT = round(log(2)/log(lambda))+1;
0063             % number of crunchings
0064             NNT = floor(log(t)/log(2));
0065             % evaluating covariance of experimental signal
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             % evaluating criterion
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             %% construct output analysis object
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 % data reduction
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 % dot product with weight S and time lag
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 % dot product with time lag
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

Generated on Mon 31-Mar-2008 13:54:54 by m2html © 2003