root/library/doc/tutorial/02userguide_estim.dox @ 1321

/*!
\page userguide_estim BDM Use - Estimation and Bayes Rule

Bayesian theory is predominantly used in system identification, or estimation problems. 
This section is concerned with recursive estimation, as implemented in prepared scenario \c estimator.

Table of contents:
\ref ug2_theory 
\ref ug2_arx_basic 
\ref ug2_model_sel 
\ref ug2_bm_composition 
\ref ug_est_ext 

The function of the \c estimator is graphically illustrated:
\dot
digraph estimation{
        node [shape=box];
        {rank="same"; "Data Source"; "Bayesian Model"}
        "Data Source" -> "Bayesian Model" [label="data"];
        "Bayesian Model" -> "Result Logger" [label="estimation\n result"];
        "Data Source" -> "Result Logger" [label="Simulated\n data"];
}
\enddot

Here,
\li Data Source is an object (class DS) providing sequential data, \f$ [d_1, d_2, \ldots d_t] \f$.
\li Bayesian Model is an object (class BM) performing Bayesian filtering,
\li Result Logger is an object (class logger) dedicated to storing important data from the experiment.

Since objects  datasource and the logger has already been introduced in section \ref userguide_sim, it remains to introduce 
object \c Bayesian \c Model (bdm::BM).

\section ug2_theory Bayes rule and estimation
The object bdm::BM is basic software image of the Bayes rule:
\f[ f(x_t|d_1\ldots d_t) \propto f(d_t|x_t,d_1\ldots d_{t-1}) f(x_t| d_1\ldots d_{t-1}) \f]

Since this operation can not be defined universally, the object is defined as abstract class with methods for:
 - <b> Bayes rule </b> as defined above, operation bdm::BM::bayes() which expects to get the current data record \c dt, \f$ d_t \f$
 - <b> evidence </b> i.e. numerical value of \f$ f(d_t|d_1\ldots d_{t-1})\f$ as a typical side-product, since it is required in denominator of the above formula.
   For some models, computation of this value may require extra effort, and can be switched off.
 - <b> prediction </b> the object has enough information to create the one-step ahead predictor, i.e. \f[ f(d_{t+1}| d_1 \ldots d_{t}), \f] 
        
Implementation of these operations is heavily dependent on the specific class of prior pdf, or its approximations. We can identify only a few principal approaches to this problem. For example, analytical estimation which is possible within sufficient the Exponential Family, or estimation when both prior and posterior are approximated by empirical densities. 
These approaches are first level of descendants of class \c BM, classes bdm::BMEF and bdm::PF, respectively.

\section ug2_arx_basic Estimation of ARX models

Autoregressive models has already been introduced in \ref ug_arx_sim where their simulator has been presented.
We will use results of simulation of the ARX datasource defined there to provide data for estimation using MemDS.

The following code is from bdmtoolbox/tutorial/userguide/arx_basic_example.m
\code
A1.class = 'ARX';
A1.rv = y;
A1.rgr = RVtimes([y,y],[-3,-1]) ;
A1.log_level = 'logbounds,logevidence';
\endcode 
This is the minimal configuration of an ARX estimator. 

The first three fields are self explanatory, they identify which data are predicted (field \c rv) and which are in regresor (field \c rgr).
The field \c log_level is a string of options passed to the object. In particular, class \c BM understand only options related to storing results:
 - logbounds - store also lower and upper bounds on estimates (obtained by calling BM::posterior().qbounds()),
 - logevidence - store also evidence of each step of the Bayes rule.
These values are stored in given logger (\ref ug_store). By default, only mean values of the estimate are stored.

Storing of the evidence is useful, e.g. in model selection task when two models are compared.

The bounds are useful e.g. for visualization of the results. Run of the example should provide result like the following:
\image latex arx_basic_example.png "Typical run of tutorial/userguide/arx_basic_example.m" width=\linewidth

\section ug2_model_sel Model selection

In Bayesian framework, model selection is done via comparison of evidence (marginal likelihood) of the recorded data. See [some theory].

A trivial example how this can be done is presented in file bdmtoolbox/tutorial/userguide/arx_selection_example.m. The code extends the basic A1 object as follows:
\code
A2=A1;
A2.constant = 0;

A3=A2;
A3.frg = 0.95;
\endcode
That is, two other ARX estimators are created, 
 - A2 which is the same as A1 except it does not model constant term in the linear regression. Note that if the constant was set to zero, then this is the correct model.
 - A3 which is the same as A2, but assumes time-variant parameters with forgetting factor 0.95.
 
Since all estimator were configured to store values of marginal log-likelihood, we can easily compare them by computing total log-likelihood for each of them and converting them to probabilities. Typically, the results should look like:
\code
Model_probabilities =

    0.0002    0.7318    0.2680
\endcode
Hence, the true model A2 was correctly identified as the most likely to produce this data. 

For this task, additional technical adjustments were needed:
\code
A1.name='A1';
A2.name='A2';
A2.rv_param = RV({'a2th', 'r'},[2,1],[0,0]);
A3.name='A3';
A3.rv_param = RV({'a3th', 'r'},[2,1],[0,0]);
\endcode
First, in order to distinguish the estimators from each other, the estimators were given names. Hence, the results will be logged with prefix given by the name, such as M.A1_evidence.

Second, if the parameters of a ARX model are not specified, they are automatically named \c theta and \c r. However, in this case, \c A1 and \c A2 differ in size, hence their random variables differ and can not use the same name. Therefore, we have explicitly used another names (RVs) of the parameters.

\section ug2_bm_composition Composition of estimators

Similarly to pdfs which could be composed via \c mprod, the Bayesian models can be composed together. However, justification of this step is less clear than in the case of epdfs.

One possible theoretical base of composition is the Marginalized particle filter, which splits the prior and the posterior in two parts:
\f[ f(x_t|d_1\ldots d_t)=f(x_{1,t}|x_{2,t},d_1\ldots d_t)f(x_{2,t}|d_1\ldots d_t) \f]
each of these parts is estimated using different approach. The first part is assumed to be analytically tractable, while the second is approximated using empirical approximation.

The whole algorithm runs by parallel evaluation of many \c BMs for estimation of \f$ x_{1,t}\f$, each of them conditioned on value of a sample of \f$ x_{2,t}\f$.

For example, the forgetting factor, \f$ \phi \f$ of an ARX model can be considered to be unknown. Then, the whole parameter space is \f$ [\theta_t, r_t, \phi_t]\f$ decomposed as follows:
\f[ f(\theta_t, r_t, \phi_t) = f(\theta_t, r_t| \phi_t) f(\phi_t) \f]
Note that for known trajectory of \f$ \phi_t \f$ the standard ARX estimator can be used if we find a way how to feed the changing \f$ \phi_t \f$ into it.
This is achieved by a trivial extension using inheritance method bdm::BM::condition().

Extension of standard ARX estimator to conditional estimator is implemented as class bdm::ARXfrg. The only difference from standard ARX is that this object will obtain its forgetting factor externally as a conditioning variable. 
Informally, the name 'ARXfrg' means: "if anybody calls your condition(0.9), it tells you new value of forgetting factor".

The MPF estimator for this case is specified as follows:
\code
%%%%%% ARX estimator conditioned on frg

A1.class = 'ARXfrg';
A1.rv = y;
A1.rgr = RVtimes([y,u],[-3,-1]) ; 
A1.log_level ='logbounds,logevidence';
A1.frg = 0.9;
A1.name = 'A1';

%%%%%% Random walk on frg - Dirichlet 
phi_pdf.class = 'mDirich';         % random walk on coefficient phi
phi_pdf.rv    = RV('phi',2);       % 2D random walk - frg is the first element
phi_pdf.k     = 0.01;              % width of the random walk
phi_pdf.betac = [0.01 0.01];       % stabilizing elememnt of random walk

%%%%%% Particle
p.class = 'MarginalizedParticle';
p.parameter_pdf = phi_pdf;         % 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 = 1.0;             % resampling parameter
E.n = 100;                         % number of particles
E.prior.class = 'eDirich';         % prior on non-linear part
E.prior.beta  = [2 1]; % 
E.log_level = 'logbounds';
E.name = 'MPF';

M=estimator(DS,{E});

\endcode
 
Here, the configuration structure \c A1 is a description of an ARX model, as used in previous examples, the only difference is in its name 'ARXfrg'. 

The configuration structure \c phi_pdf defines random walk on the forgetting factor. It was chosen as Dirichlet, hence it will produce 2-dimensional vector of \f$[\phi, 1-\phi]\f$. The class \c ARXfrg was designed to read only the first element of its condition.
The random walk of type mDirich is:
\f[ f(\phi_t|\phi_{t-1}) = Di (\phi_{t-1}/k + \beta_c) \f]
where \f$ k \f$ influences the spread of the walk and \f$ \beta_c \f$ has the role of stabilizing, to avoid traps of corner cases such as [0,1] and [1,0].
Its influence on the results is quite important.

This example is implemented as bdmtoolbox/tutorial/userguide/frg_example.m
Its typical run should look like the following:
\image html frg_example_small.png 
\image latex frg_example.png "Typical run of tutorial/userguide/frg_example.m" width=\linewidth

Note: error bars in this case are not directly comparable with those of previous examples. The MPF class implements the qbounds function as minimum and maximum of bounds in the considered set (even if its weight is extremely small). Hence, the bounds of the MPF are probably larger than it should be. Nevertheless, they provide great help when designing and tuning algorithms.

\section ug_est_ext Matlab extensions of the Bayesian estimators

Similarly to the extension of pdf, the estimators (or filters) can be extended via prepared class \c mexBM in directory bdmtoolbox/mex/mex_classes.

An example of such class is mexLaplaceBM in \<toolbox_dir\>/tutorial/userguide/laplace_example.m

Note that matlab-extended classes of mexEpdf, specifically, mexDirac and mexLaplace are used as outputs of methods posterior and epredictor, respectively.

In order to create a new extension of an estimator, copy file with class mexLaplaceBM.m and redefine the methods therein. If needed create new classes for pdfs by inheriting from mexEpdf, it the same way as in the mexLaplace.m example class.

For list of all estimators, see \ref app_base.
*/