tfar_llfa.m

用于模拟时变非平稳的ARMA过程

function l= tfar_llfa(AB, Pyy, re_im)
% function l= tfar_llfa(AB, Pyy, re_im)
%   This file is part of the TFPM toolbox v1.0 (c)
%   [e-mail] and underlies the GPL.
% 
% Computation of the NEGATIVE asymptotic log-likelihood function of
% Gaussian TFAR(M, L; N; 1/2) models with parameters AB= [Aml
% B0l]. Takes care of re_im! 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;tfpm;
N    =  64;
MAR  =   2;
LAR  =   2;
MMA  =   0;
LMA  = LAR;
re_im= 'r';
mo_no= 'm';
tfpm_file_gen;
%-------------
beta = 1/2;
y= tfarma_gen(randn(N, 1), Aml, Bml, beta);
Ryy= corr_est(y, y, -1, 1/2);
Ayy= nm_to_ml(Ryy);
[Psi, Mask]= tf_multiwin(N, MAR, LAR, 5, 1, 1);
Pyy= real(nm_to_nk(ml_to_nm(Ayy.*conj(Psi))));

Aml0= Aml;
B0l0= Bml;
Amlref= Aml;
B0lref= Bml;
if(re_im=='r')
   y= real(y);
   Aml0= param_hermite(Aml0);
   Aml0= [real(Aml0(1:LAR+1, 2:end)); imag(Aml0(LAR+2:end, 2:end))];
else
   Aml0= [real(Aml0(:, 2:end)) imag(Aml0(:, 2:end))];
end;
B0l0= param_hermite(B0l0);
B0l0= [real(B0l0(1:LAR+1, :)); imag(B0l0(LAR+2:end, :))];
AB= [Aml0 B0l0];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Dimensions:
N= size(Pyy, 1);

% Check for re_im:
if(re_im=='r')
   [MAR, LAR]= param_dim(AB);
% Split up AB, add a_{0, l}:
   Aml= [[zeros(LAR, 1); 1; zeros(LAR, 1)] AB(:, 1:end-1)];
   Aml= [Aml(1:LAR+1, :); flipud(Aml(1:LAR, :))] + j*[flipud(-Aml(LAR+2:end, :)); zeros(1, MAR+1); Aml(LAR+2:end, :)];
else
   [MAR, LAR]= param_dim(AB);MAR= MAR/2;
   Aml= [[zeros(LAR, 1); 1; zeros(LAR, 1)] AB(:, 1:MAR)] + j*[zeros(2*LAR+1, 1) AB(:, MAR+1:2*MAR)];
end;
B0l= AB(:, end);
B0l= [B0l(1:LAR+1); flipud(B0l(1:LAR))] + j*[flipud(-B0l(LAR+2:end)); 0; B0l(LAR+2:end)];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Amlref-Aml
B0lref-B0l
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Wyy= tfarma_wvsp(Aml, B0l, N, 1/2);

l=  -(N*log(2*pi) + 1/N*sum(sum(log(Wyy)+Pyy./Wyy)))/2;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Slow, exact reference:
Anm= param_expand(Aml, N);
Bn0= real(param_expand(B0l, N));
H= tvarma_impr(Anm, Bn0);
R= H*H';
lref= (N*log(2*pi) + log(real(det(R))) + real(y'*inv(R)*y))/2;
[l lref]
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%