clear all
% name random variables
x = RV({'x'},1);
y = RV({'y'},1);
u = RV({'u'},1);

% create f(x_t| x_{t-1}, u_{t})
fx.class = 'mlnorm<ldmat>';
fx.rv    = x;
fx.rvc   = RVtimes([x,u], [-1, 0]);
fx.A = [0.5 -0.9];
fx.const = 0;
fx.R = 0.01;

% create f(y_t| y_{t-1}, u_{t-1})
fy.class = 'mlnorm<ldmat>';
fy.rv    = y;
fy.rvc   = RVjoin([x,u]);
fy.A     = [1, 0.1];
fy.const = 0;
fy.R     = 1e-3;

% create f(u_t| )
fu.class = 'egauss';
fu.rv    = u;
fu.mu    = 0;
fu.R     = 1e-1;

% create DS
DS.class = 'PdfDS';
DS.pdf.class  = 'mprod';
DS.pdf.pdfs  = { fy, fx, fu};
DS.init_rv = RVtimes([x], [-1]);
DS.init_values = [0.1];

% debug DS % MMM=simulator(DS);

%%%%%% PF estimator

PF.class = 'PF';
PF.particle.class = 'BootstrapParticle';
PF.particle.parameter_pdf = fx;
PF.particle.observation_pdf = fy;
PF.log_level ='logbounds,logevidence';
PF.prior.class = 'egauss';
PF.prior.mu = 0;
PF.prior.R = 0.2;
PF.res_threshold = 1;
PF.n = 1000;

exper.ndat =100;

M=estimator(DS,{PF},exper);

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

plotestimates(M.DS_dt_x, ...
    M.Est0_apost_mean_x, ...
    M.Est0_apost_lbound_x, ...
    M.Est0_apost_ubound_x);
set(gca,'YLim',[-1.5,1]);

title('Variance parameters r')