Context Navigation

← Previous Change
Next Change →

Changeset 737 for library/bdm/estim

Timestamp:

11/25/09 12:14:38 (15 years ago)

Author:

mido

Message:

ASTYLER RUN OVER THE WHOLE LIBRARY, JUPEE

Location:

library/bdm/estim

Files:

: 10 modified

arx.cpp (modified) (6 diffs)
arx.h (modified) (10 diffs)
arx_straux.cpp (modified) (9 diffs)
arx_straux.h (modified) (2 diffs)
kalman.cpp (modified) (13 diffs)
kalman.h (modified) (12 diffs)
mixtures.cpp (modified) (1 diff)
mixtures.h (modified) (2 diffs)
particles.cpp (modified) (4 diffs)
particles.h (modified) (16 diffs)

Legend:

: Unmodified
: Added
: Removed

library/bdm/estim/arx.cpp

r723	r737
3	3
4	4	void ARX::bayes_weighted ( const vec &yt, const vec &cond, const double w ) {
5
6		bdm_assert_debug~~(yt.length()>=dimy,"ARX::bayes yt is smaller then dimc");~~
7		bdm_assert_debug~~(cond.length()>=dimc,"ARX::bayes cond is smaller then dimc");~~
	5
	6	bdm_assert_debug ( yt.length() >= dimy, "ARX::bayes yt is smaller then dimc" );
	7	bdm_assert_debug ( cond.length() >= dimc, "ARX::bayes cond is smaller then dimc" );
8	8	double lnc;
9	9	//cache
10		ldmat &V~~=est._V();~~
11		double &nu=est._nu();
12
13		dyad.set_subvector~~(0,yt~~);
14		dyad.set_subvector~~(dimy,cond~~);
	10	ldmat &V = est._V();
	11	double &nu = est._nu();
	12
	13	dyad.set_subvector ( 0, yt );
	14	dyad.set_subvector ( dimy, cond );
15	15	// possible "1" is there from the beginning
16
	16
17	17	if ( frg < 1.0 ) {
18	18	est.pow ( frg ); // multiply V and nu
19	19
20
	20
21	21	//stabilize
22		ldmat V0=alter_est._V(); //$ copy
23		double &nu0=alter_est._nu();
24
25		V0*=(1-frg);
	22	ldmat V0 = alter_est._V(); //$ copy
	23	double &nu0 = alter_est._nu();
	24
	25	V0 *= ( 1 - frg );
26	26	V += V0; //stabilization
27		nu +=(1-frg)*nu0;
28
	27	nu += ( 1 - frg ) * nu0;
	28
29	29	// recompute loglikelihood of new "prior"
30	30	if ( evalll ) {
…	…
50	50	double lll;
51	51	vec dyad_p = dyad;
52		dyad_p.set_subvector~~(0,yt~~);
	52	dyad_p.set_subvector ( 0, yt );
53	53
54	54	if ( frg < 1.0 ) {
…	…
101	101	mlstudent* ARX::predictor_student ( ) const {
102	102	const ldmat &V = posterior()._V();
103
	103
104	104	mat mu ( dimy, V.rows() - dimy );
105	105	mat R ( dimy, dimy );
…	…
114	114
115	115
116		if ( have_constant) { // no constant term
	116	if ( have_constant ) { // no constant term
117	117	//Assume the constant term is the last one:
118	118	if ( mu.cols() > 1 ) {
…	…
206	206	// rgrlen - including constant!!!
207	207	dimc = rrv->_dsize();
208
	208
209	209	yrv = *yrv_;
210	210	rvc = *rrv;
211
	211
212	212	string opt;
213		if ( UI::get~~(opt, set, "options", UI::optional~~) ) {
214		BM::set_options~~(opt~~);
	213	if ( UI::get ( opt, set, "options", UI::optional ) ) {
	214	BM::set_options ( opt );
215	215	}
216	216	int constant;
217		if (~~!UI::get(constant, set, "constant", UI::optional))~~{
218		have_constant=true;
219		} else {
220		have_constant~~=constant>~~0;
221		}
222		int rgrlen = dimc~~+int(have_constant==true~~);
	217	if ( !UI::get ( constant, set, "constant", UI::optional ) ) {
	218	have_constant = true;
	219	} else {
	220	have_constant = constant > 0;
	221	}
	222	int rgrlen = dimc + int ( have_constant == true );
223	223
224	224	//init
225		shared_ptr<egiw> pri~~=UI::build<egiw>(set, "prior", UI::optional~~);
226		if (~~pri~~) {
227		bdm_assert~~(pri->_dimx()==dimy,"prior is not compatible"~~);
228		bdm_assert~~(pri->_V().rows()==dimy+rgrlen,"prior is not compatible"~~);
229		est.set_parameters~~( pri->_dimx(),pri->_V(), pri->_nu()~~);
230		}~~else~~{
231		est.set_parameters~~( dimy, zeros(dimy+rgrlen)~~);
232		set_prior_default~~(est~~);
233		}
234
235		shared_ptr<egiw> alt~~=UI::build<egiw>(set, "alternative", UI::optional~~);
236		if (~~alt~~) {
237		bdm_assert~~(alt->_dimx()==dimy,"alternative is not compatible"~~);
238		bdm_assert~~(alt->_V().rows()==dimy+rgrlen,"alternative is not compatible"~~);
239		alter_est.set_parameters~~( alt->_dimx(),alt->_V(), alt->_nu()~~);
	225	shared_ptr<egiw> pri = UI::build<egiw> ( set, "prior", UI::optional );
	226	if ( pri ) {
	227	bdm_assert ( pri->_dimx() == dimy, "prior is not compatible" );
	228	bdm_assert ( pri->_V().rows() == dimy + rgrlen, "prior is not compatible" );
	229	est.set_parameters ( pri->_dimx(), pri->_V(), pri->_nu() );
	230	} else {
	231	est.set_parameters ( dimy, zeros ( dimy + rgrlen ) );
	232	set_prior_default ( est );
	233	}
	234
	235	shared_ptr<egiw> alt = UI::build<egiw> ( set, "alternative", UI::optional );
	236	if ( alt ) {
	237	bdm_assert ( alt->_dimx() == dimy, "alternative is not compatible" );
	238	bdm_assert ( alt->_V().rows() == dimy + rgrlen, "alternative is not compatible" );
	239	alter_est.set_parameters ( alt->_dimx(), alt->_V(), alt->_nu() );
240	240	} else {
241	241	alter_est = est;
…	…
247	247
248	248	set_parameters ( frg );
249
	249
250	250	//name results (for logging)
251		shared_ptr<RV> rv_par~~=UI::build<RV>(set, "rv_param",~~UI::optional );
252		if (~~!rv_par)~~{
	251	shared_ptr<RV> rv_par = UI::build<RV> ( set, "rv_param", UI::optional );
	252	if ( !rv_par ) {
253	253	est.set_rv ( RV ( "{theta r }", vec_2 ( dimyrgrlen, dimydimy ) ) );
254	254	} else {

library/bdm/estim/arx.h

r723	r737
54	54	//! \name Constructors
55	55	//!@{
56		ARX ( const double frg0 = 1.0 ) : BMEF ( frg0 ), have_constant~~(true~~), dyad(), est(), alter_est() {};
57		ARX ( const ARX &A0 ) : BMEF (~~A0), have_constant(A0.have_constant), dyad(A0.dyad),est(A0.est),alter_est(A0.alter_est~~) { };
	56	ARX ( const double frg0 = 1.0 ) : BMEF ( frg0 ), have_constant ( true ), dyad(), est(), alter_est() {};
	57	ARX ( const ARX &A0 ) : BMEF ( A0 ), have_constant ( A0.have_constant ), dyad ( A0.dyad ), est ( A0.est ), alter_est ( A0.alter_est ) { };
58	58	ARX* _copy_() const;
59	59	void set_parameters ( double frg0 ) {
…	…
61	61	}
62	62	void set_constant ( bool const0 ) {
63		have_constant=const0;
	63	have_constant = const0;
64	64	}
65	65	void set_statistics ( int dimy0, const ldmat V0, double nu0 = -1.0 ) {
…	…
77	77
78	78	//! Weighted Bayes \f$ dt = [y_t psi_t] \f$.
79		void bayes_weighted ( const vec &yt, const vec &cond~~=empty_vec, const double w=~~1.0 );
80		void bayes~~( const vec &yt, const vec &cond=~~empty_vec ) {
81		bayes_weighted ( yt,cond, 1.0 );
	79	void bayes_weighted ( const vec &yt, const vec &cond = empty_vec, const double w = 1.0 );
	80	void bayes ( const vec &yt, const vec &cond = empty_vec ) {
	81	bayes_weighted ( yt, cond, 1.0 );
82	82	};
83	83	double logpred ( const vec &yt ) const;
…	…
100	100	template<class sq_T>
101	101	shared_ptr<mlnorm<sq_T> > ml_predictor() const;
102		//! fast version of predicto
	102	//! fast version of predicto
103	103	template<class sq_T>
104		void ml_predictor_update~~(mlnorm<sq_T> &pred~~) const;
	104	void ml_predictor_update ( mlnorm<sq_T> &pred ) const;
105	105	mlstudent* predictor_student() const;
106	106	//! Brute force structure estimation.\return indeces of accepted regressors.
…	…
112	112	//!\name Access attributes
113	113	//!@{
114		//! return correctly typed posterior (covariant return)
115		const egiw& posterior() const {
	114	//! return correctly typed posterior (covariant return)
	115	const egiw& posterior() const {
116	116	return est;
117	117	}
…	…
129	129	prior = {class='egiw',...}; // Prior density, when given default is used instead
130	130	alternative = {class='egiw',...}; // Alternative density in stabilized estimation, when not given prior is used
131
	131
132	132	frg = 1.0; // forgetting, default frg=1.0
133	133
134		rv_param = RV({names_of_parameters}} // description of parametetr names
	134	rv_param = RV({names_of_parameters}} // description of parametetr names
135	135	// default: ["theta_i" and "r_i"]
136	136	\endcode
…	…
140	140	void validate() {
141	141	//if dimc not set set it from V
142		if (~~dimc==0)~~{
143		dimc = posterior()._V().rows()~~-dimy-int(have_constant==true~~);
	142	if ( dimc == 0 ) {
	143	dimc = posterior()._V().rows() - dimy - int ( have_constant == true );
144	144	}
145
146		if (~~have_constant~~) {
147		dyad = ones~~(dimy+dimc+1~~);
148		} else {
149		dyad = zeros~~(dimy+dimc~~);
	145
	146	if ( have_constant ) {
	147	dyad = ones ( dimy + dimc + 1 );
	148	} else {
	149	dyad = zeros ( dimy + dimc );
150	150	}
151
152		}
153		//! function sets prior and alternative density
154		void set_prior~~(const RV &drv, egiw &prior)~~{
	151
	152	}
	153	//! function sets prior and alternative density
	154	void set_prior ( const RV &drv, egiw &prior ) {
155	155	//TODO check ranges in RV and build prior
156	156	};
157	157	//! build default prior and alternative when all values are set
158		void set_prior_default~~(egiw &prior)~~{
159		//assume
160		vec dV0~~(prior._V().rows()~~);
161		dV0.set_subvector~~(0,prior._dimx()-1, 1.0~~);
162		dV0.set_subvector~~(prior._dimx(),dV0.length()-1, 1e-5~~);
163
164		prior.set_parameters~~(prior._dimx(),ldmat(dV0)~~);
	158	void set_prior_default ( egiw &prior ) {
	159	//assume
	160	vec dV0 ( prior._V().rows() );
	161	dV0.set_subvector ( 0, prior._dimx() - 1, 1.0 );
	162	dV0.set_subvector ( prior._dimx(), dV0.length() - 1, 1e-5 );
	163
	164	prior.set_parameters ( prior._dimx(), ldmat ( dV0 ) );
165	165	}
166	166	};
…	…
174	174	The symbol \f$ \phi \f$ is assumed to be the last of the conditioning variables.
175	175	*/
176		class ARXfrg : public ARX{
177		public:
178		ARXfrg():ARX(){};
179		//! copy constructor
180		ARXfrg(const ARXfrg &A0):ARX(A0){};
181		ARXfrg* _copy_() const {ARXfrg A = new ARXfrg(this); return A;}
182
183		void bayes(const vec &val, const vec &cond){
184		frg = cond(dimc-1); // last in cond is phi
185		ARX::bayes(val,cond);
	176	class ARXfrg : public ARX {
	177	public:
	178	ARXfrg() : ARX() {};
	179	//! copy constructor
	180	ARXfrg ( const ARXfrg &A0 ) : ARX ( A0 ) {};
	181	ARXfrg* _copy_() const {
	182	ARXfrg A = new ARXfrg ( this );
	183	return A;
	184	}
	185
	186	void bayes ( const vec &val, const vec &cond ) {
	187	frg = cond ( dimc - 1 ); // last in cond is phi
	188	ARX::bayes ( val, cond );
186	189	}
187	190	void validate() {
188	191	ARX::validate();
189		rvc.add~~(RV("{phi }",vec_1(1))~~);
190		dimc +=1;
	192	rvc.add ( RV ( "{phi }", vec_1 ( 1 ) ) );
	193	dimc += 1;
191	194	}
192	195	};
193		UIREGISTER~~(ARXfrg~~);
	196	UIREGISTER ( ARXfrg );
194	197
195	198
…	…
198	201	shared_ptr< mlnorm<sq_T> > ARX::ml_predictor ( ) const {
199	202	shared_ptr< mlnorm<sq_T> > tmp = new mlnorm<sq_T> ( );
200		tmp->set_rv~~(yrv~~);
201		tmp->set_rvc~~(_rvc()~~);
202
203		ml_predictor_update~~(*tmp~~);
	203	tmp->set_rv ( yrv );
	204	tmp->set_rvc ( _rvc() );
	205
	206	ml_predictor_update ( *tmp );
204	207	tmp->validate();
205	208	return tmp;
…	…
207	210
208	211	template<class sq_T>
209		void ARX::ml_predictor_update~~(mlnorm<sq_T> &pred~~) const {
	212	void ARX::ml_predictor_update ( mlnorm<sq_T> &pred ) const {
210	213	mat mu ( dimy, posterior()._V().rows() - dimy );
211	214	mat R ( dimy, dimy );
212
	215
213	216	est.mean_mat ( mu, R ); //mu =
214	217	mu = mu.T();
215	218	//correction for student-t -- TODO check if correct!!
216
217		if ( have_constant) { // constant term
	219
	220	if ( have_constant ) { // constant term
218	221	//Assume the constant term is the last one:
219	222	pred.set_parameters ( mu.get_cols ( 0, mu.cols() - 2 ), mu.get_col ( mu.cols() - 1 ), sq_T ( R ) );
220	223	} else {
221		pred.set_parameters ( mu, zeros ( dimy ), sq_T( R ) );
	224	pred.set_parameters ( mu, zeros ( dimy ), sq_T ( R ) );
222	225	}
223	226	}

library/bdm/estim/arx_straux.cpp

r706	r737
7	7
8	8
9
	9
10	10	void str_bitset ( bvec &out, ivec ns, int nbits ) {
11
12		for ( int i = 0; i < ns.length(); i++ ) {
	11
	12	for ( int i = 0; i < ns.length(); i++ ) {
13	13	out ( ns ( i ) - 2 ) = 1;
14	14	}
…	…
158	158	L.swap_cols ( i, j );
159	159
160
	160
161	161	//% We must be working with LINES of matrix L !
162	162
163		mat r = L.get_row ( i );
	163	mat r = L.get_row ( i );
164	164	r.transpose();
165	165	mat f = L.get_row ( j );
166	166	f.transpose();
167	167
168		double Dr = d ( i );
	168	double Dr = d ( i );
169	169	double Df = d ( j );
170	170
171	171	sedydr ( r, f, Dr, Df, j );
172	172
173		double r0 = r ( i, 0 );
	173	double r0 = r ( i, 0 );
174	174	Dr = Dr * r0 * r0;
175	175	r = r / r0;
176	176
177		mat pom_mat = r.transpose();
	177	mat pom_mat = r.transpose();
178	178	L.set_row ( i, pom_mat.get_row ( 0 ) );
179	179	pom_mat = f.transpose();
…	…
190	190	void str_bitres ( bvec &out, ivec ns, int nbits ) {
191	191
192		for ( int i = 0; i < ns.length(); i++ ) {
	192	for ( int i = 0; i < ns.length(); i++ ) {
193	193	out ( ns ( i ) - 2 ) = 0;
194	194	}
…	…
197	197	str_aux sestrremove ( str_aux in, ivec removed_elements ) {
198	198	//% Removes elements from regressor
199
	199
200	200	int n_strL = length ( in.strL );
201	201	str_aux out = in;
202
203		for ( int i = 0; i < removed_elements.length(); i++ ) {
	202
	203	for ( int i = 0; i < removed_elements.length(); i++ ) {
204	204
205	205	int f = removed_elements ( i );
…	…
216	216	//% END
217	217	}
218
219
	218
	219
220	220	}
221	221	out.posit1 = ( find ( out.strL == 1 ) ) ( 0 ) + 1;
…	…
252	252	//logliks = [global_best.loglik];
253	253
254		for ( int j = 1; j < global_best.length(); j++ ) {
	254	for ( int j = 1; j < global_best.length(); j++ ) {
255	255	if ( global_best ( j ).loglik < global_best ( i ).loglik ) {
256	256	i = j;
257	257	}
258	258	}
259		if ( global_best ( i ).loglik < newone.loglik ) {
	259	if ( global_best ( i ).loglik < newone.loglik ) {
260	260	// if ~any(logliks == new.loglik);
261	261	addit = 1;
262	262
263	263
264		~~for (int j = 0; j < global_best.length(); j++)~~{
265		~~if (newone.bitstr == global_best(j).bitstr)~~{
266		addit = 0;
267		break;
268		}
269		}
270		if ( addit ) {
271		global_best ( i ) = newone;
	264	for ( int j = 0; j < global_best.length(); j++ ) {
	265	if ( newone.bitstr == global_best ( j ).bitstr ) {
	266	addit = 0;
	267	break;
	268	}
	269	}
	270	if ( addit ) {
	271	global_best ( i ) = newone;
272	272	// DEBUGging print:
273	273	// fprintf('ADDED structure, add_new: %s, loglik=%g\n', strPrintstr(new), new.loglik);
274		}
275		}
	274	}
	275	}
276	276	} else
277	277	global_best = concat ( global_best, newone );
…	…
413	413	else {
414	414	//% start from random structure
415
416
417
418		~~ivec last_str = find ( to_bvec( concat( 0, floor_i ( 2 * randun ( n_data - 1 )~~) ) ) ) + 1;// this creates random vector consisting of indexes, and sorted
419		~~last = sestrremove ( full, setdiff ( all_str, concat( concat~~( 1 , last_str ), empty.strRgr ) ) );
	415
	416
	417
	418	ivec last_str = find ( to_bvec ( concat ( 0, floor_i ( 2 * randun ( n_data - 1 ) ) ) ) ) + 1;// this creates random vector consisting of indexes, and sorted
	419	last = sestrremove ( full, setdiff ( all_str, concat ( concat ( 1 , last_str ), empty.strRgr ) ) );
420	420
421	421	}
…	…
448	448	}
449	449	//% Nesting by adding elements (enrichment)
450		ivec added_items = setdiff( last.strMis, belief_out );
	450	ivec added_items = setdiff ( last.strMis, belief_out );
451	451	ivec added_item;
452	452
…	…
469	469	//% Making best structure last structure.
470	470	last = best;
471		}
	471	}
472	472
473	473

library/bdm/estim/arx_straux.h

r646	r737
16	16	namespace bdm {
17	17
18
	18
19	19	struct str_aux {
20	20	vec d0;
…	…
29	29	int posit1; // regressand position
30	30	int nbits; // number of bits available in double
31		bvec bitstr;
	31	bvec bitstr;
32	32	double loglik; // loglikelihood
33		};
34
	33	};
	34
35	35
36	36
37	37	struct str_statistics {
38		long long int allstrs;
39		int nrand;
40		int unique;
41		int to;
42		double cputime_seconds;
43		double itemspeed;
44		int muto;
45		ivec mutos;
46		vec maxmutos;
	38	long long int allstrs;
	39	int nrand;
	40	int unique;
	41	int to;
	42	double cputime_seconds;
	43	double itemspeed;
	44	int muto;
	45	ivec mutos;
	46	vec maxmutos;
47	47	};
48	48
49	49
50		//! Rplication of Ludvik Tesar original straux1 from mixtools straux1
51		ivec straux1~~(ldmat Ld, double nu, ldmat Ld0, double nu0, ivec belief, int nbest, int max_nrep, double lambda, int order_k, Array<str_aux> &rgrsout~~);
	50	//! Rplication of Ludvik Tesar original straux1 from mixtools straux1
	51	ivec straux1 ( ldmat Ld, double nu, ldmat Ld0, double nu0, ivec belief, int nbest, int max_nrep, double lambda, int order_k, Array<str_aux> &rgrsout );
52	52
53	53	}

library/bdm/estim/kalman.cpp

r686	r737
13	13	const vec &u = cond; // in this case cond=ut
14	14	const vec &y = yt;
15
	15
16	16	vec& mu = est._mu();
17	17	mat &P = est._R();
…	…
20	20	//Time update
21	21	mu = A * mu + B * u;
22		P = A * P * A.transpose() + (~~mat)~~Q;
	22	P = A * P * A.transpose() + ( mat ) Q;
23	23
24	24	//Data update
25		_Ry = C * P * C.transpose() + (~~mat)~~R;
	25	_Ry = C * P * C.transpose() + ( mat ) R;
26	26	_K = P * C.transpose() * inv ( _Ry );
27	27	P -= _K * C * P; // P = P -KCP;
…	…
30	30
31	31	if ( evalll ) {
32		ll~~=fy.evallog(y~~);
	32	ll = fy.evallog ( y );
33	33	}
34	34	};
…	…
43	43	phxu = phxu0;
44	44
45		set_dim~~( pfxu0->_dimx()~~);
	45	set_dim ( pfxu0->_dimx() );
46	46	dimy = phxu0->dimension();
47	47	dimc = pfxu0->_dimu();
48		est.set_parameters~~( zeros(dimension()), eye(dimension()~~) );
	48	est.set_parameters ( zeros ( dimension() ), eye ( dimension() ) );
49	49
50	50	A.set_size ( dimension(), dimension() );
…	…
60	60	}
61	61
62		void EKFfull::bayes ( const vec &yt, const vec &cond) {
	62	void EKFfull::bayes ( const vec &yt, const vec &cond ) {
63	63	bdm_assert_debug ( yt.length() == ( dimy ), "EKFull::bayes wrong size of dt" );
64	64	bdm_assert_debug ( cond.length() == ( dimc ), "EKFull::bayes wrong size of dt" );
65
	65
66	66	const vec &u = cond;
67	67	const vec &y = yt; //lazy to change it
…	…
70	70	mat& _Ry = fy._R();
71	71	vec& yp = fy._mu();
72
	72
73	73	pfxu->dfdx_cond ( mu, zeros ( dimc ), A, true );
74	74	phxu->dfdx_cond ( mu, zeros ( dimc ), C, true );
…	…
76	76	//Time update
77	77	mu = pfxu->eval ( mu, u );// Amu + Bu;
78		P = A * P * A.transpose() + (~~mat)~~Q;
	78	P = A * P * A.transpose() + ( mat ) Q;
79	79
80	80	//Data update
81		_Ry = C * P * C.transpose() + (~~mat)~~R;
	81	_Ry = C * P * C.transpose() + ( mat ) R;
82	82	_K = P * C.transpose() * inv ( _Ry );
83	83	P -= _K * C * P; // P = P -KCP;
…	…
86	86
87	87	if ( BM::evalll ) {
88		ll~~=fy.evallog(y~~);
	88	ll = fy.evallog ( y );
89	89	}
90	90	};
…	…
95	95
96	96	( ( StateSpace<chmat>* ) this )->set_parameters ( A0, B0, C0, D0, Q0, R0 );
97
98		_K~~=zeros(dimension(),dimy~~);
99		}
100
101		void KalmanCh::initialize(){
	97
	98	_K = zeros ( dimension(), dimy );
	99	}
	100
	101	void KalmanCh::initialize() {
102	102	preA = zeros ( dimy + dimension() + dimension(), dimy + dimension() );
103	103	// preA.clear();
…	…
107	107
108	108	void KalmanCh::bayes ( const vec &yt, const vec &cond ) {
109		bdm_assert_debug~~(yt.length()==dimy, "yt mismatch"~~);
110		bdm_assert_debug~~(cond.length()==dimc, "yt mismatch"~~);
111
	109	bdm_assert_debug ( yt.length() == dimy, "yt mismatch" );
	110	bdm_assert_debug ( cond.length() == dimc, "yt mismatch" );
	111
112	112	const vec &u = cond;
113	113	const vec &y = yt;
114	114	vec pom ( dimy );
115	115
116		chmat &_P=est._R();
	116	chmat &_P = est._R();
117	117	vec &_mu = est._mu();
118		mat _K~~(dimension(),dimy~~);
119		chmat &_Ry=fy._R();
	118	mat _K ( dimension(), dimy );
	119	chmat &_Ry = fy._R();
120	120	vec &_yp = fy._mu();
121	121	//TODO get rid of Q in qr()!
…	…
156	156	phxu = phxu0;
157	157
158		set_dim~~( pfxu0->_dimx()~~);
	158	set_dim ( pfxu0->_dimx() );
159	159	dimy = phxu0->dimension();
160	160	dimc = pfxu0->_dimu();
161
	161
162	162	vec &_mu = est._mu();
163	163	// if mu is not set, set it to zeros, just for constant terms of A and C
164		if ( _mu.length() != dimension() ) _mu = zeros ( dimension() );
	164	if ( _mu.length() != dimension() ) _mu = zeros ( dimension() );
165	165	A = zeros ( dimension(), dimension() );
166	166	C = zeros ( dimy, dimension() );
…	…
194	194	chmat &_Ry = fy._R();
195	195	vec &_yp = fy._mu();
196
	196
197	197	pfxu->dfdx_cond ( _mu, u, A, false ); //update A by a derivative of fx
198	198	phxu->dfdx_cond ( _mu, u, C, false ); //update A by a derivative of fx
…	…
244	244
245	245	void EKFCh::from_setting ( const Setting &set ) {
246		BM::from_setting~~(set~~);
	246	BM::from_setting ( set );
247	247	shared_ptr<diffbifn> IM = UI::build<diffbifn> ( set, "IM", UI::compulsory );
248	248	shared_ptr<diffbifn> OM = UI::build<diffbifn> ( set, "OM", UI::compulsory );

library/bdm/estim/kalman.h

r723	r737
21	21	//#include <../applications/pmsm/simulator_zdenek/ekf_example/pmsm_mod.h>
22	22
23		namespace bdm
24		{
	23	namespace bdm {
25	24
26	25	/*!
…	…
32	31	*/
33	32	template<class sq_T>
34		class StateSpace
35		{
36		protected:
37		//! Matrix A
38		mat A;
39		//! Matrix B
40		mat B;
41		//! Matrix C
42		mat C;
43		//! Matrix D
44		mat D;
45		//! Matrix Q in square-root form
46		sq_T Q;
47		//! Matrix R in square-root form
48		sq_T R;
49		public:
50		StateSpace() : A(), B(), C(), D(), Q(), R() {}
51		//!copy constructor
52		StateSpace(const StateSpace<sq_T> &S0) : A(S0.A), B(S0.B), C(S0.C), D(S0.D), Q(S0.Q), R(S0.R) {}
53		//! set all matrix parameters
54		void set_parameters (const mat &A0, const mat &B0, const mat &C0, const mat &D0, const sq_T &Q0, const sq_T &R0);
55		//! validation
56		void validate();
57		//! not virtual in this case
58		void from_setting (const Setting &set) {
59		UI::get (A, set, "A", UI::compulsory);
60		UI::get (B, set, "B", UI::compulsory);
61		UI::get (C, set, "C", UI::compulsory);
62		UI::get (D, set, "D", UI::compulsory);
63		mat Qtm, Rtm; // full matrices
64		if(!UI::get(Qtm, set, "Q", UI::optional)){
65		vec dq;
66		UI::get(dq, set, "dQ", UI::compulsory);
67		Qtm=diag(dq);
68		}
69		if(!UI::get(Rtm, set, "R", UI::optional)){
70		vec dr;
71		UI::get(dr, set, "dQ", UI::compulsory);
72		Rtm=diag(dr);
73		}
74		R=Rtm; // automatic conversion to square-root form
75		Q=Qtm;
76
77		validate();
78		}
79		//! access function
80		const mat& _A() const {return A;}
81		//! access function
82		const mat& _B()const {return B;}
83		//! access function
84		const mat& _C()const {return C;}
85		//! access function
86		const mat& _D()const {return D;}
87		//! access function
88		const sq_T& _Q()const {return Q;}
89		//! access function
90		const sq_T& _R()const {return R;}
91		};
92
93		//! Common abstract base for Kalman filters
	33	class StateSpace {
	34	protected:
	35	//! Matrix A
	36	mat A;
	37	//! Matrix B
	38	mat B;
	39	//! Matrix C
	40	mat C;
	41	//! Matrix D
	42	mat D;
	43	//! Matrix Q in square-root form
	44	sq_T Q;
	45	//! Matrix R in square-root form
	46	sq_T R;
	47	public:
	48	StateSpace() : A(), B(), C(), D(), Q(), R() {}
	49	//!copy constructor
	50	StateSpace ( const StateSpace<sq_T> &S0 ) : A ( S0.A ), B ( S0.B ), C ( S0.C ), D ( S0.D ), Q ( S0.Q ), R ( S0.R ) {}
	51	//! set all matrix parameters
	52	void set_parameters ( const mat &A0, const mat &B0, const mat &C0, const mat &D0, const sq_T &Q0, const sq_T &R0 );
	53	//! validation
	54	void validate();
	55	//! not virtual in this case
	56	void from_setting ( const Setting &set ) {
	57	UI::get ( A, set, "A", UI::compulsory );
	58	UI::get ( B, set, "B", UI::compulsory );
	59	UI::get ( C, set, "C", UI::compulsory );
	60	UI::get ( D, set, "D", UI::compulsory );
	61	mat Qtm, Rtm; // full matrices
	62	if ( !UI::get ( Qtm, set, "Q", UI::optional ) ) {
	63	vec dq;
	64	UI::get ( dq, set, "dQ", UI::compulsory );
	65	Qtm = diag ( dq );
	66	}
	67	if ( !UI::get ( Rtm, set, "R", UI::optional ) ) {
	68	vec dr;
	69	UI::get ( dr, set, "dQ", UI::compulsory );
	70	Rtm = diag ( dr );
	71	}
	72	R = Rtm; // automatic conversion to square-root form
	73	Q = Qtm;
	74
	75	validate();
	76	}
	77	//! access function
	78	const mat& _A() const {
	79	return A;
	80	}
	81	//! access function
	82	const mat& _B() const {
	83	return B;
	84	}
	85	//! access function
	86	const mat& _C() const {
	87	return C;
	88	}
	89	//! access function
	90	const mat& _D() const {
	91	return D;
	92	}
	93	//! access function
	94	const sq_T& _Q() const {
	95	return Q;
	96	}
	97	//! access function
	98	const sq_T& _R() const {
	99	return R;
	100	}
	101	};
	102
	103	//! Common abstract base for Kalman filters
94	104	template<class sq_T>
95		class Kalman: public BM, public StateSpace<sq_T>
96		{
97		protected:
98		//! id of output
99		RV yrv;
100		//! Kalman gain
101		mat _K;
102		//!posterior
103		enorm<sq_T> est;
104		//!marginal on data f(y\|y)
105		enorm<sq_T> fy;
106		public:
107		Kalman<sq_T>() : BM(), StateSpace<sq_T>(), yrv(), _K(), est(){}
108		//! Copy constructor
109		Kalman<sq_T>(const Kalman<sq_T> &K0) : BM(K0), StateSpace<sq_T>(K0), yrv(K0.yrv), _K(K0._K), est(K0.est), fy(K0.fy){}
110		//!set statistics of the posterior
111		void set_statistics (const vec &mu0, const mat &P0) {est.set_parameters (mu0, P0); };
112		//!set statistics of the posterior
113		void set_statistics (const vec &mu0, const sq_T &P0) {est.set_parameters (mu0, P0); };
114		//! return correctly typed posterior (covariant return)
115		const enorm<sq_T>& posterior() const {return est;}
116		//! load basic elements of Kalman from structure
117		void from_setting (const Setting &set) {
118		StateSpace<sq_T>::from_setting(set);
119
120		mat P0; vec mu0;
121		UI::get(mu0, set, "mu0", UI::optional);
122		UI::get(P0, set, "P0", UI::optional);
123		set_statistics(mu0,P0);
124		// Initial values
125		UI::get (yrv, set, "yrv", UI::optional);
126		UI::get (rvc, set, "urv", UI::optional);
127		set_yrv(concat(yrv,rvc));
128
129		validate();
130		}
131		//! validate object
132		void validate() {
133		StateSpace<sq_T>::validate();
134		dimy = this->C.rows();
135		dimc = this->B.cols();
136		set_dim(this->A.rows());
137
138		bdm_assert(est.dimension(), "Statistics and model parameters mismatch");
139		}
	105	class Kalman: public BM, public StateSpace<sq_T> {
	106	protected:
	107	//! id of output
	108	RV yrv;
	109	//! Kalman gain
	110	mat _K;
	111	//!posterior
	112	enorm<sq_T> est;
	113	//!marginal on data f(y\|y)
	114	enorm<sq_T> fy;
	115	public:
	116	Kalman<sq_T>() : BM(), StateSpace<sq_T>(), yrv(), _K(), est() {}
	117	//! Copy constructor
	118	Kalman<sq_T> ( const Kalman<sq_T> &K0 ) : BM ( K0 ), StateSpace<sq_T> ( K0 ), yrv ( K0.yrv ), _K ( K0._K ), est ( K0.est ), fy ( K0.fy ) {}
	119	//!set statistics of the posterior
	120	void set_statistics ( const vec &mu0, const mat &P0 ) {
	121	est.set_parameters ( mu0, P0 );
	122	};
	123	//!set statistics of the posterior
	124	void set_statistics ( const vec &mu0, const sq_T &P0 ) {
	125	est.set_parameters ( mu0, P0 );
	126	};
	127	//! return correctly typed posterior (covariant return)
	128	const enorm<sq_T>& posterior() const {
	129	return est;
	130	}
	131	//! load basic elements of Kalman from structure
	132	void from_setting ( const Setting &set ) {
	133	StateSpace<sq_T>::from_setting ( set );
	134
	135	mat P0;
	136	vec mu0;
	137	UI::get ( mu0, set, "mu0", UI::optional );
	138	UI::get ( P0, set, "P0", UI::optional );
	139	set_statistics ( mu0, P0 );
	140	// Initial values
	141	UI::get ( yrv, set, "yrv", UI::optional );
	142	UI::get ( rvc, set, "urv", UI::optional );
	143	set_yrv ( concat ( yrv, rvc ) );
	144
	145	validate();
	146	}
	147	//! validate object
	148	void validate() {
	149	StateSpace<sq_T>::validate();
	150	dimy = this->C.rows();
	151	dimc = this->B.cols();
	152	set_dim ( this->A.rows() );
	153
	154	bdm_assert ( est.dimension(), "Statistics and model parameters mismatch" );
	155	}
140	156	};
141	157	/*!
…	…
143	159	*/
144	160
145		class KalmanFull : public Kalman<fsqmat>
146		{
147		public:
148		//! For EKFfull;
149		KalmanFull() :Kalman<fsqmat>(){};
150		//! Here dt = [yt;ut] of appropriate dimensions
151		void bayes (const vec &yt, const vec &cond=empty_vec);
152		BM* _copy_() const {
153		KalmanFull* K = new KalmanFull;
154		K->set_parameters (A, B, C, D, Q, R);
155		K->set_statistics (est._mu(), est._R());
156		return K;
157		}
158		};
159		UIREGISTER(KalmanFull);
	161	class KalmanFull : public Kalman<fsqmat> {
	162	public:
	163	//! For EKFfull;
	164	KalmanFull() : Kalman<fsqmat>() {};
	165	//! Here dt = [yt;ut] of appropriate dimensions
	166	void bayes ( const vec &yt, const vec &cond = empty_vec );
	167	BM* _copy_() const {
	168	KalmanFull* K = new KalmanFull;
	169	K->set_parameters ( A, B, C, D, Q, R );
	170	K->set_statistics ( est._mu(), est._R() );
	171	return K;
	172	}
	173	};
	174	UIREGISTER ( KalmanFull );
160	175
161	176
…	…
167	182	Complete constructor:
168	183	*/
169		class KalmanCh : public Kalman<chmat>
170		{
171		protected:
172		//! @{ \name Internal storage - needs initialize()
173		//! pre array (triangular matrix)
174		mat preA;
175		//! post array (triangular matrix)
176		mat postA;
177		//!@}
178		public:
179		//! copy constructor
180		BM* _copy_() const {
181		KalmanCh* K = new KalmanCh;
182		K->set_parameters (A, B, C, D, Q, R);
183		K->set_statistics (est._mu(), est._R());
184		K->validate();
185		return K;
186		}
187		//! set parameters for adapt from Kalman
188		void set_parameters (const mat &A0, const mat &B0, const mat &C0, const mat &D0, const chmat &Q0, const chmat &R0);
189		//! initialize internal parametetrs
190		void initialize();
191
192		/*!\brief Here dt = [yt;ut] of appropriate dimensions
193
194		The following equality hold::\f[
195		\left[\begin{array}{cc}
196		R^{0.5}\\
197		P_{t\|t-1}^{0.5}C' & P_{t\|t-1}^{0.5}CA'\\
198		& Q^{0.5}\end{array}\right]<\mathrm{orth.oper.}>=\left[\begin{array}{cc}
199		R_{y}^{0.5} & KA'\\
200		& P_{t+1\|t}^{0.5}\\
201		\\\end{array}\right]\f]
202
203		Thus this object evaluates only predictors! Not filtering densities.
204		*/
205		void bayes (const vec &yt, const vec &cond=empty_vec);
206
207		void from_setting(const Setting &set){
208		Kalman<chmat>::from_setting(set);
209		validate();
210		}
211		void validate() {
212		Kalman<chmat>::validate();
213		initialize();
214		}
215		};
216		UIREGISTER(KalmanCh);
	184	class KalmanCh : public Kalman<chmat> {
	185	protected:
	186	//! @{ \name Internal storage - needs initialize()
	187	//! pre array (triangular matrix)
	188	mat preA;
	189	//! post array (triangular matrix)
	190	mat postA;
	191	//!@}
	192	public:
	193	//! copy constructor
	194	BM* _copy_() const {
	195	KalmanCh* K = new KalmanCh;
	196	K->set_parameters ( A, B, C, D, Q, R );
	197	K->set_statistics ( est._mu(), est._R() );
	198	K->validate();
	199	return K;
	200	}
	201	//! set parameters for adapt from Kalman
	202	void set_parameters ( const mat &A0, const mat &B0, const mat &C0, const mat &D0, const chmat &Q0, const chmat &R0 );
	203	//! initialize internal parametetrs
	204	void initialize();
	205
	206	/*!\brief Here dt = [yt;ut] of appropriate dimensions
	207
	208	The following equality hold::\f[
	209	\left[\begin{array}{cc}
	210	R^{0.5}\\
	211	P_{t\|t-1}^{0.5}C' & P_{t\|t-1}^{0.5}CA'\\
	212	& Q^{0.5}\end{array}\right]<\mathrm{orth.oper.}>=\left[\begin{array}{cc}
	213	R_{y}^{0.5} & KA'\\
	214	& P_{t+1\|t}^{0.5}\\
	215	\\\end{array}\right]\f]
	216
	217	Thus this object evaluates only predictors! Not filtering densities.
	218	*/
	219	void bayes ( const vec &yt, const vec &cond = empty_vec );
	220
	221	void from_setting ( const Setting &set ) {
	222	Kalman<chmat>::from_setting ( set );
	223	validate();
	224	}
	225	void validate() {
	226	Kalman<chmat>::validate();
	227	initialize();
	228	}
	229	};
	230	UIREGISTER ( KalmanCh );
217	231
218	232	/*!
…	…
221	235	An approximation of the exact Bayesian filter with Gaussian noices and non-linear evolutions of their mean.
222	236	*/
223		class EKFfull : public KalmanFull
224		{
225		protected:
226		//! Internal Model f(x,u)
227		shared_ptr<diffbifn> pfxu;
228
229		//! Observation Model h(x,u)
230		shared_ptr<diffbifn> phxu;
231
232		public:
233		//! Default constructor
234		EKFfull ();
235
236		//! Set nonlinear functions for mean values and covariance matrices.
237		void set_parameters (const shared_ptr<diffbifn> &pfxu, const shared_ptr<diffbifn> &phxu, const mat Q0, const mat R0);
238
239		//! Here dt = [yt;ut] of appropriate dimensions
240		void bayes (const vec &yt, const vec &cond=empty_vec);
241		//! set estimates
242		void set_statistics (const vec &mu0, const mat &P0) {
243		est.set_parameters (mu0, P0);
244		};
245		//! access function
246		const mat _R() {
247		return est._R().to_mat();
248		}
249		void from_setting (const Setting &set) {
250		BM::from_setting(set);
251		shared_ptr<diffbifn> IM = UI::build<diffbifn> ( set, "IM", UI::compulsory );
252		shared_ptr<diffbifn> OM = UI::build<diffbifn> ( set, "OM", UI::compulsory );
253
254		//statistics
255		int dim = IM->dimension();
256		vec mu0;
257		if ( !UI::get ( mu0, set, "mu0" ) )
258		mu0 = zeros ( dim );
259
260		mat P0;
261		vec dP0;
262		if ( UI::get ( dP0, set, "dP0" ) )
263		P0 = diag ( dP0 );
264		else if ( !UI::get ( P0, set, "P0" ) )
265		P0 = eye ( dim );
266
267		set_statistics ( mu0, P0 );
268
269		//parameters
270		vec dQ, dR;
271		UI::get ( dQ, set, "dQ", UI::compulsory );
272		UI::get ( dR, set, "dR", UI::compulsory );
273		set_parameters ( IM, OM, diag ( dQ ), diag ( dR ) );
274
275		string options;
276		if ( UI::get ( options, set, "options" ) )
277		set_options ( options );
	237	class EKFfull : public KalmanFull {
	238	protected:
	239	//! Internal Model f(x,u)
	240	shared_ptr<diffbifn> pfxu;
	241
	242	//! Observation Model h(x,u)
	243	shared_ptr<diffbifn> phxu;
	244
	245	public:
	246	//! Default constructor
	247	EKFfull ();
	248
	249	//! Set nonlinear functions for mean values and covariance matrices.
	250	void set_parameters ( const shared_ptr<diffbifn> &pfxu, const shared_ptr<diffbifn> &phxu, const mat Q0, const mat R0 );
	251
	252	//! Here dt = [yt;ut] of appropriate dimensions
	253	void bayes ( const vec &yt, const vec &cond = empty_vec );
	254	//! set estimates
	255	void set_statistics ( const vec &mu0, const mat &P0 ) {
	256	est.set_parameters ( mu0, P0 );
	257	};
	258	//! access function
	259	const mat _R() {
	260	return est._R().to_mat();
	261	}
	262	void from_setting ( const Setting &set ) {
	263	BM::from_setting ( set );
	264	shared_ptr<diffbifn> IM = UI::build<diffbifn> ( set, "IM", UI::compulsory );
	265	shared_ptr<diffbifn> OM = UI::build<diffbifn> ( set, "OM", UI::compulsory );
	266
	267	//statistics
	268	int dim = IM->dimension();
	269	vec mu0;
	270	if ( !UI::get ( mu0, set, "mu0" ) )
	271	mu0 = zeros ( dim );
	272
	273	mat P0;
	274	vec dP0;
	275	if ( UI::get ( dP0, set, "dP0" ) )
	276	P0 = diag ( dP0 );
	277	else if ( !UI::get ( P0, set, "P0" ) )
	278	P0 = eye ( dim );
	279
	280	set_statistics ( mu0, P0 );
	281
	282	//parameters
	283	vec dQ, dR;
	284	UI::get ( dQ, set, "dQ", UI::compulsory );
	285	UI::get ( dR, set, "dR", UI::compulsory );
	286	set_parameters ( IM, OM, diag ( dQ ), diag ( dR ) );
	287
	288	string options;
	289	if ( UI::get ( options, set, "options" ) )
	290	set_options ( options );
278	291	// pfxu = UI::build<diffbifn>(set, "IM", UI::compulsory);
279	292	// phxu = UI::build<diffbifn>(set, "OM", UI::compulsory);
280		//
	293	//
281	294	// mat R0;
282	295	// UI::get(R0, set, "R",UI::compulsory);
283	296	// mat Q0;
284	297	// UI::get(Q0, set, "Q",UI::compulsory);
285		//
286		//
	298	//
	299	//
287	300	// mat P0; vec mu0;
288	301	// UI::get(mu0, set, "mu0", UI::optional);
…	…
293	306	// UI::get (urv, set, "urv", UI::optional);
294	307	// set_drv(concat(yrv,urv));
295		//
	308	//
296	309	// // setup StateSpace
297	310	// pfxu->dfdu_cond(mu0, zeros(pfxu->_dimu()), A,true);
298	311	// phxu->dfdu_cond(mu0, zeros(pfxu->_dimu()), C,true);
299		//
300		validate();
301		}
302		void validate() {
303		// check stats and IM and OM
304		}
305		};
306		UIREGISTER~~(EKFfull~~);
	312	//
	313	validate();
	314	}
	315	void validate() {
	316	// check stats and IM and OM
	317	}
	318	};
	319	UIREGISTER ( EKFfull );
307	320
308	321
…	…
313	326	*/
314	327
315		class EKFCh : public KalmanCh
316		{
317		protected:
318		//! Internal Model f(x,u)
319		shared_ptr<diffbifn> pfxu;
320
321		//! Observation Model h(x,u)
322		shared_ptr<diffbifn> phxu;
323		public:
324		//! copy constructor duplicated - calls different set_parameters
325		BM* _copy_() const {
326		EKFCh* E = new EKFCh;
327		E->set_parameters (pfxu, phxu, Q, R);
328		E->set_statistics (est._mu(), est._R());
329		return E;
330		}
331		//! Set nonlinear functions for mean values and covariance matrices.
332		void set_parameters (const shared_ptr<diffbifn> &pfxu, const shared_ptr<diffbifn> &phxu, const chmat Q0, const chmat R0);
333
334		//! Here dt = [yt;ut] of appropriate dimensions
335		void bayes (const vec &yt, const vec &cond=empty_vec);
336
337		void from_setting (const Setting &set);
338
339		void validate(){};
340		// TODO dodelat void to_setting( Setting &set ) const;
341
342		};
343
344		UIREGISTER (EKFCh);
345		SHAREDPTR (EKFCh);
	328	class EKFCh : public KalmanCh {
	329	protected:
	330	//! Internal Model f(x,u)
	331	shared_ptr<diffbifn> pfxu;
	332
	333	//! Observation Model h(x,u)
	334	shared_ptr<diffbifn> phxu;
	335	public:
	336	//! copy constructor duplicated - calls different set_parameters
	337	BM* _copy_() const {
	338	EKFCh* E = new EKFCh;
	339	E->set_parameters ( pfxu, phxu, Q, R );
	340	E->set_statistics ( est._mu(), est._R() );
	341	return E;
	342	}
	343	//! Set nonlinear functions for mean values and covariance matrices.
	344	void set_parameters ( const shared_ptr<diffbifn> &pfxu, const shared_ptr<diffbifn> &phxu, const chmat Q0, const chmat R0 );
	345
	346	//! Here dt = [yt;ut] of appropriate dimensions
	347	void bayes ( const vec &yt, const vec &cond = empty_vec );
	348
	349	void from_setting ( const Setting &set );
	350
	351	void validate() {};
	352	// TODO dodelat void to_setting( Setting &set ) const;
	353
	354	};
	355
	356	UIREGISTER ( EKFCh );
	357	SHAREDPTR ( EKFCh );
346	358
347	359
…	…
355	367	The next step is performed with the new statistics for all models.
356	368	*/
357		class MultiModel: public BM
358		{
359		protected:
360		//! List of models between which we switch
361		Array<EKFCh*> Models;
362		//! vector of model weights
363		vec w;
364		//! cache of model lls
365		vec _lls;
366		//! type of switching policy [1=maximum,2=...]
367		int policy;
368		//! internal statistics
369		enorm<chmat> est;
370		public:
371		//! set internal parameters
372		void set_parameters (Array<EKFCh*> A, int pol0 = 1) {
373		Models = A;//TODO: test if evalll is set
374		w.set_length (A.length());
375		_lls.set_length (A.length());
376		policy = pol0;
377
378		est.set_rv (RV ("MM", A (0)->posterior().dimension(), 0));
379		est.set_parameters (A (0)->posterior().mean(), A (0)->posterior()._R());
380		}
381		void bayes (const vec &yt, const vec &cond=empty_vec) {
382		int n = Models.length();
383		int i;
384		for (i = 0; i < n; i++) {
385		Models (i)->bayes (yt);
386		_lls (i) = Models (i)->_ll();
387		}
388		double mlls = max (_lls);
389		w = exp (_lls - mlls);
390		w /= sum (w); //normalization
391		//set statistics
392		switch (policy) {
393		case 1: {
394		int mi = max_index (w);
395		const enorm<chmat> &st = Models (mi)->posterior() ;
396		est.set_parameters (st.mean(), st._R());
397		}
398		break;
399		default:
400		bdm_error ("unknown policy");
401		}
402		// copy result to all models
403		for (i = 0; i < n; i++) {
404		Models (i)->set_statistics (est.mean(), est._R());
405		}
406		}
407		//! return correctly typed posterior (covariant return)
408		const enorm<chmat>& posterior() const {
409		return est;
410		}
411
412		void from_setting (const Setting &set);
413
414		// TODO dodelat void to_setting( Setting &set ) const;
415
416		};
417
418		UIREGISTER (MultiModel);
419		SHAREDPTR (MultiModel);
420
421		//! conversion of outer ARX model (mlnorm) to state space model
	369	class MultiModel: public BM {
	370	protected:
	371	//! List of models between which we switch
	372	Array<EKFCh*> Models;
	373	//! vector of model weights
	374	vec w;
	375	//! cache of model lls
	376	vec _lls;
	377	//! type of switching policy [1=maximum,2=...]
	378	int policy;
	379	//! internal statistics
	380	enorm<chmat> est;
	381	public:
	382	//! set internal parameters
	383	void set_parameters ( Array<EKFCh*> A, int pol0 = 1 ) {
	384	Models = A;//TODO: test if evalll is set
	385	w.set_length ( A.length() );
	386	_lls.set_length ( A.length() );
	387	policy = pol0;
	388
	389	est.set_rv ( RV ( "MM", A ( 0 )->posterior().dimension(), 0 ) );
	390	est.set_parameters ( A ( 0 )->posterior().mean(), A ( 0 )->posterior()._R() );
	391	}
	392	void bayes ( const vec &yt, const vec &cond = empty_vec ) {
	393	int n = Models.length();
	394	int i;
	395	for ( i = 0; i < n; i++ ) {
	396	Models ( i )->bayes ( yt );
	397	_lls ( i ) = Models ( i )->_ll();
	398	}
	399	double mlls = max ( _lls );
	400	w = exp ( _lls - mlls );
	401	w /= sum ( w ); //normalization
	402	//set statistics
	403	switch ( policy ) {
	404	case 1: {
	405	int mi = max_index ( w );
	406	const enorm<chmat> &st = Models ( mi )->posterior() ;
	407	est.set_parameters ( st.mean(), st._R() );
	408	}
	409	break;
	410	default:
	411	bdm_error ( "unknown policy" );
	412	}
	413	// copy result to all models
	414	for ( i = 0; i < n; i++ ) {
	415	Models ( i )->set_statistics ( est.mean(), est._R() );
	416	}
	417	}
	418	//! return correctly typed posterior (covariant return)
	419	const enorm<chmat>& posterior() const {
	420	return est;
	421	}
	422
	423	void from_setting ( const Setting &set );
	424
	425	// TODO dodelat void to_setting( Setting &set ) const;
	426
	427	};
	428
	429	UIREGISTER ( MultiModel );
	430	SHAREDPTR ( MultiModel );
	431
	432	//! conversion of outer ARX model (mlnorm) to state space model
422	433	/*!
423	434	The model is constructed as:
…	…
429	440	*/
430	441	//template<class sq_T>
431		class StateCanonical: public StateSpace<fsqmat>{
432		protected:
433		//! remember connection from theta ->A
434		datalink_part th2A;
435		//! remember connection from theta ->C
436		datalink_part th2C;
437		//! remember connection from theta ->D
438		datalink_part th2D;
439		//!cached first row of A
440		vec A1row;
441		//!cached first row of C
442		vec C1row;
443		//!cached first row of D
444		vec D1row;
445
446		public:
447		//! set up this object to match given mlnorm
448		void connect_mlnorm(const mlnorm<fsqmat> &ml){
449		//get ids of yrv
450		const RV &yrv = ml._rv();
451		//need to determine u_t - it is all in _rvc that is not in ml._rv()
452		RV rgr0 = ml._rvc().remove_time();
453		RV urv = rgr0.subt(yrv);
454
455		//We can do only 1d now... :(
456		bdm_assert(yrv._dsize()==1, "Only for SISO so far..." );
457
458		// create names for
459		RV xrv; //empty
460		RV Crv; //empty
461		int td=ml._rvc().mint();
462		// assuming strictly proper function!!!
463		for (int t=-1;t>=td;t--){
464		xrv.add(yrv.copy_t(t));
465		Crv.add(urv.copy_t(t));
466		}
467
468		// get mapp
469		th2A.set_connection(xrv, ml._rvc());
470		th2C.set_connection(Crv, ml._rvc());
471		th2D.set_connection(urv, ml._rvc());
472
473		//set matrix sizes
474		this->A=zeros(xrv._dsize(),xrv._dsize());
475		for (int j=1; j<xrv._dsize(); j++){A(j,j-1)=1.0;} // off diagonal
476		this->B=zeros(xrv._dsize(),1);
477		this->B(0) = 1.0;
478		this->C=zeros(1,xrv._dsize());
479		this->D=zeros(1,urv._dsize());
480		this->Q = zeros(xrv._dsize(),xrv._dsize());
481		// R is set by update
482
483		//set cache
484		this->A1row = zeros(xrv._dsize());
485		this->C1row = zeros(xrv._dsize());
486		this->D1row = zeros(urv._dsize());
487
488		update_from(ml);
489		validate();
490		};
491		//! fast function to update parameters from ml - not checked for compatibility!!
492		void update_from(const mlnorm<fsqmat> &ml){
493
494		vec theta = ml._A().get_row(0); // this
495
496		th2A.filldown(theta,A1row);
497		th2C.filldown(theta,C1row);
498		th2D.filldown(theta,D1row);
499
500		R = ml._R();
501
502		A.set_row(0,A1row);
503		C.set_row(0,C1row+D1row(0)*A1row);
504		D.set_row(0,D1row);
505
506		}
	442	class StateCanonical: public StateSpace<fsqmat> {
	443	protected:
	444	//! remember connection from theta ->A
	445	datalink_part th2A;
	446	//! remember connection from theta ->C
	447	datalink_part th2C;
	448	//! remember connection from theta ->D
	449	datalink_part th2D;
	450	//!cached first row of A
	451	vec A1row;
	452	//!cached first row of C
	453	vec C1row;
	454	//!cached first row of D
	455	vec D1row;
	456
	457	public:
	458	//! set up this object to match given mlnorm
	459	void connect_mlnorm ( const mlnorm<fsqmat> &ml ) {
	460	//get ids of yrv
	461	const RV &yrv = ml._rv();
	462	//need to determine u_t - it is all in _rvc that is not in ml._rv()
	463	RV rgr0 = ml._rvc().remove_time();
	464	RV urv = rgr0.subt ( yrv );
	465
	466	//We can do only 1d now... :(
	467	bdm_assert ( yrv._dsize() == 1, "Only for SISO so far..." );
	468
	469	// create names for
	470	RV xrv; //empty
	471	RV Crv; //empty
	472	int td = ml._rvc().mint();
	473	// assuming strictly proper function!!!
	474	for ( int t = -1; t >= td; t-- ) {
	475	xrv.add ( yrv.copy_t ( t ) );
	476	Crv.add ( urv.copy_t ( t ) );
	477	}
	478
	479	// get mapp
	480	th2A.set_connection ( xrv, ml._rvc() );
	481	th2C.set_connection ( Crv, ml._rvc() );
	482	th2D.set_connection ( urv, ml._rvc() );
	483
	484	//set matrix sizes
	485	this->A = zeros ( xrv._dsize(), xrv._dsize() );
	486	for ( int j = 1; j < xrv._dsize(); j++ ) {
	487	A ( j, j - 1 ) = 1.0; // off diagonal
	488	}
	489	this->B = zeros ( xrv._dsize(), 1 );
	490	this->B ( 0 ) = 1.0;
	491	this->C = zeros ( 1, xrv._dsize() );
	492	this->D = zeros ( 1, urv._dsize() );
	493	this->Q = zeros ( xrv._dsize(), xrv._dsize() );
	494	// R is set by update
	495
	496	//set cache
	497	this->A1row = zeros ( xrv._dsize() );
	498	this->C1row = zeros ( xrv._dsize() );
	499	this->D1row = zeros ( urv._dsize() );
	500
	501	update_from ( ml );
	502	validate();
	503	};
	504	//! fast function to update parameters from ml - not checked for compatibility!!
	505	void update_from ( const mlnorm<fsqmat> &ml ) {
	506
	507	vec theta = ml._A().get_row ( 0 ); // this
	508
	509	th2A.filldown ( theta, A1row );
	510	th2C.filldown ( theta, C1row );
	511	th2D.filldown ( theta, D1row );
	512
	513	R = ml._R();
	514
	515	A.set_row ( 0, A1row );
	516	C.set_row ( 0, C1row + D1row ( 0 ) *A1row );
	517	D.set_row ( 0, D1row );
	518
	519	}
507	520	};
508	521	/*!
509	522	State-Space representation of multivariate autoregressive model.
510	523	The original model:
511		\f[ y_t = \theta [\ldots y_{t-k}, \ldots u_{t-l}, \ldots z_{t-m}]' + \Sigma^{-1/2} e_t \f]
	524	\f[ y_t = \theta [\ldots y_{t-k}, \ldots u_{t-l}, \ldots z_{t-m}]' + \Sigma^{-1/2} e_t \f]
512	525	where \f$ k,l,m \f$ are maximum delayes of corresponding variables in the regressor.
513	526
…	…
519	532
520	533	*/
521		class StateFromARX: public StateSpace<chmat>{
522		protected:
523		//! remember connection from theta ->A
524		datalink_part th2A;
525		//! remember connection from theta ->B
526		datalink_part th2B;
527		//!function adds n diagonal elements from given starting point r,c
528		void diagonal_part(mat &A, int r, int c, int n){
529		for(int i=0; i<n; i++){A(r,c)=1.0; r++;c++;}
530		};
531		//! similar to ARX.have_constant
532		bool have_constant;
533		public:
534		//! set up this object to match given mlnorm
535		//! Note that state-space and common mpdf use different meaning of \f$ _t \f$ in \f$ u_t \f$.
536		//!While mlnorm typically assumes that \f$ u_t \rightarrow y_t \f$ in state space it is \f$ u_{t-1} \rightarrow y_t \f$
537		//! For consequences in notation of internal variable xt see arx2statespace_notes.lyx.
538		void connect_mlnorm(const mlnorm<chmat> &ml, RV &xrv, RV &urv){
539
540		//get ids of yrv
541		const RV &yrv = ml._rv();
542		//need to determine u_t - it is all in _rvc that is not in ml._rv()
543		const RV &rgr = ml._rvc();
544		RV rgr0 = rgr.remove_time();
545		urv = rgr0.subt(yrv);
546
547		// create names for state variables
548		xrv=yrv;
549
550		int y_multiplicity = -rgr.mint(yrv);
551		int y_block_size = yrv.length()*(y_multiplicity); // current yt + all delayed yts
552		for (int m=0;m<y_multiplicity - 1;m++){ // ========= -1 is important see arx2statespace_notes
553		xrv.add(yrv.copy_t(-m-1)); //add delayed yt
	534	class StateFromARX: public StateSpace<chmat> {
	535	protected:
	536	//! remember connection from theta ->A
	537	datalink_part th2A;
	538	//! remember connection from theta ->B
	539	datalink_part th2B;
	540	//!function adds n diagonal elements from given starting point r,c
	541	void diagonal_part ( mat &A, int r, int c, int n ) {
	542	for ( int i = 0; i < n; i++ ) {
	543	A ( r, c ) = 1.0;
	544	r++;
	545	c++;
	546	}
	547	};
	548	//! similar to ARX.have_constant
	549	bool have_constant;
	550	public:
	551	//! set up this object to match given mlnorm
	552	//! Note that state-space and common mpdf use different meaning of \f$ _t \f$ in \f$ u_t \f$.
	553	//!While mlnorm typically assumes that \f$ u_t \rightarrow y_t \f$ in state space it is \f$ u_{t-1} \rightarrow y_t \f$
	554	//! For consequences in notation of internal variable xt see arx2statespace_notes.lyx.
	555	void connect_mlnorm ( const mlnorm<chmat> &ml, RV &xrv, RV &urv ) {
	556
	557	//get ids of yrv
	558	const RV &yrv = ml._rv();
	559	//need to determine u_t - it is all in _rvc that is not in ml._rv()
	560	const RV &rgr = ml._rvc();
	561	RV rgr0 = rgr.remove_time();
	562	urv = rgr0.subt ( yrv );
	563
	564	// create names for state variables
	565	xrv = yrv;
	566
	567	int y_multiplicity = -rgr.mint ( yrv );
	568	int y_block_size = yrv.length() * ( y_multiplicity ); // current yt + all delayed yts
	569	for ( int m = 0; m < y_multiplicity - 1; m++ ) { // ========= -1 is important see arx2statespace_notes
	570	xrv.add ( yrv.copy_t ( -m - 1 ) ); //add delayed yt
	571	}
	572	//! temporary RV for connection to ml.rvc, since notation of xrv and ml.rvc does not match
	573	RV xrv_ml = xrv.copy_t ( -1 );
	574
	575	// add regressors
	576	ivec u_block_sizes ( urv.length() ); // no_blocks = yt + unique rgr
	577	for ( int r = 0; r < urv.length(); r++ ) {
	578	RV R = urv.subselect ( vec_1 ( r ) ); //non-delayed element of rgr
	579	int r_size = urv.size ( r );
	580	int r_multiplicity = -rgr.mint ( R );
	581	u_block_sizes ( r ) = r_size * r_multiplicity ;
	582	for ( int m = 0; m < r_multiplicity; m++ ) {
	583	xrv.add ( R.copy_t ( -m - 1 ) ); //add delayed yt
	584	xrv_ml.add ( R.copy_t ( -m - 1 ) ); //add delayed yt
554	585	}
555		//! temporary RV for connection to ml.rvc, since notation of xrv and ml.rvc does not match
556		RV xrv_ml = xrv.copy_t(-1);
557
558		// add regressors
559		ivec u_block_sizes(urv.length()); // no_blocks = yt + unique rgr
560		for (int r=0;r<urv.length();r++){
561		RV R=urv.subselect(vec_1(r)); //non-delayed element of rgr
562		int r_size=urv.size(r);
563		int r_multiplicity=-rgr.mint(R);
564		u_block_sizes(r)= r_size*r_multiplicity ;
565		for(int m=0;m<r_multiplicity;m++){
566		xrv.add(R.copy_t(-m-1)); //add delayed yt
567		xrv_ml.add(R.copy_t(-m-1)); //add delayed yt
568		}
	586	}
	587	// add constant
	588	if ( any ( ml._mu_const() != 0.0 ) ) {
	589	have_constant = true;
	590	xrv.add ( RV ( "bdm_reserved_constant_one", 1 ) );
	591	} else {
	592	have_constant = false;
	593	}
	594
	595
	596	// get mapp
	597	th2A.set_connection ( xrv_ml, ml._rvc() );
	598	th2B.set_connection ( urv, ml._rvc() );
	599
	600	//set matrix sizes
	601	this->A = zeros ( xrv._dsize(), xrv._dsize() );
	602	//create y block
	603	diagonal_part ( this->A, yrv._dsize(), 0, y_block_size - yrv._dsize() );
	604
	605	this->B = zeros ( xrv._dsize(), urv._dsize() );
	606	//add diagonals for rgr
	607	int active_x = y_block_size;
	608	for ( int r = 0; r < urv.length(); r++ ) {
	609	diagonal_part ( this->A, active_x + urv.size ( r ), active_x, u_block_sizes ( r ) - urv.size ( r ) );
	610	this->B.set_submatrix ( active_x, 0, eye ( urv.size ( r ) ) );
	611	active_x += u_block_sizes ( r );
	612	}
	613	this->C = zeros ( yrv._dsize(), xrv._dsize() );
	614	this->C.set_submatrix ( 0, 0, eye ( yrv._dsize() ) );
	615	this->D = zeros ( yrv._dsize(), urv._dsize() );
	616	this->R.setCh ( zeros ( yrv._dsize(), yrv._dsize() ) );
	617	this->Q.setCh ( zeros ( xrv._dsize(), xrv._dsize() ) );
	618	// Q is set by update
	619
	620	update_from ( ml );
	621	validate();
	622	};
	623	//! fast function to update parameters from ml - not checked for compatibility!!
	624	void update_from ( const mlnorm<chmat> &ml ) {
	625
	626	vec Arow = zeros ( A.cols() );
	627	vec Brow = zeros ( B.cols() );
	628	// ROW- WISE EVALUATION =====
	629	for ( int i = 0; i < ml._rv()._dsize(); i++ ) {
	630
	631	vec theta = ml._A().get_row ( i );
	632
	633	th2A.filldown ( theta, Arow );
	634	if ( have_constant ) {
	635	// constant is always at the end no need for datalink
	636	Arow ( A.cols() - 1 ) = ml._mu_const() ( i );
569	637	}
570		// add constant
571		if (any(ml._mu_const()!=0.0)){
572		have_constant=true;
573		xrv.add(RV("bdm_reserved_constant_one",1));
574		} else {have_constant=false;}
575
576
577		// get mapp
578		th2A.set_connection(xrv_ml, ml._rvc());
579		th2B.set_connection(urv, ml._rvc());
580
581		//set matrix sizes
582		this->A=zeros(xrv._dsize(),xrv._dsize());
583		//create y block
584		diagonal_part(this->A, yrv._dsize(), 0, y_block_size-yrv._dsize() );
585
586		this->B=zeros(xrv._dsize(), urv._dsize());
587		//add diagonals for rgr
588		int active_x=y_block_size;
589		for (int r=0; r<urv.length(); r++){
590		diagonal_part( this->A, active_x+urv.size(r),active_x, u_block_sizes(r)-urv.size(r));
591		this->B.set_submatrix(active_x, 0, eye(urv.size(r)));
592		active_x+=u_block_sizes(r);
593		}
594		this->C=zeros(yrv._dsize(), xrv._dsize());
595		this->C.set_submatrix(0,0,eye(yrv._dsize()));
596		this->D=zeros(yrv._dsize(),urv._dsize());
597		this->R.setCh(zeros(yrv._dsize(),yrv._dsize()));
598		this->Q.setCh(zeros(xrv._dsize(),xrv._dsize()));
599		// Q is set by update
600
601		update_from(ml);
602		validate();
603		};
604		//! fast function to update parameters from ml - not checked for compatibility!!
605		void update_from(const mlnorm<chmat> &ml){
606
607		vec Arow=zeros(A.cols());
608		vec Brow=zeros(B.cols());
609		// ROW- WISE EVALUATION =====
610		for(int i=0;i<ml._rv()._dsize(); i++){
611
612		vec theta = ml._A().get_row(i);
613
614		th2A.filldown(theta,Arow);
615		if (have_constant){
616		// constant is always at the end no need for datalink
617		Arow(A.cols()-1) = ml._mu_const()(i);
618		}
619		this->A.set_row(i,Arow);
620
621		th2B.filldown(theta,Brow);
622		this->B.set_row(i,Brow);
623		}
624		this->Q._Ch().set_submatrix(0,0,ml.__R()._Ch());
625
626		};
627		//! access function
628		bool _have_constant() const {return have_constant;}
	638	this->A.set_row ( i, Arow );
	639
	640	th2B.filldown ( theta, Brow );
	641	this->B.set_row ( i, Brow );
	642	}
	643	this->Q._Ch().set_submatrix ( 0, 0, ml.__R()._Ch() );
	644
	645	};
	646	//! access function
	647	bool _have_constant() const {
	648	return have_constant;
	649	}
629	650	};
630	651
…	…
632	653
633	654	template<class sq_T>
634		void StateSpace<sq_T>::set_parameters (const mat &A0, const mat &B0, const mat &C0, const mat &D0, const sq_T &Q0, const sq_T &R0)
635		{
	655	void StateSpace<sq_T>::set_parameters ( const mat &A0, const mat &B0, const mat &C0, const mat &D0, const sq_T &Q0, const sq_T &R0 ) {
636	656
637	657	A = A0;
…	…
645	665
646	666	template<class sq_T>
647		void StateSpace<sq_T>::validate(){
648		bdm_assert (~~A.cols() == A.rows(), "KalmanFull: A is not square"~~);
649		bdm_assert (~~B.rows() == A.rows(), "KalmanFull: B is not compatible"~~);
650		bdm_assert (~~C.cols() == A.rows(), "KalmanFull: C is not compatible"~~);
651		bdm_assert ( (~~D.rows() == C.rows()) && (D.cols() == B.cols()), "KalmanFull: D is not compatible"~~);
652		bdm_assert ( (~~Q.cols() == A.rows()) && (Q.rows() == A.rows()), "KalmanFull: Q is not compatible"~~);
653		bdm_assert ( (~~R.cols() == C.rows()) && (R.rows() == C.rows()), "KalmanFull: R is not compatible"~~);
	667	void StateSpace<sq_T>::validate() {
	668	bdm_assert ( A.cols() == A.rows(), "KalmanFull: A is not square" );
	669	bdm_assert ( B.rows() == A.rows(), "KalmanFull: B is not compatible" );
	670	bdm_assert ( C.cols() == A.rows(), "KalmanFull: C is not compatible" );
	671	bdm_assert ( ( D.rows() == C.rows() ) && ( D.cols() == B.cols() ), "KalmanFull: D is not compatible" );
	672	bdm_assert ( ( Q.cols() == A.rows() ) && ( Q.rows() == A.rows() ), "KalmanFull: Q is not compatible" );
	673	bdm_assert ( ( R.cols() == C.rows() ) && ( R.rows() == C.rows() ), "KalmanFull: R is not compatible" );
654	674	}
655	675

library/bdm/estim/mixtures.cpp

r735	r737
121	121	}
122	122
123		void MixEF::bayes ( const vec &data, const vec &cond=empty_vec ) {
	123	void MixEF::bayes ( const vec &data, const vec &cond = empty_vec ) {
124	124
125	125	};
126	126
127		void MixEF::bayes ( const mat &data, const vec &cond=empty_vec ) {
	127	void MixEF::bayes ( const mat &data, const vec &cond = empty_vec ) {
128	128	this->bayes_batch ( data, cond, ones ( data.cols() ) );
129	129	};

library/bdm/estim/mixtures.h

r735	r737
111	111	//! EM algorithm
112	112	void bayes ( const mat &yt, const vec &cond );
113		//! batch weighted Bayes rule
	113	//! batch weighted Bayes rule
114	114	void bayes_batch ( const mat &yt, const mat &cond, const vec &wData );
115	115	double logpred ( const vec &yt ) const;
…	…
130	130	method = M;
131	131	}
132
133		void to_setting~~(Setting &set) const~~{
134		UI::save~~(Coms,set,"Coms"~~);
135		Setting &wei~~=set.add("weights",Setting::TypeGroup~~);
136		weights.to_setting~~(wei~~);
	132
	133	void to_setting ( Setting &set ) const {
	134	UI::save ( Coms, set, "Coms" );
	135	Setting &wei = set.add ( "weights", Setting::TypeGroup );
	136	weights.to_setting ( wei );
137	137	}
138	138	};

library/bdm/estim/particles.cpp

r700	r737
5	5	using std::endl;
6	6
7		void PF::bayes_gensmp~~(const vec &ut)~~{
8		if (~~ut.length()>0)~~{
9		vec cond~~(par->dimensionc()~~);
10		cond.set_subvector~~(par->dimension(),ut~~);
11		#pragma parallel for
12		for (int i = 0; i < n; i++ ) {
13		cond.set_subvector~~(0,_samples ( i )~~);
	7	void PF::bayes_gensmp ( const vec &ut ) {
	8	if ( ut.length() > 0 ) {
	9	vec cond ( par->dimensionc() );
	10	cond.set_subvector ( par->dimension(), ut );
	11	#pragma parallel for
	12	for ( int i = 0; i < n; i++ ) {
	13	cond.set_subvector ( 0, _samples ( i ) );
14	14	_samples ( i ) = par->samplecond ( cond );
15		lls(i) = 0;
	15	lls ( i ) = 0;
16	16	}
17		} else {
18		#pragma parallel for
19		for (int i = 0; i < n; i++ ) {
	17	} else {
	18	#pragma parallel for
	19	for ( int i = 0; i < n; i++ ) {
20	20	_samples ( i ) = par->samplecond ( _samples ( i ) );
21		lls(i) = 0;
	21	lls ( i ) = 0;
22	22	}
23	23	}
24	24	}
25	25
26		void PF::bayes_weights(){
27		//
28		double mlls~~=max(lls~~);
	26	void PF::bayes_weights() {
	27	//
	28	double mlls = max ( lls );
29	29	// compute weights
30		for (int i = 0; i < n; i++ ) {
	30	for ( int i = 0; i < n; i++ ) {
31	31	_w ( i ) *= exp ( lls ( i ) - mlls ); // multiply w by likelihood
32	32	}
33	33
34	34	//renormalize
35		double sw=sum(_w);
36		if (!std::isfinite(sw)) {
37		for (int i=0;i<n;i++){
38		if (!std::isfinite(_w(i))) {_w(i)=0;}
	35	double sw = sum ( _w );
	36	if ( !std::isfinite ( sw ) ) {
	37	for ( int i = 0; i < n; i++ ) {
	38	if ( !std::isfinite ( _w ( i ) ) ) {
	39	_w ( i ) = 0;
	40	}
39	41	}
40		sw = sum~~(_w~~);
41		if (~~!std::isfinite(sw) \|\| sw==0.0~~) {
42		bdm_error~~("Particle filter is lost; no particle is good enough."~~);
	42	sw = sum ( _w );
	43	if ( !std::isfinite ( sw ) \|\| sw == 0.0 ) {
	44	bdm_error ( "Particle filter is lost; no particle is good enough." );
43	45	}
44	46	}
…	…
46	48	}
47	49
48		void PF::bayes( const vec &yt, const vec &cond ) {
49		const vec &ut=cond; //todo
50
	50	void PF::bayes ( const vec &yt, const vec &cond ) {
	51	const vec &ut = cond; //todo
	52
51	53	int i;
52	54	// generate samples - time step
53		bayes_gensmp~~(ut~~);
	55	bayes_gensmp ( ut );
54	56	// weight them - data step
55	57	for ( i = 0; i < n; i++ ) {
…	…
59	61	bayes_weights();
60	62
61		if (~~do_resampling()~~) {
62		est.resample ( resmethod);
	63	if ( do_resampling() ) {
	64	est.resample ( resmethod );
63	65	}
64	66
…	…
73	75	// }
74	76
75		void MPF::bayes ( const vec &yt, const vec &cond ) {
76		// follows PF::bayes in most places!!!
	77	void MPF::bayes ( const vec &yt, const vec &cond ) {
	78	// follows PF::bayes in most places!!!
77	79	int i;
78		int n=pf->__w().length();
	80	int n = pf->__w().length();
79	81	vec &lls = pf->_lls();
80		Array<vec> &samples=pf->__samples();
81
	82	Array<vec> &samples = pf->__samples();
	83
82	84	// generate samples - time step
83		pf->bayes_gensmp~~(vec(0)~~);
	85	pf->bayes_gensmp ( vec ( 0 ) );
84	86	// weight them - data step
85		#pragma parallel for
	87	#pragma parallel for
86	88	for ( i = 0; i < n; i++ ) {
87		vec bm_cond~~(BMs(i)->dimensionc()~~);
88		this2bm.fill_cond~~(yt,cond, bm_cond~~);
89		pf2bm.filldown~~(samples(i), bm_cond~~);
90		BMs~~(i) -> bayes(this2bm.pushdown(yt), bm_cond~~);
91		lls ( i ) += BMs(i)->_ll();
	89	vec bm_cond ( BMs ( i )->dimensionc() );
	90	this2bm.fill_cond ( yt, cond, bm_cond );
	91	pf2bm.filldown ( samples ( i ), bm_cond );
	92	BMs ( i ) -> bayes ( this2bm.pushdown ( yt ), bm_cond );
	93	lls ( i ) += BMs ( i )->_ll();
92	94	}
93	95
94	96	pf->bayes_weights();
95
	97
96	98	ivec ind;
97		if (~~pf->do_resampling()~~) {
98		pf->resample ( ind);
99
100		#pragma omp parallel for
	99	if ( pf->do_resampling() ) {
	100	pf->resample ( ind );
	101
	102	#pragma omp parallel for
101	103	for ( i = 0; i < n; i++ ) {
102	104	if ( ind ( i ) != i ) {//replace the current Bm by a new one

library/bdm/estim/particles.h

r700	r737
41	41	//! internal structure storing loglikelihood of predictions
42	42	vec lls;
43
	43
44	44	//! which resampling method will be used
45	45	RESAMPLING_METHOD resmethod;
…	…
47	47	//! For example, for 0.5 resampling is performed when the numebr of active aprticles drops belo 50%.
48	48	double res_threshold;
49
	49
50	50	//! \name Options
51	51	//!@{
…	…
57	57	PF ( ) : est(), _w ( est._w() ), _samples ( est._samples() ) {
58	58	};
59
60		void set_parameters (~~int n0, double res_th0=~~0.5, RESAMPLING_METHOD rm = SYSTEMATIC ) {
	59
	60	void set_parameters ( int n0, double res_th0 = 0.5, RESAMPLING_METHOD rm = SYSTEMATIC ) {
61	61	n = n0;
62	62	res_threshold = res_th0;
63	63	resmethod = rm;
64	64	};
65		void set_model ( shared_ptr<pdf> par0, shared_ptr<pdf> obs0) {
	65	void set_model ( shared_ptr<pdf> par0, shared_ptr<pdf> obs0 ) {
66	66	par = par0;
67	67	obs = obs0;
68	68	// set values for posterior
69		est.set_rv~~(par->_rv()~~);
	69	est.set_rv ( par->_rv() );
70	70	};
71	71	void set_statistics ( const vec w0, const epdf &epdf0 ) {
…	…
73	73	};
74	74	void set_statistics ( const eEmp &epdf0 ) {
75		bdm_assert_debug~~(epdf0._rv().equal(par->_rv()),"Incompatibel input"~~);
76		est=epdf0;
	75	bdm_assert_debug ( epdf0._rv().equal ( par->_rv() ), "Incompatibel input" );
	76	est = epdf0;
77	77	};
78	78	//!@}
…	…
85	85	}
86	86	//! bayes I - generate samples and add their weights to lls
87		virtual void bayes_gensmp~~(const vec &ut~~);
88		//! bayes II - compute weights of the
	87	virtual void bayes_gensmp ( const vec &ut );
	88	//! bayes II - compute weights of the
89	89	virtual void bayes_weights();
90	90	//! important part of particle filtering - decide if it is time to perform resampling
91		virtual bool do_resampling(){
	91	virtual bool do_resampling() {
92	92	double eff = 1.0 / ( _w * _w );
93	93	return eff < ( res_threshold*n );
…	…
95	95	void bayes ( const vec &yt, const vec &cond );
96	96	//!access function
97		vec& __w() { return _w; }
	97	vec& __w() {
	98	return _w;
	99	}
98	100	//!access function
99		vec& _lls() { return lls; }
	101	vec& _lls() {
	102	return lls;
	103	}
100	104	//!access function
101		RESAMPLING_METHOD _resmethod() const { return resmethod; }
	105	RESAMPLING_METHOD _resmethod() const {
	106	return resmethod;
	107	}
102	108	//! return correctly typed posterior (covariant return)
103		const eEmp& posterior() const {return est;}
104
	109	const eEmp& posterior() const {
	110	return est;
	111	}
	112
105	113	/*! configuration structure for basic PF
106	114	\code
…	…
115	123	\endcode
116	124	*/
117		void from_setting~~(const Setting &set)~~{
118		BM::from_setting~~(set~~);
119		par = UI::build<pdf>~~(set,"parameter_pdf",UI::compulsory~~);
120		obs = UI::build<pdf>~~(set,"observation_pdf",UI::compulsory~~);
121
122		prior_from_set~~(set~~);
123		resmethod_from_set~~(set~~);
	125	void from_setting ( const Setting &set ) {
	126	BM::from_setting ( set );
	127	par = UI::build<pdf> ( set, "parameter_pdf", UI::compulsory );
	128	obs = UI::build<pdf> ( set, "observation_pdf", UI::compulsory );
	129
	130	prior_from_set ( set );
	131	resmethod_from_set ( set );
124	132	// set resampling method
125	133	//set drv
126	134	//find potential input - what remains in rvc when we subtract rv
127		RV u = par->_rvc().remove_time().subt~~( par->_rv() );~~
	135	RV u = par->_rvc().remove_time().subt ( par->_rv() );
128	136	//find potential input - what remains in rvc when we subtract x_t
129		RV obs_u = obs->_rvc().remove_time().subt~~( par->_rv() );~~
130
131		u.add~~(obs_u~~); // join both u, and check if they do not overlap
132
133		set_yrv(obs->_rv() );
	137	RV obs_u = obs->_rvc().remove_time().subt ( par->_rv() );
	138
	139	u.add ( obs_u ); // join both u, and check if they do not overlap
	140
	141	set_yrv ( obs->_rv() );
134	142	rvc = u;
135	143	}
136	144	//! auxiliary function reading parameter 'resmethod' from configuration file
137		void resmethod_from_set~~(const Setting &set)~~{
	145	void resmethod_from_set ( const Setting &set ) {
138	146	string resmeth;
139		if (~~UI::get(resmeth,set,"resmethod",UI::optional))~~{
140		if (~~resmeth=="systematic"~~) {
141		resmethod= SYSTEMATIC;
	147	if ( UI::get ( resmeth, set, "resmethod", UI::optional ) ) {
	148	if ( resmeth == "systematic" ) {
	149	resmethod = SYSTEMATIC;
142	150	} else {
143		if (~~resmeth=="multinomial")~~{
144		resmethod=MULTINOMIAL;
	151	if ( resmeth == "multinomial" ) {
	152	resmethod = MULTINOMIAL;
145	153	} else {
146		if (~~resmeth=="stratified")~~{
147		resmethod= STRATIFIED;
	154	if ( resmeth == "stratified" ) {
	155	resmethod = STRATIFIED;
148	156	} else {
149		bdm_error~~("Unknown resampling method"~~);
	157	bdm_error ( "Unknown resampling method" );
150	158	}
151	159	}
152	160	}
153	161	} else {
154		resmethod=SYSTEMATIC;
	162	resmethod = SYSTEMATIC;
155	163	};
156		if~~(!UI::get(res_threshold, set, "res_threshold", UI::optional))~~{
157		res_threshold=0.5;
	164	if ( !UI::get ( res_threshold, set, "res_threshold", UI::optional ) ) {
	165	res_threshold = 0.5;
158	166	}
159	167	//validate();
160	168	}
161	169	//! load prior information from set and set internal structures accordingly
162		void prior_from_set~~(const Setting & set)~~{
163		shared_ptr<epdf> pri = UI::build<epdf>~~(set,"prior",UI::compulsory~~);
164
165		eEmp test_emp~~=dynamic_cast<eEmp>(&(*pri)~~);
166		if (~~test_emp~~) { // given pdf is sampled
167		est=*test_emp;
	170	void prior_from_set ( const Setting & set ) {
	171	shared_ptr<epdf> pri = UI::build<epdf> ( set, "prior", UI::compulsory );
	172
	173	eEmp test_emp = dynamic_cast<eEmp> ( & ( *pri ) );
	174	if ( test_emp ) { // given pdf is sampled
	175	est = *test_emp;
168	176	} else {
169	177	//int n;
170		if (!UI::get(n,set,"n",UI::optional)){n=10;}
	178	if ( !UI::get ( n, set, "n", UI::optional ) ) {
	179	n = 10;
	180	}
171	181	// sample from prior
172		set_statistics~~(ones(n)/n, *pri~~);
	182	set_statistics ( ones ( n ) / n, *pri );
173	183	}
174	184	n = est._w().length();
175	185	//validate();
176	186	}
177
178		void validate(){
179		bdm_assert~~(par,"PF::parameter_pdf not specified."~~);
180		n=_w.length();
181		lls~~=zeros(n~~);
182
183		if (~~par->_rv()._dsize()>0~~) {
184		bdm_assert~~(par->_rv()._dsize()==est.dimension(),~~"Mismatch of RV and dimension of posterior" );
	187
	188	void validate() {
	189	bdm_assert ( par, "PF::parameter_pdf not specified." );
	190	n = _w.length();
	191	lls = zeros ( n );
	192
	193	if ( par->_rv()._dsize() > 0 ) {
	194	bdm_assert ( par->_rv()._dsize() == est.dimension(), "Mismatch of RV and dimension of posterior" );
185	195	}
186	196	}
187	197	//! resample posterior density (from outside - see MPF)
188		void resample~~(ivec &ind)~~{
189		est.resample~~(ind,resmethod~~);
	198	void resample ( ivec &ind ) {
	199	est.resample ( ind, resmethod );
190	200	}
191	201	//! access function
192		Array<vec>& __samples(){return _samples;}
	202	Array<vec>& __samples() {
	203	return _samples;
	204	}
193	205	};
194		UIREGISTER~~(PF~~);
	206	UIREGISTER ( PF );
195	207
196	208	/*!
…	…
202	214
203	215	class MPF : public BM {
204		protected:
205		//! particle filter on non-linear variable
	216	protected:
	217	//! particle filter on non-linear variable
206	218	shared_ptr<PF> pf;
207	219	//! Array of Bayesian models
…	…
217	229	public:
218	230	//! constructor
219		mpfepdf (~~shared_ptr<PF> &pf0, Array<BM*> &BMs0): epdf(), pf(pf0), BMs(BMs0~~) { };
	231	mpfepdf ( shared_ptr<PF> &pf0, Array<BM*> &BMs0 ) : epdf(), pf ( pf0 ), BMs ( BMs0 ) { };
220	232	//! a variant of set parameters - this time, parameters are read from BMs and pf
221		void read_parameters(){
222		rv = concat~~(pf->posterior()._rv(), BMs(0)->posterior()._rv()~~);
223		dim = pf->posterior().dimension() + BMs(0)->posterior().dimension();
224		bdm_assert_debug~~(dim == rv._dsize(), "Wrong name "~~);
	233	void read_parameters() {
	234	rv = concat ( pf->posterior()._rv(), BMs ( 0 )->posterior()._rv() );
	235	dim = pf->posterior().dimension() + BMs ( 0 )->posterior().dimension();
	236	bdm_assert_debug ( dim == rv._dsize(), "Wrong name " );
225	237	}
226	238	vec mean() const {
227	239	const vec &w = pf->posterior()._w();
228		vec pom = zeros ( BMs(0)->posterior ().dimension() );
	240	vec pom = zeros ( BMs ( 0 )->posterior ().dimension() );
229	241	//compute mean of BMs
230	242	for ( int i = 0; i < w.length(); i++ ) {
…	…
235	247	vec variance() const {
236	248	const vec &w = pf->posterior()._w();
237
238		vec pom = zeros ( BMs(0)->posterior ().dimension() );
239		vec pom2 = zeros ( BMs(0)->posterior ().dimension() );
	249
	250	vec pom = zeros ( BMs ( 0 )->posterior ().dimension() );
	251	vec pom2 = zeros ( BMs ( 0 )->posterior ().dimension() );
240	252	vec mea;
241
	253
242	254	for ( int i = 0; i < w.length(); i++ ) {
243		// save current mean
	255	// save current mean
244	256	mea = BMs ( i )->posterior().mean();
245	257	pom += mea * w ( i );
…	…
249	261	return concat ( pf->posterior().variance(), pom2 - pow ( pom, 2 ) );
250	262	}
251
	263
252	264	void qbounds ( vec &lb, vec &ub, double perc = 0.95 ) const {
253	265	//bounds on particles
…	…
305	317	//!datalink from PF part to BM
306	318	datalink_part pf2bm;
307
	319
308	320	public:
309	321	//! Default constructor.
310		MPF () : jest (~~pf,BMs~~) {};
	322	MPF () : jest ( pf, BMs ) {};
311	323	//! set all parameters at once
312	324	void set_parameters ( shared_ptr<pdf> par0, shared_ptr<pdf> obs0, int n0, RESAMPLING_METHOD rm = SYSTEMATIC ) {
313		pf->set_model ( par0, obs0);
314		pf->set_parameters(n0, rm );
	325	pf->set_model ( par0, obs0 );
	326	pf->set_parameters ( n0, rm );
315	327	BMs.set_length ( n0 );
316	328	}
…	…
318	330	void set_BM ( const BM &BMcond0 ) {
319	331
320		int n=pf->__w().length();
321		BMs.set_length(n);
	332	int n = pf->__w().length();
	333	BMs.set_length ( n );
322	334	// copy
323	335	//BMcond0 .condition ( pf->posterior()._sample ( 0 ) );
…	…
331	343	return jest;
332	344	}
333		//! Extends options understood by BM::set_options by option
334		//! \li logmeans - meaning
	345	//! Extends options understood by BM::set_options by option
	346	//! \li logmeans - meaning
335	347	void set_options ( const string &opt ) {
336		BM::set_options~~(opt~~);
	348	BM::set_options ( opt );
337	349	}
338	350
…	…
341	353	return BMs ( i );
342	354	}
343
	355
344	356	/*! configuration structure for basic PF
345	357	\code
…	…
352	364	// resampling method
353	365	\endcode
354		*/
355		void from_setting~~(const Setting &set)~~{
356		shared_ptr<pdf> par = UI::build<pdf>~~(set,"parameter_pdf",UI::compulsory~~);
357		shared_ptr<pdf> obs= new pdf(); // not used!!
	366	*/
	367	void from_setting ( const Setting &set ) {
	368	shared_ptr<pdf> par = UI::build<pdf> ( set, "parameter_pdf", UI::compulsory );
	369	shared_ptr<pdf> obs = new pdf(); // not used!!
358	370
359	371	pf = new PF;
360	372	// prior must be set before BM
361		pf->prior_from_set~~(set~~);
362		pf->resmethod_from_set~~(set~~);
363		pf->set_model~~(par,obs~~);
364
365		shared_ptr<BM> BM0 =~~UI::build<BM>(set,"BM",UI::compulsory~~);
366		set_BM~~(*BM0~~);
367
	373	pf->prior_from_set ( set );
	374	pf->resmethod_from_set ( set );
	375	pf->set_model ( par, obs );
	376
	377	shared_ptr<BM> BM0 = UI::build<BM> ( set, "BM", UI::compulsory );
	378	set_BM ( *BM0 );
	379
368	380	string opt;
369		if (~~UI::get(opt,set,"options",UI::optional))~~{
370		set_options~~(opt~~);
	381	if ( UI::get ( opt, set, "options", UI::optional ) ) {
	382	set_options ( opt );
371	383	}
372	384	//set drv
373	385	//??set_yrv(concat(BM0->_yrv(),u) );
374		set_yrv(BM0->_yrv() );
375		rvc = BM0->_rvc().subt~~(par->_rv()~~);
	386	set_yrv ( BM0->_yrv() );
	387	rvc = BM0->_rvc().subt ( par->_rv() );
376	388	//find potential input - what remains in rvc when we subtract rv
377		RV u = par->_rvc().subt~~( par->_rv().copy_t(-1) );~~
378		rvc.add(u);
	389	RV u = par->_rvc().subt ( par->_rv().copy_t ( -1 ) );
	390	rvc.add ( u );
379	391	dimc = rvc._dsize();
380	392	validate();
381	393	}
382		void validate(){
383		try{
	394	void validate() {
	395	try {
384	396	pf->validate();
385		} catch (~~std::exception &e)~~{
386		throw UIException~~("Error in PF part of MPF:"~~);
	397	} catch ( std::exception &e ) {
	398	throw UIException ( "Error in PF part of MPF:" );
387	399	}
388	400	jest.read_parameters();
389		this2bm.set_connection~~(BMs(0)->_yrv(), BMs(0)->_rvc(), yrv, rvc~~);
390		this2pf.set_connection~~(pf->_yrv(), pf->_rvc(), yrv, rvc~~);
391		pf2bm.set_connection~~(BMs(0)->_rvc(), pf->posterior()._rv()~~);
	401	this2bm.set_connection ( BMs ( 0 )->_yrv(), BMs ( 0 )->_rvc(), yrv, rvc );
	402	this2pf.set_connection ( pf->_yrv(), pf->_rvc(), yrv, rvc );
	403	pf2bm.set_connection ( BMs ( 0 )->_rvc(), pf->posterior()._rv() );
392	404	}
393	405
394	406	};
395		UIREGISTER~~(MPF~~);
	407	UIREGISTER ( MPF );
396	408
397	409	}

Download in other formats: