% load data created by the MpdfDS_example
clear all
load pdfds_results

DS.class   = 'MemDS';
DS.Data    = Data;
DS.drv     = drv;

%%%%%% ARX estimator

A1.class = 'ARX';
A1.yrv = y;
A1.rv = RV({'theta','r'},[3,1]);
A1.rgr = RVtimes([y,u],[-3,-1]) ; % correct structure is {y,y}
A1.log_level = 'logbounds,logevidence';

A2.class = 'ARX';
A2.yrv = u;
%A2.rv = RV({'theta','r'},[3,1]);
A2.rgr = RVtimes([y,u],[-3,-1]) ; % correct structure is {y,y}
A2.log_level = 'logbounds,logevidence';

P.class='ProdBM';
P.BMs={A1,A2};
P.log_level = 'logbounds';

[M,Apost]=estimator(DS,{P});

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

subplot(1,2,1);
hold off
plotestimates(true_theta, ...
    M.Est0_apost_mean_theta, ...
    M.Est0_apost_lbound_theta, ...
    M.Est0_apost_ubound_theta);
set(gca,'YLim',[-1.5,1]);

subplot(1,2,2);
hold off
plotestimates(true_R, ...
    M.Est0_apost_mean_r, ...
    M.Est0_apost_lbound_r, ...
    M.Est0_apost_ubound_r);

title('Variance parameters r')
