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.03 -0.02; -0.02 0.05];
fx.g = g;
fx.rv = x;
fx.rvc = RVtimes(x,-1);

fy.class = 'mgnorm<chmat>';
fy.R = 0.01*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.rv = RV('R',16);
A.rgr = RV({});
A.dimx=4;
A.constant = 0;
A.frg=1.0;
A.prior.class='egiw';
A.prior.dimx=4;
A.prior.V =1e-3*eye(4);
%A.prior.nu = 10;


M.class = 'NoiseParticle';
M.g = g;
M.h = h;
M.yrv = y;
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 = 0.9;
PF.prior.class = 'enorm<ldmat>';
PF.prior.mu = [0.2;0.3];
PF.prior.R = 0.1*eye(2);


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

% plot

figure(1);
hold off
plot(O.Est0_apost_mean_R);

figure(2)
hold off
plot(O.Est0_apost_mean_x);
hold on
plot(O.DS_dt_x,'--');