root/applications/bdmtoolbox/tutorial/mpdm/dist_estim_example.m @ 760

% name random variables
y = RV({'y'},1);
z = RV({'z'},1);
u = RV({'u'},1);

% create f(y_t| z_{t-1}, z_{t-2})
fy.class = 'mlnorm<ldmat>';
fy.rv    = y;
fy.rvc   = RVtimes([z,z], [0, -1]);
fy.A     = [1.8, -0.9];
fy.const = 0.2;
fy.R     = 1e-2;

fz.class = 'mlnorm<ldmat>';
fz.rv    = z;
fz.rvc   = RVtimes([u,u], [0, -1]);
fz.A     = [1.8, -0.9];
fz.const = -0.2;
fz.R     = 1e-2;

fu.class = 'enorm<ldmat>';
fu.rv    = u;
fu.mu    = [0];
fu.R     = 0.001;

f.class = 'mprod';
f.pdfs  = {fy,fz,fu};

DS.class = 'PdfDS';
DS.pdf = f;

% create ARX estimator
A1.class = 'ARX';
A1.rv = y;
A1.rgr = RVtimes([z,z],[0,-1]) ; % correct structure is {y,y}
A1.options ='logfull,logll';
A1.frg = 0.95;
A1.constant = 1;

A2=A1;
A2.rv  = z;
A2.rgr = RVtimes([u,u],[0,-1]) ; % correct structure is {y,y}

[M,Set]=estimator(DS,{A1,A2},struct('ndat',100));

%% post-merging

Merger.class='merger_mix';
Merger.method='lognormal';
Merger.beta=1.2;
Merger.ncoms=20;
Merger.stop_niter=10;
Merger.effss_coef=0.9;

support.class='discrete_support';
support.epdf= struct('class','enorm<ldmat>','mu',[1,1,1,1,0.1],'R',eye(5));
support.npoints=[100];



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% plot results
ndat = size(M.DS_u,1);

true_theta1 = [fy.A fy.const];
true_theta2 = [fz.A fy.const];
true_R = fy.R;

A1_mean = zeros(4,ndat);
A2_mean = zeros(4,ndat);
MG_mean = zeros(5,ndat);

for t=1:ndat
    f1=Set.Est0_apost{t};
    f1.rv = RV({'a','b','c','r'});
    f2=Set.Est1_apost{t};
    f2.rv = RV({'a','b','d','r'});
    
    [res,ftilde] = merger({f1,f2},support,Merger);
    
    A1_mean(:,t) = epdf_mean(f1);
    A2_mean(:,t) = epdf_mean(f2);
    MG_mean(:,t) = epdf_mean(ftilde);
end;


figure(1);
subplot(3,1,1);
plot(M.DS_y);
title('y');

subplot(3,1,2);
plot(M.DS_z);
title('z')

subplot(3,1,2);
plot(M.DS_u);
title('u')


figure(2)
subplot(2,2,1);
plot(A1_mean');
subplot(2,2,2);
plot(A2_mean');
subplot(2,2,3);
plot(MG_mean');