root/applications/bdmtoolbox/tutorial/userguide/partfrg_estim.m @ 1100

clear all
% load data created by the MpdfDS_example
load pdfds_results

DS.class   = 'MemDS';
DS.Data    = Data;
DS.drv     = drv;


%%%%%% ARX estimator conditioned on frg

A1.class = 'ARXpartialforg';
A1.yrv = y;
A1.rv = RV({'theta','r'},[2,1]);
A1.rgr = RVtimes([y,u],[-3,-1]) ; 
A1.log_level = 'logbounds';
A1.constant = 0;
A1.name = 'A1';

Apri =A1;
Apri.class = 'ARX';
DSpri = DS;
DSpri.Data = DS.Data(:,1:6);
% get decent prior -- estimate with data first
[Dum,post]=estimator(DSpri,{Apri});

A1.prior = post.estimators{1}.posterior;

% we have 2 parameters - i.e. 4 hypotheses

%%%%%% Random walk on frg - Dirichlet 
walk.class = 'mDirich';         % random walk on coefficient phi
walk.rv    = RV({'phi'},4);       % 2D random walk - frg is the first element
walk.k     = 0.001;              % width of the random walk
walk.betac = 0.1*ones(1,4);         % stabilizing elememnt of random walk

%%%%%% Particle
p.class = 'MarginalizedParticle';
p.parameter_pdf = walk;         % Random walk is the parameter evolution model
p.bm    = A1;

% prior on ARX
%%%%%% Combining estimators in Marginalized particle filter
E.class = 'PF';
E.particle = p;                    % ARX is the analytical part
E.res_threshold = 0.7;             % resampling parameter
E.n = 100;                           % number of particles
E.prior.class = 'eBeta';         % prior on non-linear part
E.prior.alpha = 5*ones(1,4); % 
E.prior.beta  = ones(1,4); % 
E.log_level = 'logbounds,logweights';
E.name = 'MPF';

[M,Str]=estimator(DS,{E});

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

figure(1);
subplot(2,2,1);
plotestimates(true_theta, M.MPF_apost_mean_theta, M.MPF_apost_lbound_theta, M.MPF_apost_ubound_theta);
title(' Regression parameters \theta')
set(gca,'YLim',[-1.5,1]);

subplot(2,2,2);
plotestimates(true_R, M.MPF_apost_mean_r,M.MPF_apost_lbound_r,M.MPF_apost_ubound_r);
title('Variance parameters r')

subplot(2,2,3);
plotestimates(1, M.MPF_apost_mean_phi(:,1),M.MPF_apost_lbound_phi(:,1),M.MPF_apost_ubound_phi(:,1));
title('Forgetting factor')