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'; 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'; 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')