Context Navigation

particles.h @ 883

Revision 870, 11.7 kB (checked in by mido, 14 years ago)
LOG_LEVEL final cut (or rather semifinal? I hope to improve work with ids soon) and also it rests to adapt applications, work is in progress
Property svn:eol-style set to `native`

Line
1	/*!
2	\file
3	\brief Bayesian Filtering using stochastic sampling (Particle Filters)
4	\author Vaclav Smidl.
5
6	-----------------------------------
7	BDM++ - C++ library for Bayesian Decision Making under Uncertainty
8
9	Using IT++ for numerical operations
10	-----------------------------------
11	*/
12
13	#ifndef PARTICLES_H
14	#define PARTICLES_H
15
16
17	#include "../estim/arx_ext.h"
18
19	namespace bdm {
20
21	/*!
22	* \brief Trivial particle filter with proposal density equal to parameter evolution model.
23
24	Posterior density is represented by a weighted empirical density (\c eEmp ).
25	*/
26
27	class PF : public BM {
28	//! \var log_level_enums weights
29	//! all weightes will be logged
30
31	//! \var log_level_enums samples
32	//! all samples will be logged
33	LOG_LEVEL(PF,weights,samples);
34
35	protected:
36	//!number of particles;
37	int n;
38	//!posterior density
39	eEmp est;
40	//! pointer into \c eEmp
41	vec &_w;
42	//! pointer into \c eEmp
43	Array<vec> &_samples;
44	//! Parameter evolution model
45	shared_ptr<pdf> par;
46	//! Observation model
47	shared_ptr<pdf> obs;
48	//! internal structure storing loglikelihood of predictions
49	vec lls;
50
51	//! which resampling method will be used
52	RESAMPLING_METHOD resmethod;
53	//! resampling threshold; in this case its meaning is minimum ratio of active particles
54	//! For example, for 0.5 resampling is performed when the numebr of active aprticles drops belo 50%.
55	double res_threshold;
56
57	//! \name Options
58	//!@{
59	//!@}
60
61	public:
62	//! \name Constructors
63	//!@{
64	PF ( ) : est(), _w ( est._w() ), _samples ( est._samples() ) {
65	};
66
67	void set_parameters ( int n0, double res_th0 = 0.5, RESAMPLING_METHOD rm = SYSTEMATIC ) {
68	n = n0;
69	res_threshold = res_th0;
70	resmethod = rm;
71	};
72	void set_model ( shared_ptr<pdf> par0, shared_ptr<pdf> obs0 ) {
73	par = par0;
74	obs = obs0;
75	// set values for posterior
76	est.set_rv ( par->_rv() );
77	};
78	void set_statistics ( const vec w0, const epdf &epdf0 ) {
79	est.set_statistics ( w0, epdf0 );
80	};
81	void set_statistics ( const eEmp &epdf0 ) {
82	bdm_assert_debug ( epdf0._rv().equal ( par->_rv() ), "Incompatible input" );
83	est = epdf0;
84	};
85	//!@}
86
87	//! Set posterior density by sampling from epdf0
88	//! bayes I - generate samples and add their weights to lls
89	virtual void bayes_gensmp ( const vec &ut );
90	//! bayes II - compute weights of the
91	virtual void bayes_weights();
92	//! important part of particle filtering - decide if it is time to perform resampling
93	virtual bool do_resampling() {
94	double eff = 1.0 / ( _w * _w );
95	return eff < ( res_threshold*n );
96	}
97	void bayes ( const vec &yt, const vec &cond );
98	//!access function
99	vec& __w() {
100	return _w;
101	}
102	//!access function
103	vec& _lls() {
104	return lls;
105	}
106	//!access function
107	RESAMPLING_METHOD _resmethod() const {
108	return resmethod;
109	}
110	//! return correctly typed posterior (covariant return)
111	const eEmp& posterior() const {
112	return est;
113	}
114
115	/*! configuration structure for basic PF
116	\code
117	parameter_pdf = pdf_class; // parameter evolution pdf
118	observation_pdf = pdf_class; // observation pdf
119	prior = epdf_class; // prior probability density
120	--- optional ---
121	n = 10; // number of particles
122	resmethod = 'systematic', or 'multinomial', or 'stratified'
123	// resampling method
124	res_threshold = 0.5; // resample when active particles drop below 50%
125	\endcode
126	*/
127	void from_setting ( const Setting &set ) {
128	BM::from_setting ( set );
129	par = UI::build<pdf> ( set, "parameter_pdf", UI::compulsory );
130	obs = UI::build<pdf> ( set, "observation_pdf", UI::compulsory );
131
132	prior_from_set ( set );
133	resmethod_from_set ( set );
134	// set resampling method
135	//set drv
136	//find potential input - what remains in rvc when we subtract rv
137	RV u = par->_rvc().remove_time().subt ( par->_rv() );
138	//find potential input - what remains in rvc when we subtract x_t
139	RV obs_u = obs->_rvc().remove_time().subt ( par->_rv() );
140
141	u.add ( obs_u ); // join both u, and check if they do not overlap
142
143	set_yrv ( obs->_rv() );
144	rvc = u;
145	}
146	//! auxiliary function reading parameter 'resmethod' from configuration file
147	void resmethod_from_set ( const Setting &set ) {
148	string resmeth;
149	if ( UI::get ( resmeth, set, "resmethod", UI::optional ) ) {
150	if ( resmeth == "systematic" ) {
151	resmethod = SYSTEMATIC;
152	} else {
153	if ( resmeth == "multinomial" ) {
154	resmethod = MULTINOMIAL;
155	} else {
156	if ( resmeth == "stratified" ) {
157	resmethod = STRATIFIED;
158	} else {
159	bdm_error ( "Unknown resampling method" );
160	}
161	}
162	}
163	} else {
164	resmethod = SYSTEMATIC;
165	};
166	if ( !UI::get ( res_threshold, set, "res_threshold", UI::optional ) ) {
167	res_threshold = 0.5;
168	}
169	//validate();
170	}
171	//! load prior information from set and set internal structures accordingly
172	void prior_from_set ( const Setting & set ) {
173	shared_ptr<epdf> pri = UI::build<epdf> ( set, "prior", UI::compulsory );
174
175	eEmp test_emp = dynamic_cast<eEmp> ( & ( *pri ) );
176	if ( test_emp ) { // given pdf is sampled
177	est = *test_emp;
178	} else {
179	//int n;
180	if ( !UI::get ( n, set, "n", UI::optional ) ) {
181	n = 10;
182	}
183	// sample from prior
184	set_statistics ( ones ( n ) / n, *pri );
185	}
186	n = est._w().length();
187	lls=zeros(n);
188	}
189
190	void validate() {
191	BM::validate();
192	bdm_assert ( par, "PF::parameter_pdf not specified." );
193	n = _w.length();
194	lls = zeros ( n );
195
196	if ( par->_rv()._dsize() > 0 ) {
197	bdm_assert ( par->_rv()._dsize() == est.dimension(), "Mismatch of RV and dimension of posterior" );
198	}
199	}
200	//! resample posterior density (from outside - see MPF)
201	void resample ( ivec &ind ) {
202	est.resample ( ind, resmethod );
203	}
204	//! access function
205	Array<vec>& __samples() {
206	return _samples;
207	}
208
209	virtual double logpred ( const vec &yt ) const NOT_IMPLEMENTED(0);
210
211	virtual epdf* epredictor() const NOT_IMPLEMENTED(NULL);
212
213	virtual pdf* predictor() const NOT_IMPLEMENTED(NULL);
214	};
215	UIREGISTER ( PF );
216
217	/*!
218	\brief Marginalized Particle filter
219
220	A composition of particle filter with exact (or almost exact) bayesian models (BMs).
221	The Bayesian models provide marginalized predictive density. Internaly this is achieved by virtual class MPFpdf.
222	*/
223
224	class MPF : public BM {
225	//! \var log_level_enums means
226	//! TODO DOPLNIT
227	LOG_LEVEL(MPF,means);
228
229	protected:
230	//! particle filter on non-linear variable
231	shared_ptr<PF> pf;
232	//! Array of Bayesian models
233	Array<BM*> BMs;
234
235	//! internal class for pdf providing composition of eEmp with external components
236
237	class mpfepdf : public epdf {
238	//! pointer to particle filter
239	shared_ptr<PF> &pf;
240	//! pointer to Array of BMs
241	Array<BM*> &BMs;
242	public:
243	//! constructor
244	mpfepdf ( shared_ptr<PF> &pf0, Array<BM*> &BMs0 ) : epdf(), pf ( pf0 ), BMs ( BMs0 ) { };
245	//! a variant of set parameters - this time, parameters are read from BMs and pf
246	void read_parameters() {
247	rv = concat ( pf->posterior()._rv(), BMs ( 0 )->posterior()._rv() );
248	dim = pf->posterior().dimension() + BMs ( 0 )->posterior().dimension();
249	bdm_assert_debug ( dim == rv._dsize(), "Wrong name " );
250	}
251	vec mean() const;
252
253	vec variance() const;
254
255	void qbounds ( vec &lb, vec &ub, double perc = 0.95 ) const;
256
257	vec sample() const NOT_IMPLEMENTED(0);
258
259	double evallog ( const vec &val ) const NOT_IMPLEMENTED(0);
260	};
261
262	//! Density joining PF.est with conditional parts
263	mpfepdf jest;
264
265	//! datalink from global yt and cond (Up) to BMs yt and cond (Down)
266	datalink_m2m this2bm;
267	//! datalink from global yt and cond (Up) to PFs yt and cond (Down)
268	datalink_m2m this2pf;
269	//!datalink from PF part to BM
270	datalink_part pf2bm;
271
272	public:
273	//! Default constructor.
274	MPF () : jest ( pf, BMs ) {};
275	//! set all parameters at once
276	void set_pf ( shared_ptr<pdf> par0, int n0, RESAMPLING_METHOD rm = SYSTEMATIC ) {
277	if (!pf) pf=new PF;
278	pf->set_model ( par0, par0 ); // <=== nasty!!!
279	pf->set_parameters ( n0, rm );
280	pf->set_rv(par0->_rv());
281	BMs.set_length ( n0 );
282	}
283	//! set a prototype of BM, copy it to as many times as there is particles in pf
284	void set_BM ( const BM &BMcond0 ) {
285
286	int n = pf->__w().length();
287	BMs.set_length ( n );
288	// copy
289	//BMcond0 .condition ( pf->posterior()._sample ( 0 ) );
290	for ( int i = 0; i < n; i++ ) {
291	BMs ( i ) = (BM*) BMcond0._copy();
292	}
293	};
294
295	void bayes ( const vec &yt, const vec &cond );
296
297	const epdf& posterior() const {
298	return jest;
299	}
300
301	//!Access function
302	const BM* _BM ( int i ) {
303	return BMs ( i );
304	}
305	PF& _pf() {return *pf;}
306
307
308	virtual double logpred ( const vec &yt ) const NOT_IMPLEMENTED(0);
309
310	virtual epdf* epredictor() const NOT_IMPLEMENTED(NULL);
311
312	virtual pdf* predictor() const NOT_IMPLEMENTED(NULL);
313
314
315	/*! configuration structure for basic PF
316	\code
317	BM = BM_class; // Bayesian filtr for analytical part of the model
318	parameter_pdf = pdf_class; // transitional pdf for non-parametric part of the model
319	prior = epdf_class; // prior probability density
320	--- optional ---
321	n = 10; // number of particles
322	resmethod = 'systematic', or 'multinomial', or 'stratified'
323	// resampling method
324	\endcode
325	*/
326	void from_setting ( const Setting &set ) {
327	BM::from_setting( set );
328
329	shared_ptr<pdf> par = UI::build<pdf> ( set, "parameter_pdf", UI::compulsory );
330
331	pf = new PF;
332	// prior must be set before BM
333	pf->prior_from_set ( set );
334	pf->resmethod_from_set ( set );
335	pf->set_model ( par, par ); // too hackish!
336
337	shared_ptr<BM> BM0 = UI::build<BM> ( set, "BM", UI::compulsory );
338	set_BM ( *BM0 );
339
340	//set drv
341	//??set_yrv(concat(BM0->_yrv(),u) );
342	set_yrv ( BM0->_yrv() );
343	rvc = BM0->_rvc().subt ( par->_rv() );
344	//find potential input - what remains in rvc when we subtract rv
345	RV u = par->_rvc().subt ( par->_rv().copy_t ( -1 ) );
346	rvc.add ( u );
347	dimc = rvc._dsize();
348	validate();
349	}
350
351	void validate() {
352	BM::validate();
353	try {
354	pf->validate();
355	} catch ( std::exception ) {
356	throw UIException ( "Error in PF part of MPF:" );
357	}
358	jest.read_parameters();
359	this2bm.set_connection ( BMs ( 0 )->_yrv(), BMs ( 0 )->_rvc(), yrv, rvc );
360	this2pf.set_connection ( pf->_yrv(), pf->_rvc(), yrv, rvc );
361	pf2bm.set_connection ( BMs ( 0 )->_rvc(), pf->posterior()._rv() );
362	}
363	};
364	UIREGISTER ( MPF );
365
366	/*! ARXg for estimation of state-space variances
367	*/
368	class MPF_ARXg :public BM{
369	protected:
370	shared_ptr<PF> pf;
371	//! pointer to Array of BMs
372	Array<ARX*> BMso;
373	//! pointer to Array of BMs
374	Array<ARX*> BMsp;
375	//!parameter evolution
376	shared_ptr<fnc> g;
377	//!observation function
378	shared_ptr<fnc> h;
379
380	public:
381	void bayes(const vec &yt, const vec &cond );
382	void from_setting(const Setting &set) ;
383	void validate() {
384	bdm_assert(g->dimensionc()==g->dimension(),"not supported yet");
385	bdm_assert(h->dimensionc()==g->dimension(),"not supported yet");
386	}
387
388	double logpred(const vec &cond) const NOT_IMPLEMENTED(0.0);
389	epdf* epredictor() const NOT_IMPLEMENTED(NULL);
390	pdf* predictor() const NOT_IMPLEMENTED(NULL);
391	const epdf& posterior() const {return pf->posterior();};
392
393	void log_register( logger &L, const string &prefix ){
394	BM::log_register(L,prefix);
395	logrec->ids.set_size ( 3 );
396	logrec->ids(1)= L.add_vector(RV("Q",dimension()*dimension()), prefix+L.prefix_sep()+"Q");
397	logrec->ids(2)= L.add_vector(RV("R",dimensiony()*dimensiony()), prefix+L.prefix_sep()+"R");
398
399	};
400	void log_write() const {
401	BM::log_write();
402	mat mQ=zeros(dimension(),dimension());
403	mat pom=zeros(dimension(),dimension());
404	mat mR=zeros(dimensiony(),dimensiony());
405	mat pom2=zeros(dimensiony(),dimensiony());
406	mat dum;
407	const vec w=pf->posterior()._w();
408	for (int i=0; i<w.length(); i++){
409	BMsp(i)->posterior().mean_mat(dum,pom);
410	mQ += w(i) * pom;
411	BMso(i)->posterior().mean_mat(dum,pom2);
412	mR += w(i) * pom2;
413
414	}
415	logrec->L.log_vector ( logrec->ids ( 1 ), cvectorize(mQ) );
416	logrec->L.log_vector ( logrec->ids ( 2 ), cvectorize(mR) );
417
418	}
419	};
420	UIREGISTER(MPF_ARXg);
421
422
423	}
424	#endif // KF_H
425

Note: See TracBrowser for help on using the browser.

Context Navigation

root/library/bdm/estim/particles.h @ 883

Download in other formats: