clear all;
x = RV('x',2);
y = RV('y',2);

g.class = 'mexFnc';
g.dim = 2;
g.dimc = 2;
g.function = 'test_function';

g.class = 'linfn';
g.A     = eye(2);
g.B     = [1;0];

h.class = 'linfn';
h.A     = eye(2);
h.B     = [1;0];

fx.class = 'mgnorm<chmat>';
fx.R = [0.3 -0.2; -0.2 0.5];
fx.g = g;
fx.rv = x;
fx.rvc = RVtimes(x,-1);

fy.class = 'mgnorm<chmat>';
fy.R = 0.1*eye(2);
fy.g = h;
fy.rv = y;
fy.rvc = x;


%%%%%% Data generator
DS.class   = 'PdfDS';
DS.pdf.class     = 'mprod';
DS.pdf.pdfs      = {fy,fx};
DS.init_rv = RVtimes(x,-1);
DS.init_values = [.2,.3]';

%%%%% Estimator
A.class = 'ARX';
A.yrv = RV('vw',4);
A.rgr = RV({});
A.dimx=4;
A.constant = 0;
A.frg=0.99;

M.class = 'NoiseParticle';
M.g = g;
M.h = h;
M.rvx = x;
M.rvxc = RVtimes(x,-1);
M.rvyc = x;
M.bm = A;

PF.class='PF';
PF.particle = M;
PF.n = 100;
PF.res_threshold = 1.0;
PF.prior.class = 'enorm<ldmat>';
PF.prior.mu = [0.2;0.3];
PF.prior.R = 0.1*eye(2);


exper.ndat = 2000;
O = estimator(DS,{PF},exper);
%%%%%% ARX estimator conditioned on frg

% plot

figure(1);
hold off
plot(M.Est0_Q_Q);
hold on
plot(ones(size(M.Est0_Q_Q,1),1)* DS.pdf.pdfs{2}.R(:)','--');

figure(2)
hold off
plot(M.Est0_apost_mean_);
hold on
plot(M.DS_x,'--');