root/applications/bdmtoolbox/doc/local/20devguide_matlab.dox @ 1192

Revision 1054, 11.2 kB (checked in by smidl, 14 years ago)

doc

Line 
1/*!
2\page devguide_mat  BDMToolbox Development - your first bdm class in Matlab
3
4bdmtoolbox provides a prodefined points of extensions in directory applications/bdmtoolbox/mex/mex_classes.
5
6The easiest way to extend bdmtoolbox is to use clasess starting with "mex". The simplies class of this type is mexFnc, which just contains name of an arbitrary matlab function to call.
7See \ref ug_pdf_fnc example.
8
9This tutorial assumes you are familiar with basic classes of BDM: pdfs, random variables (\ref userguide_pdf) and estimators (\ref userguide_estim).
10
11Previous part of the tutorial introduced ways how to use existing classes. If no premade class is suitable for your mneeds, you can create your own class in Matlab via predefined interfaces to class epdf (mexEpdf) and BM (mexBM).
12These classes are treated as regular C++ classes in prepared mex files. Creation of a new estimator (offspring of class BM) can be achieved by inheriting from class mexBM.
13
14In practice this means creation of a new file, e.g. my_bm.m, with header:
15\code
16class my_bm < mexBM
17\endcode
18This file must be located in directories known to matlab (must be in its path). This can be any directory if you use addpath command.
19See relevant documentation in Matlab.
20
21For list of existing mex* classes see <a href="annotated.html"> list </a>.
22
23\section dev_mat_pdf Creating your own probability density
24
25Creation of your own density is achieved by inheriting from mexEpdf:
26\code
27class my_pdf < mexEpdf
28\endcode
29in a file called my_pdf.m.
30
31Calling methods from mexEpdf will result in an error that tells you which functionality is missing.
32It is necessary to give them a meaning.
33
34Epdf is a class representing probability density function of a multivariate (vector) variable.
35
36The easiest example of a density is dirac Delta function
37\f[
38f(x_t) = \delta( x_t - \bar{x}_t )
39\f]
40with properties:
41 - \f$ \bar{x}_t \f$ is the statistics of this density.
42 - dimension of the density is dimension of the vector \f$ \bar{x}_t \f$
43 - the mean value of the density is \f$ E(x_t) = \bar{x}_t \f$
44 - the variance value of the density is \f$ var(x_t) = zeros() \f$
45 - all samples from this density are equal to \f$ \bar{x}_t \f$
46 
47 These properties are encoded in software as follows:
48 
49\code
50classdef mexDirac < mexEpdf
51    % Dirac delta probability distribution
52    properties                                                 % DATA structures
53        point                                                  % point of the support, \bar{x}_t
54    end
55    methods
56        function m=mean(obj)                                   % compute mean values of the density
57            m = obj.point;
58        end
59        function obj=validate(obj)                             % check if data structures are consistent
60            % point should be a column
61            if (size(obj.point,2)>1)
62                if (size(obj.point,1)==1) % it is row
63                    obj.point = obj.point';
64                end
65            else
66                error('Point in mexDirac is not a vector');
67            end
68        end
69        function dim=dimension(obj)                            % inform others about your dimensions
70            dim = size(obj.point,1);
71        end
72        function v=variance(obj)                               % compute variance
73            v=zeros(size(obj.point));
74        end
75        function l=evallog(obj,x)                              % return logarithm of your density at point x
76            if obj.point==x
77                l = inf;
78            else
79                l=0;
80            end
81        end
82        function s=sample(obj);                                % return random sample from your density
83            s = obj.point;
84        end
85    end
86end
87\endcode
88
89An additional compulsory function is called \c validate which is called always before use of the class in an algorithm. The task of this function is to detect any inconsistency in the object and give the user an appropriate warning.
90For example, in this case, the function checks if the property point is indeed a vector and not a matrix or other structure.
91
92Once you define the methdos above, this object can be plugged into (almost) an arbitrary algoritm, such as simulator, etc.
93
94For seamless use with other objects, you need to define its random variables - via attribute rv, \ref ug_pdf_create.
95
96For inspiration see example class mexLaplace.m and its use in simulation in tutorial/userguide/epdfds_example.m
97
98\section dev_mat_bm Defining your own Bayesian Estimator (BM)
99If you look into the mexBM.m file you see the basic attributes and default functions. mexBM as such will not work as a valid estimator.
100In order to make it work, you need to re-define at least functions <b>dimensions</b> and <b>bayes</b>, and attribute <b>apost_pdf</b>:
101 
102This function announces sizes of random variables to other objects.
103
104These sizes can be fixed for a specific problem, or derived from size of internal objects. E.g. for a Kalman filter, the function would be:
105\code
106function dims=dimensions(p)
107   %please fill
108   %dims = [size_of_posterior size_of_data size_of_condition]
109   dims = [size(A,1), size(C,1), size(B,2)] %
110end
111\endcode
112Which indicates that dimension of \f$x_t\f$ is the rows of matrix A, size of \f$y_t\f$ is the rows of matrix C and the size of \f$u_t\f$ is given by columns of matrix B.
113 
114The argument apost_pdf should be a valid object of class mexEpdf, see previous section how to achieve it.
115
116The second function to extend is
117\code
118function obj=bayes(obj,dt,cond)
119   % transform old estimate into new estimate
120end
121\endcode
122
123Here, you need to define an update (both time- and data-update) of your Bayesian filter.
124
125Note, that statistics of the posterior are stored in attribute apost_pdf. So, the task of this function is to update fileds in apost_pdf.
126
127Any additional attributes needed for this task shoul be added to properties of the object.
128
129
130\section dev_exa Example task: Laplace estimation
131
132Consider a systems, where data are generated by a Laplace model
133\f[
134f(y_t|a,b) = \frac{1}{2b} \exp \left( -\frac{x-a}{b} \right)
135\f]
136with unknown time-varying parameters \f$ a \f$ and \f$ b \f$.
137
138We decided to estimate parameters of the system using ML estimates of parameters on a sliding window. The estimates are then:
139\f{eqnarray}
140\hat{a}& =& \mathrm{median}([y_{t},y_{t-1},\ldots , y_{t-d}),\\
141\hat{b}& = & \frac{1}{d} \sum_{\tau=t-d}^{t} | y_\tau - \hat{a} |
142\f}
143
144Implementing this system in BDM requires the following steps:
145 - realize that ML estimation is an extreme approximation of Bayes rule with posterior in the form of Dirac delat pdf, this pdf is ready from the previous step.
146 - under this choice, the posterior is a 2-dimensional dirac, which is encoded as:
147 \code
148 apost_pdf = mexDirac;
149 apost_pdf.point = [0;0];
150 \endcode
151 - the class mexDirac as introduced has only one attribute: point, which is a vector of the support. If we wish to distinguish the meaning of these values, we use attribute rv. For example,
152 \code
153 apost_pdf.rv = RV({'a','b'});
154 \endcode
155 - the class representing an estimator of this type, needs to remember internal statistics in the form window of data, which will be stored in properties of the class
156
157 - dimensionality of the estimator is then: 2-dimensional parameters, 1-dimensional observation, 0-dimensional conditioning variables. It is encoded as:
158 \code
159  function dims=dimensions(p)
160      %please fill
161      %dims = [size_of_posterior size_of_data size_of_condition]
162      dims = [2,1,0] %
163  end
164  \endcode
165 - the method \c bayes then copies the obtained data into the window and updates estimates \f$ \hat{a}, \hat{b}\f$
166
167This estimator has been implemented as class mexLaplaceBM as follows:
168 
169This is implemented as follows:
170\code
171classdef mexLaplaceBM < mexBM
172    % Approximate Bayesian estimator of parameters of Laplace distributed observation.
173    % Maximum likelihood approximation of the Bayes rule is used, posterior is in the form of dirac.
174    properties
175        max_window_length = 10;                               % max window length (default = 10)
176        data_window =[];                                      % sliding window of data
177    end
178    methods
179        function obj=validate(obj)                            % prepare all internal objects for use
180            obj.apost_pdf = mexDirac;
181            obj.apost_pdf.point = [0;0];
182            obj.log_evidence = 0;                             % evidence is not computed!
183        end
184
185        function dims=dimensions(obj)
186            %please fill: dims = [size_of_posterior size_of_data size_of_condition]
187            dims = [2,1,0]                                    % we have: [2d parameters, 1d observations, 0d condition]
188        end
189        function obj=bayes(obj,dt,cond)                       % approximate bayes rule
190            if size(obj.data_window,2)>=obj.max_window_length
191                obj.data_window = [dt obj.data_window(1:end-1)];
192            else
193                obj.data_window = [dt obj.data_window];
194            end
195            % transform old estimate into new estimate
196            m_hat = mean(obj.data_window);
197            b_hat = sum(abs(obj.data_window-m_hat))/ size(obj.data_window,2);
198            obj.apost_pdf.point = [m_hat; b_hat];             % store result in psoterior pdf
199        end
200    end
201
202end
203\endcode
204
205\section dev_mat_alt Alternative representations of this estimator
206
207Note, that the estimator above is rather coarse approximation of the Bayes rule. Indeed, in this case, the sliding window represnts a statistics on the posterior parameters,
208and thus it should be part of the posterior, which can thus evaluate more moments that the current Dirac density.
209
210Such a design is certainly possible, however, it would result in much more complex calculations and more complex coding. If that is the aim, a new class can be designed.
211
212\section dev_mat_pred Predictors and simulation
213
214Note, that the Laplacian likelihood model of the example is not needed for its estimation. The very basic estimator introduced above does not need to know about existence of the likelihood function.
215
216However, existence of this distribution is required by related tasks of data simulation and data prediction.
217Therefore, a class mexLaplace was designed, which implements all operations of pdf in the same manner as the in \ref dev_mat_pdf.
218
219This pdf is then used in the mexLaplaceBM in function that creates predictors:
220\code
221function p=epredictor(obj,cond)                       % when predictive density is needed approximate it by Laplace with point estimates
222   % return predictive density (max likelihood)
223   p = mexLaplace;
224   p.mu = obj.apost_pdf.point(1);
225   p.b  = obj.apost_pdf.point(2);
226        % do not forget to validate
227   p=p.validate;
228        % assign descriptions
229        p.rv = yrv;
230        % rvc is empty be default
231end
232\endcode
233 
234Also, the class mexLaplace can be used to generate data, see bdmtoolbox/mex/epdfds_example.m
235
236Since it provides all required methods, it can be combined with other pdfs in chain rules, as in pdfds_example.m
237
238Similar extensions to mexEpdf and mexBM can be made to other classes. See library/mex/mex_Epdf.h and library/mex/mex_BM.h how it is done.
239As always, feel free to contact the maintainer for help.
240*/
Note: See TracBrowser for help on using the browser.