root/library/doc/tutorial/040devguide_matlab.dox @ 1064

Revision 988, 11.1 kB (checked in by smidl, 14 years ago)

doc

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