x = RV('x',2);
fx.class = 'mlnorm<ldmat>';
fx.A = [1 -0.9; 1 -0.3];
fx.const = [0;0];
fx.R = [0.1 0.05; 0.05 0.1];
fx.rv = x;
fx.rvc = RVtimes(x,-1);

y = RV('y',2);
fy.class = 'mlnorm<ldmat>';
fy.A = [1 0.5; -0.3 0];
fy.const = [0;0];
fy.R = 0.1*eye(2);
fy.rv = y;
fy.rvc = x;


%%%%%% Data generator
DS.class   = 'PdfDS';
DS.pdf.class  = 'mprod';
DS.pdf.pdfs   = {fy,fx};

Dt = simulator(DS);
%%%%%% ARX estimator conditioned on frg

A1.class = 'ARX';
A1.rv = y;
A1.rgr = x ; 
A1.constant = 0;
A1.options ='logbounds,logll';
A1.frg = 0.98;
A1.name = 'A1';

%%%%%% Combining estimators in Marginalized particle filter
E.class = 'MPF';
E.BM = A1;                         % ARX is the analytical part
E.parameter_pdf = fx;              % Random walk is the parameter evolution model
E.res_threshold = 1.0;             % resampling parameter
E.n = 1000;                         % number of particles
E.prior.class = 'enorm<ldmat>';    % prior on non-linear part
E.prior.mu    = [0 0];             % 
E.prior.R     = 10*eye(2);         % 
E.options ='logbounds,logll';
E.name = 'MPF';

exp.ndat = 1000;

M=estimator(DS,{E},exp);

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


figure(1);
subplot(4,1,1);
plotestimates(fy.A(1,1)*ones(ndat,1), M.MPF_apost_mean_theta(:,1), M.MPF_apost_lb_theta(:,1), M.MPF_apost_ub_theta(:,1));
title(' Parameters C(1,1)')
subplot(4,1,2);
plotestimates(fy.A(1,2)*ones(ndat,1), M.MPF_apost_mean_theta(:,2), M.MPF_apost_lb_theta(:,2), M.MPF_apost_ub_theta(:,2));
title(' Parameters C(1,2)')
subplot(4,1,3);
plotestimates(fy.A(2,1)*ones(ndat,1), M.MPF_apost_mean_theta(:,3), M.MPF_apost_lb_theta(:,3), M.MPF_apost_ub_theta(:,3));
title(' Parameters C(2,1)')
subplot(4,1,4);
plotestimates(fy.A(2,2)*ones(ndat,1), M.MPF_apost_mean_theta(:,4), M.MPF_apost_lb_theta(:,4), M.MPF_apost_ub_theta(:,4));
title(' Parameters C(2,2)')