root/doc/latex/manual.tex @ 275

\index{User Manual@{User Manual}}

This page contains some tutorial examples that will help you getting started using BDM.

Basic linear algebra in BDM (using IT++):\begin{itemize}
\item \hyperlink{vector_and_matrix}{A very simple tutorial about vectors and matrixes}\item \hyperlink{itfile}{Writing and reading data from files}\item \hyperlink{timer}{Using timers to measure execution time}\end{itemize}


Basic concepts and philosophy of BDM:\begin{itemize}
\item \hyperlink{intro}{Introduction to Bayesian Decision Making Toolbox BDM}\item \hyperlink{ui}{User Infos and their use}\item \hyperlink{mexfiles}{How to write and use mex files for Matlab}\end{itemize}


Library of predefined tasks:

Examples for Various mathematical models (BDM):\begin{itemize}
\item \hyperlink{arx_ui}{Running experiment {\tt estimator} with ARX data fields}\item \hyperlink{kalman}{Examples of (extended) Kalman filtering}\end{itemize}


Examples of control (BDM):\begin{itemize}
\item To be done... \end{itemize}
\hypertarget{vector_and_matrix}{}\section{A very simple tutorial about vectors and matrixes}\label{vector_and_matrix}
Let's start with a really simple example. Try to complile the following program:



\begin{DocInclude}\begin{verbatim}#include <itpp/itbase.h>

using namespace itpp;

//These lines are needed for use of cout and endl
using std::cout;
using std::endl;

int main()
{
  //Declare vectors and matricies:
  vec a, b, c;
  mat A, B;

  //Use the function linspace to define a vector:
  a = linspace(1.0, 2.0, 10);

  //Use a string of values to define a vector:
  b = "0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0";

  //Add two vectors:
  c = a + b;

  //Print results:
  cout << "a = " << a << endl;
  cout << "b = " << b << endl;
  cout << "c = " << c << endl;

  //Use a string to define a matrix:
  A = "1.0 2.0;3.0 4.0";

  //Calculate the inverse of matrix A:
  B = inv(A);

  //Print results:
  cout << "A = " << A << endl;
  cout << "B = " << B << endl;

  //Exit program:
  return 0;

}
\end{verbatim}
\end{DocInclude}


When you run this program, the output shall look like this



\begin{DocInclude}\begin{verbatim}a = [1 1.11111 1.22222 1.33333 1.44444 1.55556 1.66667 1.77778 1.88889 2]
b = [0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1]
c = [1.1 1.31111 1.52222 1.73333 1.94444 2.15556 2.36667 2.57778 2.78889 3]
A = [[1 2]
 [3 4]]
B = [[-2 1]
 [1.5 -0.5]]
\end{verbatim}
\end{DocInclude}


If this is what you see, then congratulations! You have managed to compile your first it++ program! \hypertarget{itfile}{}\section{Writing and reading data from files}\label{itfile}
Here we will use the {\tt it\_\-file} class to store some data. The program {\tt write\_\-it\_\-file.cpp} looks as follows:



\begin{DocInclude}\begin{verbatim}#include <itpp/itcomm.h>

using namespace itpp;

int main()
{
  // Declare the it_file class
  it_file ff;

  // Open a file with the name "it_file_test.it"
  ff.open("it_file_test.it");

  // Create some data to put into the file
  vec a = linspace(1, 20, 20);

  // Put the variable a into the file. The Name("a") tells the file class
  // that the next variable shall be named "a".
  ff << Name("a") << a;

  // Force the file to be written to disc. This is useful when performing
  // iterations and ensures that the information is not stored in any cache
  // memory. In this simple example it is not necessary to flush the file.
  ff.flush();

  // Close the file
  ff.close();

  // Exit program
  return 0;
}
\end{verbatim}
\end{DocInclude}


When you run this program you will obtain a file called {\tt it\_\-file\_\-test.it} in your current directory. You can read the file into Matlab/Octave to view the data by using the following commands:



\begin{Code}\begin{verbatim}itload('it_file_test.it')
figure(1); clf;
plot(a)
\end{verbatim}
\end{Code}



Note: Make sure that {\tt \$PREFIX/share/itpp} is in your Matlab/Octave path and that you run the code above from the directory where {\tt it\_\-file\_\-test.it} is located ({\tt \$PREFIX} is the IT++ installation prefix; {\tt /usr/local} by default).

The IT++ program {\tt read\_\-it\_\-file.cpp} that reads the file and prints its content can look like this:



\begin{DocInclude}\begin{verbatim}\end{verbatim}
\end{DocInclude}


Here is the output of the program:



\footnotesize\begin{verbatim}
a = [1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20]
\end{verbatim}
\normalsize
 \hypertarget{timer}{}\section{Using timers to measure execution time}\label{timer}
In this example we are using the Real\_\-Timer class to measure the execution time of a simple program. The Real\_\-Timer class is included in the itmisc library.



\begin{DocInclude}\begin{verbatim}#include <itpp/itbase.h>

using namespace itpp;

//These lines are needed for use of cout and endl
using std::cout;
using std::endl;

int main()
{
  //Declare the scalars used:
  long i, sum, N;

  //Declare tt as an instance of the timer class:
  Real_Timer tt;

  //Initiate the variables:
  N = 1000000;
  sum = 0;

  //Start and reset the timer:
  tt.tic();

  //Do some processing
  for (i = 0; i < N; i++) {
    sum += i;
  }

  // Print the elapsed time
  tt.toc_print();

  //Print the result of the processing:
  cout << "The sum of all integers from 0 to " << N - 1 << " equals " << sum << endl;

  //Exit program:
  return 0;

}
\end{verbatim}
\end{DocInclude}


When you run this program, the output will look something like this:



\begin{DocInclude}\begin{verbatim}Elapsed time = 0.000797055 seconds
The sum of all integers from 0 to 999999 equals 1783293664
\end{verbatim}
\end{DocInclude}
 \hypertarget{intro}{}\section{Introduction to Bayesian Decision Making Toolbox BDM}\label{intro}
This is a brief introduction into elements used in the BDM. The toolbox was designed for two principle tasks:

\begin{itemize}
\item Design of Bayesian decisions-making startegies,  \item Bayesian system identification for on-line and off-line scenarios.  \end{itemize}
Theoretically, the latter is a special case of the former, however we list it separately to highlight its importance in practical applications.

Here, we describe basic objects that are required for implementation of the Bayesian parameter estimation.

Key objects are: \begin{description}
\item[Bayesian Model: class {\tt BM}  ]which is an encapsulation of the likelihood function, the prior and methodology of evaluation of the Bayes rule. This methodology may be either exact or approximate. \item[Posterior density of the parameter: class {\tt epdf}  ]representing posterior density of the parameter. Methods defined on this class allow any manipulation of the posterior, such as moment evaluation, marginalization and conditioning.  \end{description}
\hypertarget{intro_bm}{}\section{Class BM}\label{intro_bm}
The class BM is designed for both on-line and off-line estimation. We make the following assumptions about data: \begin{itemize}
\item an individual data record is stored in a vector, {\tt vec} {\tt dt},  \item a set of data records is stored in a matrix,{\tt mat} {\tt D}, where each column represent one individual data record  \end{itemize}


On-line estimation is implemented by method 

\begin{Code}\begin{verbatim} void bayes(vec dt)
\end{verbatim}
\end{Code}

 Off-line estimation is implemented by method 

\begin{Code}\begin{verbatim} void bayesB(mat D)
\end{verbatim}
\end{Code}



As an intermediate product, the bayes rule computes marginal likelihood of the data records $ f(D) $. Numerical value of this quantity which is important e.g. for model selection can be obtained by calling method {\tt \_\-ll()}.\hypertarget{intro_epdf}{}\section{Getting results from BM}\label{intro_epdf}
Class {\tt BM} offers several ways how to obtain results: \begin{itemize}
\item generation of posterior or predictive pdfs, methods {\tt \_\-epdf()} and {\tt predictor()}  \item direct evaluation of predictive likelihood, method {\tt logpred()}  \end{itemize}
Underscore in the name of method {\tt \_\-epdf()} indicate that the method returns a pointer to the internal posterior density of the model. On the other hand, {\tt predictor} creates a new structure of type {\tt epdf()}.

Direct evaluation of predictive pdfs via logpred offers a shortcut for more efficient implementation.\hypertarget{intro_epdf}{}\section{Getting results from BM}\label{intro_epdf}
As introduced above, the results of parameter estimation are in the form of probability density function conditioned on numerical values. This type of information is represented by class {\tt epdf}.

This class allows such as moment evaluation via methods {\tt mean()} and {\tt variance()}, marginalization via method {\tt marginal()}, and conditioning via method {\tt condition()}.

Also, it allows generation of a sample via {\tt sample()} and evaluation of one value of the posterior parameter likelihood via {\tt evallog()}. Multivariate versions of these operations are also available by adding suffix {\tt \_\-m}, i.e. {\tt sample\_\-m()} and {\tt evallog\_\-m()}. These methods providen multiple samples and evaluation of likelihood in multiple points respectively.\hypertarget{intro_pc}{}\section{Classes for probability calculus}\label{intro_pc}
When a more demanding task then generation of point estimate of the parameter is required, the power of general probability claculus can be used. The following classes (together with {\tt epdf} introduced above) form the basis of the calculus: \begin{itemize}
\item {\tt mpdf} a pdf conditioned on another symbolic variable, \item {\tt RV} a symbolic variable on which pdfs are defined. \end{itemize}
The former class is an extension of mpdf that allows conditioning on a symbolic variable. Hence, when numerical results - such as samples - are required, numericla values of the condition must be provided. The names of methods of the {\tt epdf} are used extended by suffix {\tt cond}, i.e. {\tt samplecond()}, {\tt evallogcond()}, where {\tt cond} precedes matrix estension, i.e. {\tt samplecond\_\-m()} and {\tt evallogcond\_\-m()}.

The latter class is used to identify how symbolic variables are to be combined together. For example, consider the task of composition of pdfs via the chain rule: \[ f(a,b,c) = f(a|b,c) f(b) f(c) \] In our setup, $ f(a|b,c) $ is represented by an {\tt mpdf} while $ f(b) $ and $ f(c) $ by two {\tt epdfs}. We need to distinguish the latter two from each other and to deside in which order they should be added to the mpdf. This distinction is facilitated by the class {\tt RV} which uniquely identify a random varibale.

Therefore, each pdf keeps record on which RVs it represents; {\tt epdf} needs to know only one {\tt RV} stored in the attribute {\tt rv}; {\tt mpdf} needs to keep two {\tt RVs}, one for variable on which it is defined ({\tt rv}) and one for variable incondition which is stored in attribute {\tt rvc}. \hypertarget{ui}{}\section{User Infos and their use}\label{ui}
For easier interaction with users, experiments can be configured via structures called User Info. These structures contain information about details of the experiment that is to be performed. Since experiments involve the use of basic BDM classes and their compositions, the experiment description is also hierarchical with specific user info for each object or class.

The User Infos are designed using customized version of the libconfig library (link!). It is specialized for BDM so as to recognize basic mathematical objects, such as vectors and matrices, see ... for details.

(Technically it can be made compatible with matlab structures!)

For example a simple experiment can be configures in a following way: 

\begin{Code}\begin{verbatim}ndat = 100;                //number of data points
prior = {type="enorm";
        mu = [1, 2, 3];
        R = [1, 0, 0,
             0, 1, 0,
             0, 0, 1];
};
\end{verbatim}
\end{Code}

 Exact meaning of root fields in this structure (i.e. ndat and prior) is defined by the application (or mex file) that is using this configuration file. It will look for expected fields and it will ignore any other structures. When it does not find what it is looking for, it terminates with an appropriate error message.

A structure with field {\tt type=\char`\"{}identifier\char`\"{}} is special. Such a structure will be parsed by an appropriate class \hyperlink{classbdm_1_1UIbuilder}{bdm::UIbuilder} which will construct the desired object, in this instance of an object of the class \hyperlink{classbdm_1_1enorm}{bdm::enorm}. For a detailed example how this mechanism works in practice see \hyperlink{arx_ui}{Running experiment {\tt estimator} with ARX data fields}. \hypertarget{mexfiles}{}\section{How to write and use mex files for Matlab}\label{mexfiles}
\hypertarget{mexfiles_use}{}\section{Howto use predefined mexfiles}\label{mexfiles_use}
A range of mexfiles is predefined in directory {\tt library/mex}. Many of these mexfile process ui files (see \hyperlink{ui}{User Infos and their use}) examples of these files are in directory {\tt library/tutorial}. Note that in order to run these files you need to let matlab know where to find them: 

\begin{Code}\begin{verbatim}>> addpath path_to_bdm/library/mex
\end{verbatim}
\end{Code}



Then, you can go to {\tt library/tutorial} and run e.g. {\tt arx\_\-test\_\-mex}.\hypertarget{mexfiles_write}{}\section{Howto write custom mex file}\label{mexfiles_write}
Due to special nature of the mex files, the mex file can be split in three parts: \begin{itemize}
\item input conversion, where input arguments are converted to IT++ structures, \item main body of algorithm, where any C++ bdm constructions can be used, \item output conversion, where resulting IT++ structures are converted to mxArrays.\end{itemize}
The first and the third part is achieved using prepared IT++ routines, see IT++ documentation.

Script {\tt }./buildmex is prepared to compile and link the mexfile with bdm. 

\begin{Code}\begin{verbatim}$ ./buildmex.sh my_mex_file.cpp
\end{verbatim}
\end{Code}

 on Linux, or 

\begin{Code}\begin{verbatim}$ ./buildmex.bat my_mex_file.cpp
\end{verbatim}
\end{Code}

 on Windows.

Example of a mexfile: 

\begin{DocInclude}\begin{verbatim}#include <itpp/itmex.h>
#include <estim/arx.h>

using namespace bdm;

void mexFunction(int n_output, mxArray *output[], int n_input, const mxArray *input[])
{
    // Check the number of inputs and output arguments
        if(n_output!=1) mexErrMsgTxt("Wrong number of output variables!");
        if(n_input!=2) mexErrMsgTxt("Usage: arx1d(ysize, Data)!");

    // Convert input variables to IT++ format
        int ysize = mxArray2int(input[0]);
        mat Data = mxArray2mat(input[1]);

    // ------------------ Start of routine ---------------------------
        ARX Ar;
        Ar.set_statistics(ysize, 1e-5*eye(Data.rows()) );
        Ar.bayesB(Data);
        // ------------------ End of routine -----------------------------

    // Create output vectors
        output[0] = mxCreateDoubleMatrix(1,Data.rows(), mxREAL);

    // Convert the IT++ format to Matlab format for output
        vec2mxArray(Ar.posterior().mean(), output[0]);
}
\end{verbatim}
\end{DocInclude}
 \hypertarget{arx_ui}{}\section{Running experiment $\backslash$c estimator with ARX data fields}\label{arx_ui}
The experiment estimator::cpp can be run either on command line, or as a mex file in Matlab.\hypertarget{arx_ui_cmd}{}\section{Command-line usage}\label{arx_ui_cmd}
In order to use it for estimation of an ARX model, we can define the following \hyperlink{ui}{User Infos and their use} structure: 

\begin{DocInclude}\begin{verbatim}//Data generating system
system = {
        type = "ArxDS";
        y = {type="rv"; names=["y", "u"];};
        u = {type="rv"; names=[]; };
        rgr = {type="rv";
                names = ["y","y","y","u"];
                times = [-1, -2, -3, -1];
        };
        //AR parameters
        theta = [0.8, -0.3, 0.4, 1.0,
                 0.0, 0.0, 0.0, 0.0];
        // offset
        offset = [0.0, 0.0];
        //variance
        r = [0.1, 0.0,
             0.0, 1.0];
        // log also theta
        opt="L_theta";
};

//store results
logger = {
        type= "dirfilelog";
        dirname = "exp/arx_ui";
        maxlen = 1000; //
};

//estimation
estimator = {
        type = "ARXest";
        y = {type="rv"; names=["y"]; };
        rgr = {type="rv";
                names = ["y","y","y","u"];
                times = [-1, -2, -3, -1];
        };

        //optional fields
        dV0 = [1e-3, 1e-5, 1e-5, 1e-5, 1e-5]; //default: 1e-3 for y, 1e-5 for rgr
        nu0 = 8.;      //default: rgrlen + 2
        frg = 1.0;    // forgetting, default frg=1.0
};

//experiment description
experiment:{
        ndat = 9000;
};
\end{verbatim}
\end{DocInclude}


The structure is interpreted by application {\tt estimator}, which looks for fields: \begin{description}
\item[system]description of a Data Source generating Data. The structure must by UI with {\tt type=\char`\"{}DS\_\-offspring\char`\"{}}. In our example, it is of type \char`\"{}ArxDS\char`\"{} which is parsed by bdm::UIArxDS UIbuilder generating a Data Source simulating ARX process. \item[estimator ]description of a Baysian model used to estimate parameters of the data model. In this case, it is of type \char`\"{}ARXest\char`\"{} which is parsed by bdm::UIARX UIbuilder generating Bayesian estimator of autoregressive processess. \item[logger]description of a way how to store results. UI is of {\tt type=\char`\"{}logger\_\-offspring\char`\"{}}. In this case, it is of class \char`\"{}dirfilelog\char`\"{} which is parsed by bdm::UIdirfilelog which generates object storing data in directory specified by dirname=\char`\"{}\char`\"{} field in fileformat understood by program kst. \end{description}
When the application estimator is run with the above code, it produces a directory of data files, which can be displayed by program kst. Expected results are:  \begin{Image}
\begin{center}
\includegraphics[width=\linewidth]{arx_ui_kst.png}\caption{Typical run of estimator arx\_\-test.cfg}
\end{center}
\end{Image}
\hypertarget{arx_ui_mex}{}\section{Matlab mex file}\label{arx_ui_mex}
The matlab mex file can be run with exactly the same configuration as above. However, when we wish to see the results in Matlab, we may wish to change the logger object to {\tt type=\char`\"{}mexlog\char`\"{}} which will store the results in a matlab structure.

The exact configuration file may look as follows: 

\begin{DocInclude}\begin{verbatim}//Data generating system
system = { type="external"; filename="arx_test.cfg"; path="system";};

//store results
logger = {
        type= "mexlog";
        maxlen = 90; //
    //dirname = "exp/ax";
};

//estimation
estimator = {
   type = "ARXest";
        y = {type="rv"; names=["y"]; };
        rgr = {type="rv";
                names = ["y","y","y","u"];
                times = [-1, -2, -3, -1];
        };

        //optional fields
        dV0 = [1e-3, 1e-5, 1e-5, 1e-5, 1e-5]; //default: 1e-3 for y, 1e-5 for rgr
        //nu0 = 8.;      //default: rgrlen + 2
        frg = .9991;    // forgetting, default frg=1.0
};

//experiment description
experiment:{
        ndat = 90;
};
\end{verbatim}
\end{DocInclude}


The resulting structure can be displayed using matlab script arx\_\-test\_\-disp.m, typically producing the following results:  \begin{Image}
\begin{center}
\includegraphics[width=\linewidth]{arx_ui_mex.png}\caption{Typical run of matlab file arx\_\-test\_\-disp}
\end{center}
\end{Image}
 \hypertarget{kalman}{}\section{Examples of (extended) Kalman filtering}\label{kalman}
Kalman filtering and Extended Kalman filtering are special cases of Bayesian filtering. The Kalman filter is optimal for linear state space model with Gaussian disturbances, the extended Kalman filter is derived as linearization of non-linear state space models with Gaussian noises. Hence it is only sub-optimal filter.

More advanced filtering algorithms for non-linear non-Gaussian models can be derived, see ...\hypertarget{kalman_klm}{}\section{Kalman Filtering}\label{kalman_klm}
Kalman filtering is optimal estimation procedure for linear state space model: \begin{eqnarray} x_t &= &A x_{t-1} + B u_{t} + v_t,\\ y_t &= &C x_{t} + D u_{t} + w_t, \end{eqnarray} where $ x_t $ is the state, $ y_t $ is the system output, $ A, B, C, D$ are state matrices of appropriate dimensions, $v_t, w_t$ are zero mean Gaussian noises with covariance matrices $Q, R$, respectively.

Both prior and posterior densities on the state are Gaussian, i.e. of the class enorm.

There is a range of classes that implements this functionality, namely:\begin{itemize}
\item \hyperlink{classbdm_1_1KalmanFull}{bdm::KalmanFull} which implements the estimation algorithm on full matrices,\item \hyperlink{classbdm_1_1KalmanCh}{bdm::KalmanCh} which implements the estimation algorithm using choleski decompositions and QR algorithm.\end{itemize}
\hypertarget{kalman_ekf}{}\section{Extended Kalman Filtering}\label{kalman_ekf}
Extended Kalman filtering arise by linearization of non-linear state space model: \begin{eqnarray} x_t &= &g( x_{t-1}, u_{t}) + v_t,\\ y_t &= &h( x_{t} , u_{t}) + w_t, \end{eqnarray} where $ g(), h() $ are general non-linear functions which have finite derivatives. Remaining variables have the same meaning as in the Kalman Filter.

In order to use this class, the non-linear functions and their derivatives must be defined as an instance of class {\tt diffbifn}.

Two classes are defined:\begin{itemize}
\item \hyperlink{classbdm_1_1EKFfull}{bdm::EKFfull} on full size matrices,\item \hyperlink{classbdm_1_1EKFCh}{bdm::EKFCh} on Choleski decompositions and using QR algorithm.\end{itemize}
\hypertarget{kalman_exa}{}\section{Examples of Use}\label{kalman_exa}
The classes can be used directly in C++ or via User Info. The latter example is illustrated in file estimator. A very short example of the former follows:



\begin{DocInclude}\begin{verbatim}#include <estim/libKF.h>
using namespace bdm;
        
// estimation of AR(0) model
int main() {
        //dimensions
        int dx=3, dy=3, du=1;
        // matrices
        mat A = eye(dx);
        mat B = zeros(dx,du);
        mat C = eye(dx);
        mat D = zeros(dy,du);
        mat Q = eye(dx);
        mat R = 0.1*eye(dy);
        //prior
        mat P0 = 100*eye(dx);
        vec mu0 = zeros(dx);
        // Estimator
        KalmanCh KF;
        KF.set_parameters(A,B,C,D,/*covariances*/ Q,R);
        KF.set_statistics(mu0,P0);
        // Estimation loop
        for (int i=0;i<100;i++){
                KF.bayes(randn(dx+du));
        }
        //print results
        cout << "Posterior estimate of x is: "  << endl;
        cout << "mean: "<< KF.posterior().mean()<< endl;
        cout << "variance: "<< KF.posterior().variance()<< endl;
}
\end{verbatim}
\end{DocInclude}