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

% 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.log_level ='logfull,logevidence';
A1.rv_param = RV({'a','b','c','r'});
A1.frg = 0.95;
A1.constant = 1;

A2=A1;
A2.rv  = z;
A2.rv_param = RV({'a','b','c','r'});
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;
%Merger.dbg_file = 'm.it';

M2.class = 'merger_base';
M2.method='lognormal';
M2.beta=1.2;

support.class='discrete_support';
support.epdf= struct('class','enorm<ldmat>','mu',[2,-1,0,0.1],'R',eye(4));
support.npoints=[300];

support2.class='rectangular_support';
support2.ranges= {[0,5],[-2,2],[-2,2],[0.001,3]};
support2.gridsizes = [10,10,10,10];


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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(4,ndat);

for t=1:ndat
    f1=Set.Est0_apost{t};
    f2=Set.Est1_apost{t};
    if t>5
%        support.epdf = f1;
    end

    [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)
plot(A1_mean',':');
hold on;
plot(A2_mean','--');
plot(MG_mean');