Context Navigation

← Previous Change
Next Change →

exp_family.h

Timestamp:

06/09/10 14:00:40 (14 years ago)

Author:

mido

Message:

astyle applied all over the library

Files:

: 1 modified

library/bdm/stat/exp_family.h (modified) (39 diffs)

Legend:

: Unmodified
: Added
: Removed

library/bdm/stat/exp_family.h

r1063	r1064
37	37	public:
38	38	// eEF() :epdf() {};
39		//! default constructor
40		eEF () : epdf () {};
41		//! logarithm of the normalizing constant, \f$\mathcal{I}\f$
42		virtual double lognc() const = 0;
43
44		//!Evaluate normalized log-probability
45		virtual double evallog_nn ( const vec &val ) const NOT_IMPLEMENTED(0);
46
47		//!Evaluate normalized log-probability
48		virtual double evallog ( const vec &val ) const {
49		double tmp;
50		tmp = evallog_nn ( val ) - lognc();
51		return tmp;
52		}
53		//!Evaluate normalized log-probability for many samples
54		virtual vec evallog_mat ( const mat &Val ) const {
55		vec x ( Val.cols() );
56		for ( int i = 0; i < Val.cols(); i++ ) {
57		x ( i ) = evallog_nn ( Val.get_col ( i ) ) ;
58		}
59		return x - lognc();
60		}
61		//!Evaluate normalized log-probability for many samples
62		virtual vec evallog_mat ( const Array<vec> &Val ) const {
63		vec x ( Val.length() );
64		for ( int i = 0; i < Val.length(); i++ ) {
65		x ( i ) = evallog_nn ( Val ( i ) ) ;
66		}
67		return x - lognc();
68		}
69
70		//!Power of the density, used e.g. to flatten the density
71		virtual void pow ( double p ) NOT_IMPLEMENTED_VOID;
	39	//! default constructor
	40	eEF () : epdf () {};
	41	//! logarithm of the normalizing constant, \f$\mathcal{I}\f$
	42	virtual double lognc() const = 0;
	43
	44	//!Evaluate normalized log-probability
	45	virtual double evallog_nn ( const vec &val ) const NOT_IMPLEMENTED(0);
	46
	47	//!Evaluate normalized log-probability
	48	virtual double evallog ( const vec &val ) const {
	49	double tmp;
	50	tmp = evallog_nn ( val ) - lognc();
	51	return tmp;
	52	}
	53	//!Evaluate normalized log-probability for many samples
	54	virtual vec evallog_mat ( const mat &Val ) const {
	55	vec x ( Val.cols() );
	56	for ( int i = 0; i < Val.cols(); i++ ) {
	57	x ( i ) = evallog_nn ( Val.get_col ( i ) ) ;
	58	}
	59	return x - lognc();
	60	}
	61	//!Evaluate normalized log-probability for many samples
	62	virtual vec evallog_mat ( const Array<vec> &Val ) const {
	63	vec x ( Val.length() );
	64	for ( int i = 0; i < Val.length(); i++ ) {
	65	x ( i ) = evallog_nn ( Val ( i ) ) ;
	66	}
	67	return x - lognc();
	68	}
	69
	70	//!Power of the density, used e.g. to flatten the density
	71	virtual void pow ( double p ) NOT_IMPLEMENTED_VOID;
72	72	};
73	73
…	…
75	75	//! Estimator for Exponential family
76	76	class BMEF : public BM {
77		public:
78		//! forgetting factor
79		double frg;
80		protected:
81		//! cached value of lognc() in the previous step (used in evaluation of \c ll )
82		double last_lognc;
83		//! factor k = [0..1] for scheduling of forgetting factor: \f$ frg_t = (1-k) * frg_{t-1} + k \f$, default 0
84		double frg_sched_factor;
85		public:
86		//! Default constructor (=empty constructor)
87		BMEF ( double frg0 = 1.0 ) : BM (), frg ( frg0 ), last_lognc(0.0),frg_sched_factor(0.0) {}
88		//! Copy constructor
89		BMEF ( const BMEF &B ) : BM ( B ), frg ( B.frg ), last_lognc ( B.last_lognc ),frg_sched_factor(B.frg_sched_factor) {}
90		//!get statistics from another model
91		virtual void set_statistics ( const BMEF* BM0 ) NOT_IMPLEMENTED_VOID;
92
93		//! Weighted update of sufficient statistics (Bayes rule)
94		virtual void bayes_weighted ( const vec &data, const vec &cond = empty_vec, const double w = 1.0 ) {
95		if (frg_sched_factor>0) {frg = frg*(1-frg_sched_factor)+frg_sched_factor;}
96		};
97		//original Bayes
98		void bayes ( const vec &yt, const vec &cond = empty_vec );
99
100		//!Flatten the posterior according to the given BMEF (of the same type!)
101		virtual void flatten ( const BMEF * B, double weight=1.0 ) NOT_IMPLEMENTED_VOID;;
102
103
104		void to_setting ( Setting &set ) const
105		{
106		BM::to_setting( set );
107		UI::save(frg, set, "frg");
108		UI::save( frg_sched_factor, set, "frg_sched_factor" );
109		}
110
111		void from_setting( const Setting &set) {
112		BM::from_setting(set);
113		if ( !UI::get ( frg, set, "frg" ) )
114		frg = 1.0;
115		UI::get ( frg_sched_factor, set, "frg_sched_factor",UI::optional );
116		}
117
118		void validate() {
119		BM::validate();
120		}
	77	public:
	78	//! forgetting factor
	79	double frg;
	80	protected:
	81	//! cached value of lognc() in the previous step (used in evaluation of \c ll )
	82	double last_lognc;
	83	//! factor k = [0..1] for scheduling of forgetting factor: \f$ frg_t = (1-k) * frg_{t-1} + k \f$, default 0
	84	double frg_sched_factor;
	85	public:
	86	//! Default constructor (=empty constructor)
	87	BMEF ( double frg0 = 1.0 ) : BM (), frg ( frg0 ), last_lognc(0.0),frg_sched_factor(0.0) {}
	88	//! Copy constructor
	89	BMEF ( const BMEF &B ) : BM ( B ), frg ( B.frg ), last_lognc ( B.last_lognc ),frg_sched_factor(B.frg_sched_factor) {}
	90	//!get statistics from another model
	91	virtual void set_statistics ( const BMEF* BM0 ) NOT_IMPLEMENTED_VOID;
	92
	93	//! Weighted update of sufficient statistics (Bayes rule)
	94	virtual void bayes_weighted ( const vec &data, const vec &cond = empty_vec, const double w = 1.0 ) {
	95	if (frg_sched_factor>0) {
	96	frg = frg*(1-frg_sched_factor)+frg_sched_factor;
	97	}
	98	};
	99	//original Bayes
	100	void bayes ( const vec &yt, const vec &cond = empty_vec );
	101
	102	//!Flatten the posterior according to the given BMEF (of the same type!)
	103	virtual void flatten ( const BMEF * B, double weight=1.0 ) NOT_IMPLEMENTED_VOID;;
	104
	105
	106	void to_setting ( Setting &set ) const
	107	{
	108	BM::to_setting( set );
	109	UI::save(frg, set, "frg");
	110	UI::save( frg_sched_factor, set, "frg_sched_factor" );
	111	}
	112
	113	void from_setting( const Setting &set) {
	114	BM::from_setting(set);
	115	if ( !UI::get ( frg, set, "frg" ) )
	116	frg = 1.0;
	117	UI::get ( frg_sched_factor, set, "frg_sched_factor",UI::optional );
	118	}
	119
	120	void validate() {
	121	BM::validate();
	122	}
121	123
122	124	};
…	…
127	129	where \f$ x_t \f$ is the \c rv, \f$ y_t \f$ is the \c rvc and g is a deterministic transformation of class fn.
128	130	*/
129		class mgdirac: public pdf{
130		protected:
131		shared_ptr<fnc> g;
132		public:
133		vec samplecond(const vec &cond) {
134		bdm_assert_debug(cond.length()==g->dimensionc(),"given cond in not compatible with g");
135		vec tmp = g->eval(cond);
136		return tmp;
137		}
138		~~double evallogcond ( const vec &yt, const vec &cond )~~{
139		return std::numeric_limits< double >::max();
140		}
141		void from_setting(const Setting& set);
142		void to_setting(Setting &set) const;
143		void validate();
	131	class mgdirac: public pdf {
	132	protected:
	133	shared_ptr<fnc> g;
	134	public:
	135	vec samplecond(const vec &cond) {
	136	bdm_assert_debug(cond.length()==g->dimensionc(),"given cond in not compatible with g");
	137	vec tmp = g->eval(cond);
	138	return tmp;
	139	}
	140	double evallogcond ( const vec &yt, const vec &cond ) {
	141	return std::numeric_limits< double >::max();
	142	}
	143	void from_setting(const Setting& set);
	144	void to_setting(Setting &set) const;
	145	void validate();
144	146	};
145	147	UIREGISTER(mgdirac);
…	…
157	159	class enorm : public eEF {
158	160	protected:
159		//! mean value
160		vec mu;
161		//! Covariance matrix in decomposed form
162		sq_T R;
163		public:
164		//!\name Constructors
165		//!@{
166
167		enorm () : eEF (), mu (), R () {};
168		enorm ( const vec &mu, const sq_T &R ) {
169		set_parameters ( mu, R );
170		}
171		void set_parameters ( const vec &mu, const sq_T &R );
172		/*! Create Normal density
173		\f[ f(rv) = N(\mu, R) \f]
174		from structure
175		\code
176		class = 'enorm<ldmat>', (OR) 'enorm<chmat>', (OR) 'enorm<fsqmat>';
177		mu = []; // mean value
178		R = []; // variance, square matrix of appropriate dimension
179		\endcode
180		*/
181		void from_setting ( const Setting &root );
182		void to_setting ( Setting &root ) const ;
183
184		void validate();
185		//!@}
186
187		//! \name Mathematical operations
188		//!@{
189
190		//! dupdate in exponential form (not really handy)
191		void dupdate ( mat &v, double nu = 1.0 );
192
193		//! evaluate bhattacharya distance
194		~~double bhattacharyya(const enorm<sq_T> &e2)~~{
195		bdm_assert(dim == e2.dimension(), "enorms of differnt dimensions");
196		sq_T P=R;
197		P.add(e2._R());
198
199		double tmp = 0.125P.invqform(mu - e2._mu()) + 0.5(P.logdet() - 0.5*(R.logdet() + e2._R().logdet()));
200		return tmp;
201		}
202
203		vec sample() const;
204
205		double evallog_nn ( const vec &val ) const;
206		double lognc () const;
207		vec mean() const {
208		return mu;
209		}
210		vec variance() const {
211		return diag ( R.to_mat() );
212		}
213		mat covariance() const {
214		return R.to_mat();
215		}
216		// mlnorm<sq_T>* condition ( const RV &rvn ) const ; <=========== fails to cmpile. Why?
217		shared_ptr<pdf> condition ( const RV &rvn ) const;
218
219		// target not typed to mlnorm<sq_T, enorm<sq_T> > &
220		// because that doesn't compile (perhaps because we
221		// haven't finished defining enorm yet), but the type
222		// is required
223		void condition ( const RV &rvn, pdf &target ) const;
224
225		shared_ptr<epdf> marginal ( const RV &rvn ) const;
226		void marginal ( const RV &rvn, enorm<sq_T> &target ) const;
227		//!@}
228
229		//! \name Access to attributes
230		//!@{
231
232		vec& _mu() {
233		return mu;
234		}
235		const vec& _mu() const {
236		return mu;
237		}
238		void set_mu ( const vec mu0 ) {
239		mu = mu0;
240		}
241		sq_T& _R() {
242		return R;
243		}
244		const sq_T& _R() const {
245		return R;
246		}
247		//!@}
	161	//! mean value
	162	vec mu;
	163	//! Covariance matrix in decomposed form
	164	sq_T R;
	165	public:
	166	//!\name Constructors
	167	//!@{
	168
	169	enorm () : eEF (), mu (), R () {};
	170	enorm ( const vec &mu, const sq_T &R ) {
	171	set_parameters ( mu, R );
	172	}
	173	void set_parameters ( const vec &mu, const sq_T &R );
	174	/*! Create Normal density
	175	\f[ f(rv) = N(\mu, R) \f]
	176	from structure
	177	\code
	178	class = 'enorm<ldmat>', (OR) 'enorm<chmat>', (OR) 'enorm<fsqmat>';
	179	mu = []; // mean value
	180	R = []; // variance, square matrix of appropriate dimension
	181	\endcode
	182	*/
	183	void from_setting ( const Setting &root );
	184	void to_setting ( Setting &root ) const ;
	185
	186	void validate();
	187	//!@}
	188
	189	//! \name Mathematical operations
	190	//!@{
	191
	192	//! dupdate in exponential form (not really handy)
	193	void dupdate ( mat &v, double nu = 1.0 );
	194
	195	//! evaluate bhattacharya distance
	196	double bhattacharyya(const enorm<sq_T> &e2) {
	197	bdm_assert(dim == e2.dimension(), "enorms of differnt dimensions");
	198	sq_T P=R;
	199	P.add(e2._R());
	200
	201	double tmp = 0.125P.invqform(mu - e2._mu()) + 0.5(P.logdet() - 0.5*(R.logdet() + e2._R().logdet()));
	202	return tmp;
	203	}
	204
	205	vec sample() const;
	206
	207	double evallog_nn ( const vec &val ) const;
	208	double lognc () const;
	209	vec mean() const {
	210	return mu;
	211	}
	212	vec variance() const {
	213	return diag ( R.to_mat() );
	214	}
	215	mat covariance() const {
	216	return R.to_mat();
	217	}
	218	// mlnorm<sq_T>* condition ( const RV &rvn ) const ; <=========== fails to cmpile. Why?
	219	shared_ptr<pdf> condition ( const RV &rvn ) const;
	220
	221	// target not typed to mlnorm<sq_T, enorm<sq_T> > &
	222	// because that doesn't compile (perhaps because we
	223	// haven't finished defining enorm yet), but the type
	224	// is required
	225	void condition ( const RV &rvn, pdf &target ) const;
	226
	227	shared_ptr<epdf> marginal ( const RV &rvn ) const;
	228	void marginal ( const RV &rvn, enorm<sq_T> &target ) const;
	229	//!@}
	230
	231	//! \name Access to attributes
	232	//!@{
	233
	234	vec& _mu() {
	235	return mu;
	236	}
	237	const vec& _mu() const {
	238	return mu;
	239	}
	240	void set_mu ( const vec mu0 ) {
	241	mu = mu0;
	242	}
	243	sq_T& _R() {
	244	return R;
	245	}
	246	const sq_T& _R() const {
	247	return R;
	248	}
	249	//!@}
248	250
249	251	};
…	…
255	257	SHAREDPTR2 ( enorm, fsqmat );
256	258
257		//! \class bdm::egauss
	259	//! \class bdm::egauss
258	260	//!\brief Gaussian (Normal) distribution. Same as enorm<fsqmat>.
259	261	typedef enorm<ldmat> egauss;
…	…
266	268	/*! distribution of multivariate Student t density
267	269
268		Based on article by Genest and Zidek,
	270	Based on article by Genest and Zidek,
269	271	*/
270	272	template<class sq_T>
271		class estudent : public eEF{
272		protected:
273		//! mena value
274		vec mu;
275		//! matrix H
276		sq_T H;
277		//! degrees of freedom
278		double delta;
279		public:
280		double evallog_nn(const vec &val) const{
281		double tmp = -0.5H.logdet() - 0.5(delta + dim) * log(1+ H.invqform(val - mu)/delta);
282		return tmp;
283		}
284		double lognc() const {
285		//log(pi) = 1.14472988584940
286		double tmp = -lgamma(0.5(delta+dim))+lgamma(0.5delta) + 0.5dim(log(delta) + 1.14472988584940);
287		return tmp;
288		}
289		void marginal (const RV &rvm, estudent<sq_T> &marg) const {
290		ivec ind = rvm.findself_ids(rv); // indices of rvm in rv
291		marg._mu() = mu(ind);
292		marg._H() = sq_T(H,ind);
293		marg._delta() = delta;
294		marg.validate();
295		}
296		shared_ptr<epdf> marginal(const RV &rvm) const {
297		shared_ptr<estudent<sq_T> > tmp = new estudent<sq_T>;
298		marginal(rvm, *tmp);
299		return tmp;
300		}
301		vec sample() const NOT_IMPLEMENTED(vec(0))
302
303		vec mean() const {return mu;}
304		mat covariance() const {
305		return delta/(delta-2)*H.to_mat();
306		}
307		vec variance() const {return diag(covariance());}
308		//! \name access
309		//! @{
310		//! access function
311		vec& _mu() {return mu;}
312		//! access function
313		sq_T& _H() {return H;}
314		//! access function
315		double& _delta() {return delta;}
316		//!@}
317		//! todo
318		void from_setting(const Setting &set){
319		epdf::from_setting(set);
320		mat H0;
321		UI::get(H0,set, "H");
322		H= H0; // conversion!!
323		UI::get(delta,set,"delta");
324		UI::get(mu,set,"mu");
325		}
326		void to_setting(Setting &set) const{
327		epdf::to_setting(set);
328		UI::save(H.to_mat(), set, "H");
329		UI::save(delta, set, "delta");
330		UI::save(mu, set, "mu");
331		}
332		void validate() {
333		eEF::validate();
334		dim = H.rows();
335		}
	273	class estudent : public eEF {
	274	protected:
	275	//! mena value
	276	vec mu;
	277	//! matrix H
	278	sq_T H;
	279	//! degrees of freedom
	280	double delta;
	281	public:
	282	double evallog_nn(const vec &val) const {
	283	double tmp = -0.5H.logdet() - 0.5(delta + dim) * log(1+ H.invqform(val - mu)/delta);
	284	return tmp;
	285	}
	286	double lognc() const {
	287	//log(pi) = 1.14472988584940
	288	double tmp = -lgamma(0.5(delta+dim))+lgamma(0.5delta) + 0.5dim(log(delta) + 1.14472988584940);
	289	return tmp;
	290	}
	291	void marginal (const RV &rvm, estudent<sq_T> &marg) const {
	292	ivec ind = rvm.findself_ids(rv); // indices of rvm in rv
	293	marg._mu() = mu(ind);
	294	marg._H() = sq_T(H,ind);
	295	marg._delta() = delta;
	296	marg.validate();
	297	}
	298	shared_ptr<epdf> marginal(const RV &rvm) const {
	299	shared_ptr<estudent<sq_T> > tmp = new estudent<sq_T>;
	300	marginal(rvm, *tmp);
	301	return tmp;
	302	}
	303	vec sample() const NOT_IMPLEMENTED(vec(0))
	304
	305	vec mean() const {
	306	return mu;
	307	}
	308	mat covariance() const {
	309	return delta/(delta-2)*H.to_mat();
	310	}
	311	vec variance() const {
	312	return diag(covariance());
	313	}
	314	//! \name access
	315	//! @{
	316	//! access function
	317	vec& _mu() {
	318	return mu;
	319	}
	320	//! access function
	321	sq_T& _H() {
	322	return H;
	323	}
	324	//! access function
	325	double& _delta() {
	326	return delta;
	327	}
	328	//!@}
	329	//! todo
	330	void from_setting(const Setting &set) {
	331	epdf::from_setting(set);
	332	mat H0;
	333	UI::get(H0,set, "H");
	334	H= H0; // conversion!!
	335	UI::get(delta,set,"delta");
	336	UI::get(mu,set,"mu");
	337	}
	338	void to_setting(Setting &set) const {
	339	epdf::to_setting(set);
	340	UI::save(H.to_mat(), set, "H");
	341	UI::save(delta, set, "delta");
	342	UI::save(mu, set, "mu");
	343	}
	344	void validate() {
	345	eEF::validate();
	346	dim = H.rows();
	347	}
336	348	};
337	349	UIREGISTER2(estudent,fsqmat);
…	…
346	358	*/
347	359	class egiw : public eEF {
348		//! \var log_level_enums logvartheta
349		//! Log variance of the theta part
350
351		LOG_LEVEL(egiw,logvartheta);
352
353		protected:
354		//! Extended information matrix of sufficient statistics
355		ldmat V;
356		//! Number of data records (degrees of freedom) of sufficient statistics
357		double nu;
358		//! Dimension of the output
359		int dimx;
360		//! Dimension of the regressor
361		int nPsi;
362		public:
363		//!\name Constructors
364		//!@{
365		egiw() : eEF(),dimx(0) {};
366		egiw ( int dimx0, ldmat V0, double nu0 = -1.0 ) : eEF(),dimx(0) {
367		set_parameters ( dimx0, V0, nu0 );
368		validate();
369		};
370
371		void set_parameters ( int dimx0, ldmat V0, double nu0 = -1.0 );
372		//!@}
373
374		vec sample() const;
375		mat sample_mat ( int n ) const;
376		vec mean() const;
377		vec variance() const;
378		//mat covariance() const;
379		void sample_mat ( mat &Mi, chmat &Ri ) const;
380
381		void factorize ( mat &M, ldmat &Vz, ldmat &Lam ) const;
382		//! LS estimate of \f$\theta\f$
383		vec est_theta() const;
384
385		//! Covariance of the LS estimate
386		ldmat est_theta_cov() const;
387
388		//! expected values of the linear coefficient and the covariance matrix are written to \c M and \c R , respectively
389		void mean_mat ( mat &M, mat&R ) const;
390		//! In this instance, val= [theta, r]. For multivariate instances, it is stored columnwise val = [theta_1 theta_2 ... r_1 r_2 ]
391		double evallog_nn ( const vec &val ) const;
392		double lognc () const;
393		void pow ( double p ) {
394		V *= p;
395		nu *= p;
396		};
397
398		//! marginal density (only student for now)
399		shared_ptr<epdf> marginal(const RV &rvm) const {
400		bdm_assert(dimx==1, "Not supported");
401		//TODO - this is too trivial!!!
402		ivec ind = rvm.findself_ids(rv);
403		if (min(ind)==0) { //assume it si
404		shared_ptr<estudent<ldmat> > tmp = new estudent<ldmat>;
405		mat M;
406		ldmat Vz;
407		ldmat Lam;
408		factorize(M,Vz,Lam);
409
410		tmp->_mu() = M.get_col(0);
411		ldmat H;
412		Vz.inv(H);
413		H *=Lam._D()(0)/nu;
414		tmp->_H() = H;
415		tmp->_delta() = nu;
416		tmp->validate();
417		return tmp;
418		}
419		return NULL;
420		}
421		//! \name Access attributes
422		//!@{
423
424		ldmat& _V() {
425		return V;
426		}
427		const ldmat& _V() const {
428		return V;
429		}
430		double& _nu() {
431		return nu;
432		}
433		const double& _nu() const {
434		return nu;
435		}
436		const int & _dimx() const {
437		return dimx;
438		}
439
440		/*! Create Gauss-inverse-Wishart density
441		\f[ f(rv) = GiW(V,\nu) \f]
442		from structure
443		\code
444		class = 'egiw';
445		V.L = []; // L part of matrix V
446		V.D = []; // D part of matrix V
447		-or- V = [] // full matrix V
448		-or- dV = []; // vector of diagonal of V (when V not given)
449		nu = []; // scalar \nu ((almost) degrees of freedom)
450		// when missing, it will be computed to obtain proper pdf
451		dimx = []; // dimension of the wishart part
452		rv = RV({'name'}) // description of RV
453		rvc = RV({'name'}) // description of RV in condition
454		\endcode
455
456		\sa log_level_enums
457		*/
458		void from_setting ( const Setting &set );
459		//! see egiw::from_setting
460		void to_setting ( Setting& set ) const;
461		void validate();
462		void log_register ( bdm::logger& L, const string& prefix );
463
464		void log_write() const;
465		//!@}
	360	//! \var log_level_enums logvartheta
	361	//! Log variance of the theta part
	362
	363	LOG_LEVEL(egiw,logvartheta);
	364
	365	protected:
	366	//! Extended information matrix of sufficient statistics
	367	ldmat V;
	368	//! Number of data records (degrees of freedom) of sufficient statistics
	369	double nu;
	370	//! Dimension of the output
	371	int dimx;
	372	//! Dimension of the regressor
	373	int nPsi;
	374	public:
	375	//!\name Constructors
	376	//!@{
	377	egiw() : eEF(),dimx(0) {};
	378	egiw ( int dimx0, ldmat V0, double nu0 = -1.0 ) : eEF(),dimx(0) {
	379	set_parameters ( dimx0, V0, nu0 );
	380	validate();
	381	};
	382
	383	void set_parameters ( int dimx0, ldmat V0, double nu0 = -1.0 );
	384	//!@}
	385
	386	vec sample() const;
	387	mat sample_mat ( int n ) const;
	388	vec mean() const;
	389	vec variance() const;
	390	//mat covariance() const;
	391	void sample_mat ( mat &Mi, chmat &Ri ) const;
	392
	393	void factorize ( mat &M, ldmat &Vz, ldmat &Lam ) const;
	394	//! LS estimate of \f$\theta\f$
	395	vec est_theta() const;
	396
	397	//! Covariance of the LS estimate
	398	ldmat est_theta_cov() const;
	399
	400	//! expected values of the linear coefficient and the covariance matrix are written to \c M and \c R , respectively
	401	void mean_mat ( mat &M, mat&R ) const;
	402	//! In this instance, val= [theta, r]. For multivariate instances, it is stored columnwise val = [theta_1 theta_2 ... r_1 r_2 ]
	403	double evallog_nn ( const vec &val ) const;
	404	double lognc () const;
	405	void pow ( double p ) {
	406	V *= p;
	407	nu *= p;
	408	};
	409
	410	//! marginal density (only student for now)
	411	shared_ptr<epdf> marginal(const RV &rvm) const {
	412	bdm_assert(dimx==1, "Not supported");
	413	//TODO - this is too trivial!!!
	414	ivec ind = rvm.findself_ids(rv);
	415	if (min(ind)==0) { //assume it si
	416	shared_ptr<estudent<ldmat> > tmp = new estudent<ldmat>;
	417	mat M;
	418	ldmat Vz;
	419	ldmat Lam;
	420	factorize(M,Vz,Lam);
	421
	422	tmp->_mu() = M.get_col(0);
	423	ldmat H;
	424	Vz.inv(H);
	425	H *=Lam._D()(0)/nu;
	426	tmp->_H() = H;
	427	tmp->_delta() = nu;
	428	tmp->validate();
	429	return tmp;
	430	}
	431	return NULL;
	432	}
	433	//! \name Access attributes
	434	//!@{
	435
	436	ldmat& _V() {
	437	return V;
	438	}
	439	const ldmat& _V() const {
	440	return V;
	441	}
	442	double& _nu() {
	443	return nu;
	444	}
	445	const double& _nu() const {
	446	return nu;
	447	}
	448	const int & _dimx() const {
	449	return dimx;
	450	}
	451
	452	/*! Create Gauss-inverse-Wishart density
	453	\f[ f(rv) = GiW(V,\nu) \f]
	454	from structure
	455	\code
	456	class = 'egiw';
	457	V.L = []; // L part of matrix V
	458	V.D = []; // D part of matrix V
	459	-or- V = [] // full matrix V
	460	-or- dV = []; // vector of diagonal of V (when V not given)
	461	nu = []; // scalar \nu ((almost) degrees of freedom)
	462	// when missing, it will be computed to obtain proper pdf
	463	dimx = []; // dimension of the wishart part
	464	rv = RV({'name'}) // description of RV
	465	rvc = RV({'name'}) // description of RV in condition
	466	\endcode
	467
	468	\sa log_level_enums
	469	*/
	470	void from_setting ( const Setting &set );
	471	//! see egiw::from_setting
	472	void to_setting ( Setting& set ) const;
	473	void validate();
	474	void log_register ( bdm::logger& L, const string& prefix );
	475
	476	void log_write() const;
	477	//!@}
466	478	};
467	479	UIREGISTER ( egiw );
…	…
478	490	class eDirich: public eEF {
479	491	protected:
480		//!sufficient statistics
481		vec beta;
482		public:
483		//!\name Constructors
484		//!@{
485
486		eDirich () : eEF () {};
487		eDirich ( const eDirich &D0 ) : eEF () {
488		set_parameters ( D0.beta );
489		validate();
490		};
491		eDirich ( const vec &beta0 ) {
492		set_parameters ( beta0 );
493		validate();
494		};
495		void set_parameters ( const vec &beta0 ) {
496		beta = beta0;
497		dim = beta.length();
498		}
499		//!@}
500
501		//! using sampling procedure from wikipedia
502		vec sample() const {
503		vec y ( beta.length() );
504		for ( int i = 0; i < beta.length(); i++ ) {
505		GamRNG.setup ( beta ( i ), 1 );
	492	//!sufficient statistics
	493	vec beta;
	494	public:
	495	//!\name Constructors
	496	//!@{
	497
	498	eDirich () : eEF () {};
	499	eDirich ( const eDirich &D0 ) : eEF () {
	500	set_parameters ( D0.beta );
	501	validate();
	502	};
	503	eDirich ( const vec &beta0 ) {
	504	set_parameters ( beta0 );
	505	validate();
	506	};
	507	void set_parameters ( const vec &beta0 ) {
	508	beta = beta0;
	509	dim = beta.length();
	510	}
	511	//!@}
	512
	513	//! using sampling procedure from wikipedia
	514	vec sample() const {
	515	vec y ( beta.length() );
	516	for ( int i = 0; i < beta.length(); i++ ) {
	517	GamRNG.setup ( beta ( i ), 1 );
506	518	#pragma omp critical
507		y ( i ) = GamRNG();
508		}
509		return y / sum ( y );
510		}
511
512		vec mean() const {
513		return beta / sum ( beta );
514		};
515		vec variance() const {
516		double gamma = sum ( beta );
517		return elem_mult ( beta, ( gamma - beta ) ) / ( gammagamma ( gamma + 1 ) );
518		}
519		//! In this instance, val is ...
520		double evallog_nn ( const vec &val ) const {
521		double tmp;
522		tmp = ( beta - 1 ) * log ( val );
523		return tmp;
524		}
525
526		double lognc () const {
527		double tmp;
528		double gam = sum ( beta );
529		double lgb = 0.0;
530		for ( int i = 0; i < beta.length(); i++ ) {
531		lgb += lgamma ( beta ( i ) );
532		}
533		tmp = lgb - lgamma ( gam );
534		return tmp;
535		}
536
537		//!access function
538		vec& _beta() {
539		return beta;
540		}
541
542		/*! Create object from the following structure
543		\code
544		class = 'eDirich';
545		beta = [...]; % vector parameter beta
546		--- inherited fields ---
547		bdm::eEF::from_setting
548		\endcode
549		*/
550		void from_setting ( const Setting &set );
551
552		void validate();
553
554		void to_setting ( Setting &set ) const;
	519	y ( i ) = GamRNG();
	520	}
	521	return y / sum ( y );
	522	}
	523
	524	vec mean() const {
	525	return beta / sum ( beta );
	526	};
	527	vec variance() const {
	528	double gamma = sum ( beta );
	529	return elem_mult ( beta, ( gamma - beta ) ) / ( gammagamma ( gamma + 1 ) );
	530	}
	531	//! In this instance, val is ...
	532	double evallog_nn ( const vec &val ) const {
	533	double tmp;
	534	tmp = ( beta - 1 ) * log ( val );
	535	return tmp;
	536	}
	537
	538	double lognc () const {
	539	double tmp;
	540	double gam = sum ( beta );
	541	double lgb = 0.0;
	542	for ( int i = 0; i < beta.length(); i++ ) {
	543	lgb += lgamma ( beta ( i ) );
	544	}
	545	tmp = lgb - lgamma ( gam );
	546	return tmp;
	547	}
	548
	549	//!access function
	550	vec& _beta() {
	551	return beta;
	552	}
	553
	554	/*! Create object from the following structure
	555	\code
	556	class = 'eDirich';
	557	beta = [...]; % vector parameter beta
	558	--- inherited fields ---
	559	bdm::eEF::from_setting
	560	\endcode
	561	*/
	562	void from_setting ( const Setting &set );
	563
	564	void validate();
	565
	566	void to_setting ( Setting &set ) const;
555	567	};
556	568	UIREGISTER ( eDirich );
…	…
560	572	Continuous Dirichlet density of \f$n\f$-dimensional variable \f$x\f$
561	573	\f[
562		f(x\|\alpha,\beta) = \frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)}\, x^{\alpha-1}(1-x)^{\beta-1}
	574	f(x\|\alpha,\beta) = \frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)}\, x^{\alpha-1}(1-x)^{\beta-1}
563	575	\f]
564	576	is a simplification of Dirichlet to univariate case.
565	577	*/
566	578	class eBeta: public eEF {
567		public:
568		//!sufficient statistics
569		vec alpha;
570		//!sufficient statistics
571		vec beta;
572		public:
573		//!\name Constructors
574		//!@{
575
576		eBeta () : eEF () {};
577		eBeta ( const eBeta &B0 ) : eEF (), alpha(B0.alpha),beta(B0.beta) {
578		validate();
579		};
580		//!@}
581
582		//! using sampling procedure from wikipedia
583		vec sample() const {
584		vec y ( beta.length() ); // all vectors
585		for ( int i = 0; i < beta.length(); i++ ) {
586		GamRNG.setup ( alpha ( i ), 1 );
587		#pragma omp critical
588		double Ga = GamRNG();
589
590		GamRNG.setup ( beta ( i ), 1 );
591		#pragma omp critical
592		double Gb = GamRNG();
593
594		y ( i ) = Ga/(Ga+Gb);
595		}
596		return y;
597		}
598
599		vec mean() const {
600		return elem_div(alpha, alpha + beta); // dot-division
601		};
602		vec variance() const {
603		vec apb=alpha+beta;
604		return elem_div (elem_mult ( alpha, beta) ,
605		elem_mult ( elem_mult(apb,apb), apb+1 ) );
606		}
607		//! In this instance, val is ...
608		double evallog_nn ( const vec &val ) const {
609		double tmp;
610		tmp = ( alpha - 1 ) * log ( val ) + (beta-1)*log(1-val);
611		return tmp;
612		}
613
614		double lognc () const {
615		double lgb = 0.0;
616		for ( int i = 0; i < beta.length(); i++ ) {
617		lgb += -lgamma ( alpha(i)+beta(i) ) + lgamma(alpha(i)) + lgamma(beta(i));
618		}
619		return lgb;
620		}
621
622		/*! Create object from the following structure
623
624		\code
625		class = 'eBeta';
626		alpha = [...]; % vector parameter alpha
627		beta = [...]; % vector parameter beta of the same length as alpha
628		\endcode
629
630		Class does not call bdm::eEF::from_setting
631		*/
632		~~void from_setting ( const Setting &set )~~{
633		UI::get(alpha, set, "alpha", UI::compulsory);
634		UI::get(beta, set, "beta", UI::compulsory);
635		}
636
637		~~void validate()~~{
638		bdm_assert(alpha.length()==beta.length(), "eBeta:: alpha and beta length do not match");
639		dim = alpha.length();
640		}
641
642		~~void to_setting ( Setting &set ) const~~{
643		UI::save(alpha, set, "alpha");
644		UI::save(beta, set, "beta");
645		}
	579	public:
	580	//!sufficient statistics
	581	vec alpha;
	582	//!sufficient statistics
	583	vec beta;
	584	public:
	585	//!\name Constructors
	586	//!@{
	587
	588	eBeta () : eEF () {};
	589	eBeta ( const eBeta &B0 ) : eEF (), alpha(B0.alpha),beta(B0.beta) {
	590	validate();
	591	};
	592	//!@}
	593
	594	//! using sampling procedure from wikipedia
	595	vec sample() const {
	596	vec y ( beta.length() ); // all vectors
	597	for ( int i = 0; i < beta.length(); i++ ) {
	598	GamRNG.setup ( alpha ( i ), 1 );
	599	#pragma omp critical
	600	double Ga = GamRNG();
	601
	602	GamRNG.setup ( beta ( i ), 1 );
	603	#pragma omp critical
	604	double Gb = GamRNG();
	605
	606	y ( i ) = Ga/(Ga+Gb);
	607	}
	608	return y;
	609	}
	610
	611	vec mean() const {
	612	return elem_div(alpha, alpha + beta); // dot-division
	613	};
	614	vec variance() const {
	615	vec apb=alpha+beta;
	616	return elem_div (elem_mult ( alpha, beta) ,
	617	elem_mult ( elem_mult(apb,apb), apb+1 ) );
	618	}
	619	//! In this instance, val is ...
	620	double evallog_nn ( const vec &val ) const {
	621	double tmp;
	622	tmp = ( alpha - 1 ) * log ( val ) + (beta-1)*log(1-val);
	623	return tmp;
	624	}
	625
	626	double lognc () const {
	627	double lgb = 0.0;
	628	for ( int i = 0; i < beta.length(); i++ ) {
	629	lgb += -lgamma ( alpha(i)+beta(i) ) + lgamma(alpha(i)) + lgamma(beta(i));
	630	}
	631	return lgb;
	632	}
	633
	634	/*! Create object from the following structure
	635
	636	\code
	637	class = 'eBeta';
	638	alpha = [...]; % vector parameter alpha
	639	beta = [...]; % vector parameter beta of the same length as alpha
	640	\endcode
	641
	642	Class does not call bdm::eEF::from_setting
	643	*/
	644	void from_setting ( const Setting &set ) {
	645	UI::get(alpha, set, "alpha", UI::compulsory);
	646	UI::get(beta, set, "beta", UI::compulsory);
	647	}
	648
	649	void validate() {
	650	bdm_assert(alpha.length()==beta.length(), "eBeta:: alpha and beta length do not match");
	651	dim = alpha.length();
	652	}
	653
	654	void to_setting ( Setting &set ) const {
	655	UI::save(alpha, set, "alpha");
	656	UI::save(beta, set, "beta");
	657	}
646	658	};
647	659	UIREGISTER ( eBeta );
…	…
660	672	class mDirich: public pdf_internal<eDirich> {
661	673	protected:
662		//! constant \f$ k \f$ of the random walk
663		double k;
664		//! cache of beta_i
665		vec &_beta;
666		//! stabilizing coefficient \f$ \beta_c \f$
667		vec betac;
668		public:
669		mDirich() : pdf_internal<eDirich>(), _beta ( iepdf._beta() ) {};
670		void condition ( const vec &val ) {
671		_beta = val / k + betac;
672		};
673		/*! Create Dirichlet random walk
674		\f[ f(rv\|rvc) = Di(rvc*k) \f]
675		from structure
676		\code
677		class = 'mDirich';
678		k = 1; // multiplicative constant k
679		--- optional ---
680		rv = RV({'name'},size) // description of RV
681		beta0 = []; // initial value of beta
682		betac = []; // initial value of beta
683		\endcode
684		*/
685		void from_setting ( const Setting &set );
686		void to_setting (Setting &set) const;
687		void validate();
	674	//! constant \f$ k \f$ of the random walk
	675	double k;
	676	//! cache of beta_i
	677	vec &_beta;
	678	//! stabilizing coefficient \f$ \beta_c \f$
	679	vec betac;
	680	public:
	681	mDirich() : pdf_internal<eDirich>(), _beta ( iepdf._beta() ) {};
	682	void condition ( const vec &val ) {
	683	_beta = val / k + betac;
	684	};
	685	/*! Create Dirichlet random walk
	686	\f[ f(rv\|rvc) = Di(rvc*k) \f]
	687	from structure
	688	\code
	689	class = 'mDirich';
	690	k = 1; // multiplicative constant k
	691	--- optional ---
	692	rv = RV({'name'},size) // description of RV
	693	beta0 = []; // initial value of beta
	694	betac = []; // initial value of beta
	695	\endcode
	696	*/
	697	void from_setting ( const Setting &set );
	698	void to_setting (Setting &set) const;
	699	void validate();
688	700	};
689	701	UIREGISTER ( mDirich );
…	…
691	703	/*! \brief Random Walk with vector Beta distribution
692	704	Using simple assignment
693		\f{eqnarray*}
	705	\f{eqnarray*}
694	706	\alpha & = & rvc / k + \beta_c \\
695	707	\beta & = &(1-rvc) / k + \beta_c \\
…	…
702	714	By default is it set to 0.1;
703	715	*/
704		class mBeta: public pdf_internal<eBeta>{
705		//! vector of constants \f$ k \f$ of the random walk
706		vec k;
707		//! stabilizing coefficient \f$ \beta_c \f$
708		vec betac;
709
710		public:
711		void condition ( const vec &val ) {
712		this->iepdf.alpha = elem_div(val , k) + betac;
713		this->iepdf.beta = elem_div (1-val , k) + betac;
714		};
715
716		/*! Create Beta random walk
717		\f[ f(rv\|rvc) = \prod Beta(rvc,k) \f]
718		from structure
719		\code
720		class = 'mBeta';
721		k = 1; // multiplicative constant k
722		--- optional ---
723		rv = RV({'name'},size) // description of RV
724		beta = []; // initial value of beta
725		betac = []; // initial value of beta stabilizing constant
726		\endcode
727		*/
728		void from_setting ( const Setting &set );
729
730		void to_setting (Setting &set) const;
731
732		void validate(){
733		pdf_internal<eBeta>::validate();
734		bdm_assert(betac.length()==dimension(),"Incomaptible betac");
735		bdm_assert(k.length()==dimension(),"Incomaptible k");
736		dimc = iepdf.dimension();
737		}
738		//!
	716	class mBeta: public pdf_internal<eBeta> {
	717	//! vector of constants \f$ k \f$ of the random walk
	718	vec k;
	719	//! stabilizing coefficient \f$ \beta_c \f$
	720	vec betac;
	721
	722	public:
	723	void condition ( const vec &val ) {
	724	this->iepdf.alpha = elem_div(val , k) + betac;
	725	this->iepdf.beta = elem_div (1-val , k) + betac;
	726	};
	727
	728	/*! Create Beta random walk
	729	\f[ f(rv\|rvc) = \prod Beta(rvc,k) \f]
	730	from structure
	731	\code
	732	class = 'mBeta';
	733	k = 1; // multiplicative constant k
	734	--- optional ---
	735	rv = RV({'name'},size) // description of RV
	736	beta = []; // initial value of beta
	737	betac = []; // initial value of beta stabilizing constant
	738	\endcode
	739	*/
	740	void from_setting ( const Setting &set );
	741
	742	void to_setting (Setting &set) const;
	743
	744	void validate() {
	745	pdf_internal<eBeta>::validate();
	746	bdm_assert(betac.length()==dimension(),"Incomaptible betac");
	747	bdm_assert(k.length()==dimension(),"Incomaptible k");
	748	dimc = iepdf.dimension();
	749	}
	750	//!
739	751	};
740	752	UIREGISTER(mBeta);
…	…
743	755	class multiBM : public BMEF {
744	756	protected:
745		//! Conjugate prior and posterior
746		eDirich est;
747		//! Pointer inside est to sufficient statistics
748		vec β
749		public:
750		//!Default constructor
751		multiBM () : BMEF (), est (), beta ( est._beta() ) {
752		if ( beta.length() > 0 ) {
753		last_lognc = est.lognc();
754		} else {
755		last_lognc = 0.0;
756		}
757		}
758		//!Copy constructor
759		multiBM ( const multiBM &B ) : BMEF ( B ), est ( B.est ), beta ( est._beta() ) {}
760		//! Sets sufficient statistics to match that of givefrom mB0
761		void set_statistics ( const BM* mB0 ) {
762		const multiBM* mB = dynamic_cast<const multiBM*> ( mB0 );
763		beta = mB->beta;
764		}
765		void bayes ( const vec &yt, const vec &cond = empty_vec );
766
767		double logpred ( const vec &yt ) const;
768
769		void flatten ( const BMEF* B , double weight);
770
771		//! return correctly typed posterior (covariant return)
772		const eDirich& posterior() const {
773		return est;
774		};
775		//! constructor function
776		void set_parameters ( const vec &beta0 ) {
777		est.set_parameters ( beta0 );
778		est.validate();
779		if ( evalll ) {
780		last_lognc = est.lognc();
781		}
782		}
783
784		void to_setting ( Setting &set ) const {
785		BMEF::to_setting ( set );
786		UI::save( &est, set, "prior" );
787		}
788		void from_setting (const Setting &set ) {
789		BMEF::from_setting ( set );
790		UI::get( est, set, "prior" );
791		}
	757	//! Conjugate prior and posterior
	758	eDirich est;
	759	//! Pointer inside est to sufficient statistics
	760	vec β
	761	public:
	762	//!Default constructor
	763	multiBM () : BMEF (), est (), beta ( est._beta() ) {
	764	if ( beta.length() > 0 ) {
	765	last_lognc = est.lognc();
	766	} else {
	767	last_lognc = 0.0;
	768	}
	769	}
	770	//!Copy constructor
	771	multiBM ( const multiBM &B ) : BMEF ( B ), est ( B.est ), beta ( est._beta() ) {}
	772	//! Sets sufficient statistics to match that of givefrom mB0
	773	void set_statistics ( const BM* mB0 ) {
	774	const multiBM* mB = dynamic_cast<const multiBM*> ( mB0 );
	775	beta = mB->beta;
	776	}
	777	void bayes ( const vec &yt, const vec &cond = empty_vec );
	778
	779	double logpred ( const vec &yt ) const;
	780
	781	void flatten ( const BMEF* B , double weight);
	782
	783	//! return correctly typed posterior (covariant return)
	784	const eDirich& posterior() const {
	785	return est;
	786	};
	787	//! constructor function
	788	void set_parameters ( const vec &beta0 ) {
	789	est.set_parameters ( beta0 );
	790	est.validate();
	791	if ( evalll ) {
	792	last_lognc = est.lognc();
	793	}
	794	}
	795
	796	void to_setting ( Setting &set ) const {
	797	BMEF::to_setting ( set );
	798	UI::save( &est, set, "prior" );
	799	}
	800	void from_setting (const Setting &set ) {
	801	BMEF::from_setting ( set );
	802	UI::get( est, set, "prior" );
	803	}
792	804	};
793	805	UIREGISTER( multiBM );
…	…
804	816	class egamma : public eEF {
805	817	protected:
806		//! Vector \f$\alpha\f$
807		vec alpha;
808		//! Vector \f$\beta\f$
809		vec beta;
	818	//! Vector \f$\alpha\f$
	819	vec alpha;
	820	//! Vector \f$\beta\f$
	821	vec beta;
810	822	public :
811		//! \name Constructors
812		//!@{
813		egamma () : eEF (), alpha ( 0 ), beta ( 0 ) {};
814		egamma ( const vec &a, const vec &b ) {
815		set_parameters ( a, b );
816		validate();
817		};
818		void set_parameters ( const vec &a, const vec &b ) {
819		alpha = a, beta = b;
820		};
821		//!@}
822
823		vec sample() const;
824		double evallog ( const vec &val ) const;
825		double lognc () const;
826		//! Returns pointer to internal alpha. Potentially dengerous: use with care!
827		vec& _alpha() {
828		return alpha;
829		}
830		//! Returns pointer to internal beta. Potentially dengerous: use with care!
831		vec& _beta() {
832		return beta;
833		}
834		vec mean() const {
835		return elem_div ( alpha, beta );
836		}
837		vec variance() const {
838		return elem_div ( alpha, elem_mult ( beta, beta ) );
839		}
840
841		/*! Create object from the following structure
842
843		\code
844		class = 'egamma';
845		alpha = [...]; % vector of alpha
846		beta = [...]; % vector of beta
847		--- inherited fields ---
848		bdm::eEF::from_setting
849		\endcode
850		fulfilling formula \f[ f(rv\|rvc) = \Gamma(\alpha, \beta) \f]
851		*/
852		void from_setting ( const Setting &set );
853
854		void to_setting ( Setting &set ) const;
855		void validate();
	823	//! \name Constructors
	824	//!@{
	825	egamma () : eEF (), alpha ( 0 ), beta ( 0 ) {};
	826	egamma ( const vec &a, const vec &b ) {
	827	set_parameters ( a, b );
	828	validate();
	829	};
	830	void set_parameters ( const vec &a, const vec &b ) {
	831	alpha = a, beta = b;
	832	};
	833	//!@}
	834
	835	vec sample() const;
	836	double evallog ( const vec &val ) const;
	837	double lognc () const;
	838	//! Returns pointer to internal alpha. Potentially dengerous: use with care!
	839	vec& _alpha() {
	840	return alpha;
	841	}
	842	//! Returns pointer to internal beta. Potentially dengerous: use with care!
	843	vec& _beta() {
	844	return beta;
	845	}
	846	vec mean() const {
	847	return elem_div ( alpha, beta );
	848	}
	849	vec variance() const {
	850	return elem_div ( alpha, elem_mult ( beta, beta ) );
	851	}
	852
	853	/*! Create object from the following structure
	854
	855	\code
	856	class = 'egamma';
	857	alpha = [...]; % vector of alpha
	858	beta = [...]; % vector of beta
	859	--- inherited fields ---
	860	bdm::eEF::from_setting
	861	\endcode
	862	fulfilling formula \f[ f(rv\|rvc) = \Gamma(\alpha, \beta) \f]
	863	*/
	864	void from_setting ( const Setting &set );
	865
	866	void to_setting ( Setting &set ) const;
	867	void validate();
856	868	};
857	869	UIREGISTER ( egamma );
…	…
877	889	protected:
878	890	public :
879		//! \name Constructors
880		//! All constructors are inherited
881		//!@{
882		//!@}
883
884		vec sample() const {
885		return 1.0 / egamma::sample();
886		};
887		//! Returns poiter to alpha and beta. Potentially dangerous: use with care!
888		vec mean() const {
889		return elem_div ( beta, alpha - 1 );
890		}
891		vec variance() const {
892		vec mea = mean();
893		return elem_div ( elem_mult ( mea, mea ), alpha - 2 );
894		}
	891	//! \name Constructors
	892	//! All constructors are inherited
	893	//!@{
	894	//!@}
	895
	896	vec sample() const {
	897	return 1.0 / egamma::sample();
	898	};
	899	//! Returns poiter to alpha and beta. Potentially dangerous: use with care!
	900	vec mean() const {
	901	return elem_div ( beta, alpha - 1 );
	902	}
	903	vec variance() const {
	904	vec mea = mean();
	905	return elem_div ( elem_mult ( mea, mea ), alpha - 2 );
	906	}
895	907	};
896	908	/*
…	…
914	926	protected:
915	927	//! lower bound on support
916		vec low;
	928	vec low;
917	929	//! upper bound on support
918		vec high;
	930	vec high;
919	931	//! internal
920		vec distance;
	932	vec distance;
921	933	//! normalizing coefficients
922		double nk;
	934	double nk;
923	935	//! cache of log( \c nk )
924		double lnk;
925		public:
926		//! \name Constructors
927		//!@{
928		euni () : epdf () {}
929		euni ( const vec &low0, const vec &high0 ) {
930		set_parameters ( low0, high0 );
931		}
932		void set_parameters ( const vec &low0, const vec &high0 ) {
933		distance = high0 - low0;
934		low = low0;
935		high = high0;
936		nk = prod ( 1.0 / distance );
937		lnk = log ( nk );
938		}
939		//!@}
940
941		double evallog ( const vec &val ) const {
942		if ( any ( val < low ) && any ( val > high ) ) {
943		return -inf;
944		} else return lnk;
945		}
946		vec sample() const {
947		vec smp ( dim );
	936	double lnk;
	937	public:
	938	//! \name Constructors
	939	//!@{
	940	euni () : epdf () {}
	941	euni ( const vec &low0, const vec &high0 ) {
	942	set_parameters ( low0, high0 );
	943	}
	944	void set_parameters ( const vec &low0, const vec &high0 ) {
	945	distance = high0 - low0;
	946	low = low0;
	947	high = high0;
	948	nk = prod ( 1.0 / distance );
	949	lnk = log ( nk );
	950	}
	951	//!@}
	952
	953	double evallog ( const vec &val ) const {
	954	if ( any ( val < low ) && any ( val > high ) ) {
	955	return -inf;
	956	} else return lnk;
	957	}
	958	vec sample() const {
	959	vec smp ( dim );
948	960	#pragma omp critical
949		UniRNG.sample_vector ( dim , smp );
950		return low + elem_mult ( distance, smp );
951		}
952		//! set values of \c low and \c high
953		vec mean() const {
954		return ( high - low ) / 2.0;
955		}
956		vec variance() const {
957		return ( pow ( high, 2 ) + pow ( low, 2 ) + elem_mult ( high, low ) ) / 3.0;
958		}
959		/*! Create Uniform density
960		\f[ f(rv) = U(low,high) \f]
961		from structure
962		\code
963		class = 'euni'
964		high = [...]; // vector of upper bounds
965		low = [...]; // vector of lower bounds
966		rv = RV({'name'}); // description of RV
967		\endcode
968		*/
969		void from_setting ( const Setting &set );
970		void to_setting (Setting &set) const;
971		void validate();
	961	UniRNG.sample_vector ( dim , smp );
	962	return low + elem_mult ( distance, smp );
	963	}
	964	//! set values of \c low and \c high
	965	vec mean() const {
	966	return ( high - low ) / 2.0;
	967	}
	968	vec variance() const {
	969	return ( pow ( high, 2 ) + pow ( low, 2 ) + elem_mult ( high, low ) ) / 3.0;
	970	}
	971	/*! Create Uniform density
	972	\f[ f(rv) = U(low,high) \f]
	973	from structure
	974	\code
	975	class = 'euni'
	976	high = [...]; // vector of upper bounds
	977	low = [...]; // vector of lower bounds
	978	rv = RV({'name'}); // description of RV
	979	\endcode
	980	*/
	981	void from_setting ( const Setting &set );
	982	void to_setting (Setting &set) const;
	983	void validate();
972	984	};
973	985	UIREGISTER ( euni );
…	…
975	987	//! Uniform density with conditional mean value
976	988	class mguni : public pdf_internal<euni> {
977		//! function of the mean value
978		shared_ptr<fnc> mean;
979		//! distance from mean to both sides
980		vec delta;
981		public:
982		void condition ( const vec &cond ) {
983		vec mea = mean->eval ( cond );
984		iepdf.set_parameters ( mea - delta, mea + delta );
985		}
986		//! load from
987		void from_setting ( const Setting &set ) {
988		pdf::from_setting ( set ); //reads rv and rvc
989		UI::get ( delta, set, "delta", UI::compulsory );
990		mean = UI::build<fnc> ( set, "mean", UI::compulsory );
991		iepdf.set_parameters ( -delta, delta );
992		}
993		void to_setting (Setting &set) const {
994		pdf::to_setting ( set );
995		UI::save( iepdf.mean(), set, "delta");
996		UI::save(mean, set, "mean");
997		}
998		~~void validate()~~{
999		pdf_internal<euni>::validate();
1000		dimc = mean->dimensionc();
1001
1002		}
	989	//! function of the mean value
	990	shared_ptr<fnc> mean;
	991	//! distance from mean to both sides
	992	vec delta;
	993	public:
	994	void condition ( const vec &cond ) {
	995	vec mea = mean->eval ( cond );
	996	iepdf.set_parameters ( mea - delta, mea + delta );
	997	}
	998	//! load from
	999	void from_setting ( const Setting &set ) {
	1000	pdf::from_setting ( set ); //reads rv and rvc
	1001	UI::get ( delta, set, "delta", UI::compulsory );
	1002	mean = UI::build<fnc> ( set, "mean", UI::compulsory );
	1003	iepdf.set_parameters ( -delta, delta );
	1004	}
	1005	void to_setting (Setting &set) const {
	1006	pdf::to_setting ( set );
	1007	UI::save( iepdf.mean(), set, "delta");
	1008	UI::save(mean, set, "mean");
	1009	}
	1010	void validate() {
	1011	pdf_internal<euni>::validate();
	1012	dimc = mean->dimensionc();
	1013
	1014	}
1003	1015
1004	1016	};
…	…
1012	1024	class mlnorm : public pdf_internal< TEpdf<sq_T> > {
1013	1025	protected:
1014		//! Internal epdf that arise by conditioning on \c rvc
1015		mat A;
1016		//! Constant additive term
1017		vec mu_const;
	1026	//! Internal epdf that arise by conditioning on \c rvc
	1027	mat A;
	1028	//! Constant additive term
	1029	vec mu_const;
1018	1030	// vec& _mu; //cached epdf.mu; !!!!!! WHY NOT?
1019	1031	public:
1020		//! \name Constructors
1021		//!@{
1022		mlnorm() : pdf_internal< TEpdf<sq_T> >() {};
1023		mlnorm ( const mat &A, const vec &mu0, const sq_T &R ) : pdf_internal< TEpdf<sq_T> >() {
1024		set_parameters ( A, mu0, R );
1025		validate();
1026		}
1027
1028		//! Set \c A and \c R
1029		void set_parameters ( const mat &A0, const vec &mu0, const sq_T &R0 ) {
1030		this->iepdf.set_parameters ( zeros ( A0.rows() ), R0 );
1031		A = A0;
1032		mu_const = mu0;
1033		}
1034
1035		//!@}
1036		//! Set value of \c rvc . Result of this operation is stored in \c epdf use function \c _ep to access it.
1037		void condition ( const vec &cond ) {
1038		this->iepdf._mu() = A * cond + mu_const;
	1032	//! \name Constructors
	1033	//!@{
	1034	mlnorm() : pdf_internal< TEpdf<sq_T> >() {};
	1035	mlnorm ( const mat &A, const vec &mu0, const sq_T &R ) : pdf_internal< TEpdf<sq_T> >() {
	1036	set_parameters ( A, mu0, R );
	1037	validate();
	1038	}
	1039
	1040	//! Set \c A and \c R
	1041	void set_parameters ( const mat &A0, const vec &mu0, const sq_T &R0 ) {
	1042	this->iepdf.set_parameters ( zeros ( A0.rows() ), R0 );
	1043	A = A0;
	1044	mu_const = mu0;
	1045	}
	1046
	1047	//!@}
	1048	//! Set value of \c rvc . Result of this operation is stored in \c epdf use function \c _ep to access it.
	1049	void condition ( const vec &cond ) {
	1050	this->iepdf._mu() = A * cond + mu_const;
1039	1051	//R is already assigned;
1040		}
1041
1042		//!access function
1043		const vec& _mu_const() const {
1044		return mu_const;
1045		}
1046		//!access function
1047		const mat& _A() const {
1048		return A;
1049		}
1050		//!access function
1051		mat _R() const {
1052		return this->iepdf._R().to_mat();
1053		}
1054		//!access function
1055		sq_T __R() const {
1056		return this->iepdf._R();
1057		}
1058
1059		//! Debug stream
1060		template<typename sq_M>
1061		friend std::ostream &operator<< ( std::ostream &os, mlnorm<sq_M, enorm> &ml );
1062
1063		/*! Create Normal density with linear function of mean value
1064		\f[ f(rv\|rvc) = N(A*rvc+const, R) \f]
1065		from structure
1066		\code
1067		class = 'mlnorm<ldmat>', (OR) 'mlnorm<chmat>', (OR) 'mlnorm<fsqmat>';
1068		A = []; // matrix or vector of appropriate dimension
1069		R = []; // square matrix of appropriate dimension
1070		--- optional ---
1071		const = zeros(A.rows); // vector of constant term
1072		\endcode
1073		*/
1074		void from_setting ( const Setting &set ) {
1075		pdf::from_setting ( set );
1076
1077		UI::get ( A, set, "A", UI::compulsory );
1078		UI::get ( mu_const, set, "const", UI::optional);
1079		mat R0;
1080		UI::get ( R0, set, "R", UI::compulsory );
1081		set_parameters ( A, mu_const, R0 );
1082		}
1083
1084		void to_setting (Setting &set) const {
1085		pdf::to_setting(set);
1086		UI::save ( A, set, "A");
1087		UI::save ( mu_const, set, "const");
1088		UI::save ( _R(), set, "R");
1089		}
1090
1091		void validate() {
1092		pdf_internal<TEpdf<sq_T> >::validate();
1093		if (mu_const.length()==0) { // default in from_setting
1094		mu_const=zeros(A.rows());
1095		}
1096		bdm_assert ( A.rows() == mu_const.length(), "mlnorm: A vs. mu mismatch" );
1097		bdm_assert ( A.rows() == _R().rows(), "mlnorm: A vs. R mismatch" );
1098		this->dimc = A.cols();
1099
1100		}
	1052	}
	1053
	1054	//!access function
	1055	const vec& _mu_const() const {
	1056	return mu_const;
	1057	}
	1058	//!access function
	1059	const mat& _A() const {
	1060	return A;
	1061	}
	1062	//!access function
	1063	mat _R() const {
	1064	return this->iepdf._R().to_mat();
	1065	}
	1066	//!access function
	1067	sq_T __R() const {
	1068	return this->iepdf._R();
	1069	}
	1070
	1071	//! Debug stream
	1072	template<typename sq_M>
	1073	friend std::ostream &operator<< ( std::ostream &os, mlnorm<sq_M, enorm> &ml );
	1074
	1075	/*! Create Normal density with linear function of mean value
	1076	\f[ f(rv\|rvc) = N(A*rvc+const, R) \f]
	1077	from structure
	1078	\code
	1079	class = 'mlnorm<ldmat>', (OR) 'mlnorm<chmat>', (OR) 'mlnorm<fsqmat>';
	1080	A = []; // matrix or vector of appropriate dimension
	1081	R = []; // square matrix of appropriate dimension
	1082	--- optional ---
	1083	const = zeros(A.rows); // vector of constant term
	1084	\endcode
	1085	*/
	1086	void from_setting ( const Setting &set ) {
	1087	pdf::from_setting ( set );
	1088
	1089	UI::get ( A, set, "A", UI::compulsory );
	1090	UI::get ( mu_const, set, "const", UI::optional);
	1091	mat R0;
	1092	UI::get ( R0, set, "R", UI::compulsory );
	1093	set_parameters ( A, mu_const, R0 );
	1094	}
	1095
	1096	void to_setting (Setting &set) const {
	1097	pdf::to_setting(set);
	1098	UI::save ( A, set, "A");
	1099	UI::save ( mu_const, set, "const");
	1100	UI::save ( _R(), set, "R");
	1101	}
	1102
	1103	void validate() {
	1104	pdf_internal<TEpdf<sq_T> >::validate();
	1105	if (mu_const.length()==0) { // default in from_setting
	1106	mu_const=zeros(A.rows());
	1107	}
	1108	bdm_assert ( A.rows() == mu_const.length(), "mlnorm: A vs. mu mismatch" );
	1109	bdm_assert ( A.rows() == _R().rows(), "mlnorm: A vs. R mismatch" );
	1110	this->dimc = A.cols();
	1111
	1112	}
1101	1113	};
1102	1114	UIREGISTER2 ( mlnorm, ldmat );
…	…
1107	1119	SHAREDPTR2 ( mlnorm, chmat );
1108	1120
1109		//! \class mlgauss
	1121	//! \class mlgauss
1110	1122	//!\brief Normal distribution with linear function of mean value. Same as mlnorm<fsqmat>.
1111	1123	typedef mlnorm<fsqmat> mlgauss;
…	…
1117	1129	private:
1118	1130	// vec μ WHY NOT?
1119		shared_ptr<fnc> g;
1120
1121		public:
1122		//!default constructor
1123		mgnorm() : pdf_internal<enorm<sq_T> >() { }
1124		//!set mean function
1125		inline void set_parameters ( const shared_ptr<fnc> &g0, const sq_T &R0 );
1126		inline void condition ( const vec &cond );
1127
1128
1129		/*! Create Normal density with given function of mean value
1130		\f[ f(rv\|rvc) = N( g(rvc), R) \f]
1131		from structure
1132		\code
1133		class = 'mgnorm';
1134		g.class = 'fnc'; // function for mean value evolution
1135		g._fields_of_fnc = ...;
1136
1137		R = [1, 0; // covariance matrix
1138		0, 1];
1139		--OR --
1140		dR = [1, 1]; // diagonal of cavariance matrix
1141
1142		rv = RV({'name'}) // description of RV
1143		rvc = RV({'name'}) // description of RV in condition
1144		\endcode
1145		*/
1146
1147
1148		void from_setting ( const Setting &set ) {
1149		pdf::from_setting ( set );
1150		shared_ptr<fnc> g = UI::build<fnc> ( set, "g", UI::compulsory );
1151
1152		mat R;
1153		vec dR;
1154		if ( UI::get ( dR, set, "dR" ) )
1155		R = diag ( dR );
1156		else
1157		UI::get ( R, set, "R", UI::compulsory );
1158
1159		set_parameters ( g, R );
1160		//validate();
1161		}
1162
1163
1164		void to_setting (Setting &set) const {
1165		UI::save( g,set, "g");
1166		UI::save(this->iepdf._R().to_mat(),set, "R");
1167
1168		}
1169
1170
1171
1172		void validate() {
1173		this->iepdf.validate();
1174		bdm_assert ( g->dimension() == this->iepdf.dimension(), "incompatible function" );
1175		this->dim = g->dimension();
1176		this->dimc = g->dimensionc();
1177		this->iepdf.validate();
1178		}
	1131	shared_ptr<fnc> g;
	1132
	1133	public:
	1134	//!default constructor
	1135	mgnorm() : pdf_internal<enorm<sq_T> >() { }
	1136	//!set mean function
	1137	inline void set_parameters ( const shared_ptr<fnc> &g0, const sq_T &R0 );
	1138	inline void condition ( const vec &cond );
	1139
	1140
	1141	/*! Create Normal density with given function of mean value
	1142	\f[ f(rv\|rvc) = N( g(rvc), R) \f]
	1143	from structure
	1144	\code
	1145	class = 'mgnorm';
	1146	g.class = 'fnc'; // function for mean value evolution
	1147	g._fields_of_fnc = ...;
	1148
	1149	R = [1, 0; // covariance matrix
	1150	0, 1];
	1151	--OR --
	1152	dR = [1, 1]; // diagonal of cavariance matrix
	1153
	1154	rv = RV({'name'}) // description of RV
	1155	rvc = RV({'name'}) // description of RV in condition
	1156	\endcode
	1157	*/
	1158
	1159
	1160	void from_setting ( const Setting &set ) {
	1161	pdf::from_setting ( set );
	1162	shared_ptr<fnc> g = UI::build<fnc> ( set, "g", UI::compulsory );
	1163
	1164	mat R;
	1165	vec dR;
	1166	if ( UI::get ( dR, set, "dR" ) )
	1167	R = diag ( dR );
	1168	else
	1169	UI::get ( R, set, "R", UI::compulsory );
	1170
	1171	set_parameters ( g, R );
	1172	//validate();
	1173	}
	1174
	1175
	1176	void to_setting (Setting &set) const {
	1177	UI::save( g,set, "g");
	1178	UI::save(this->iepdf._R().to_mat(),set, "R");
	1179
	1180	}
	1181
	1182
	1183
	1184	void validate() {
	1185	this->iepdf.validate();
	1186	bdm_assert ( g->dimension() == this->iepdf.dimension(), "incompatible function" );
	1187	this->dim = g->dimension();
	1188	this->dimc = g->dimensionc();
	1189	this->iepdf.validate();
	1190	}
1179	1191
1180	1192	};
…	…
1194	1206	class mlstudent : public mlnorm<ldmat, enorm> {
1195	1207	protected:
1196		//! Variable \f$ \Lambda \f$ from theory
1197		ldmat Lambda;
1198		//! Reference to variable \f$ R \f$
1199		ldmat &_R;
1200		//! Variable \f$ R_e \f$
1201		ldmat Re;
1202		public:
1203		mlstudent () : mlnorm<ldmat, enorm> (),
1204		Lambda (), _R ( iepdf._R() ) {}
1205		//! constructor function
1206		void set_parameters ( const mat &A0, const vec &mu0, const ldmat &R0, const ldmat& Lambda0 ) {
1207		iepdf.set_parameters ( mu0, R0 );// was Lambda, why?
1208		A = A0;
1209		mu_const = mu0;
1210		Re = R0;
1211		Lambda = Lambda0;
1212		}
1213
1214		void condition ( const vec &cond );
1215
1216		void validate() {
1217		mlnorm<ldmat, enorm>::validate();
1218		bdm_assert ( A.rows() == mu_const.length(), "mlstudent: A vs. mu mismatch" );
1219		bdm_assert ( _R.rows() == A.rows(), "mlstudent: A vs. R mismatch" );
1220
1221		}
	1208	//! Variable \f$ \Lambda \f$ from theory
	1209	ldmat Lambda;
	1210	//! Reference to variable \f$ R \f$
	1211	ldmat &_R;
	1212	//! Variable \f$ R_e \f$
	1213	ldmat Re;
	1214	public:
	1215	mlstudent () : mlnorm<ldmat, enorm> (),
	1216	Lambda (), _R ( iepdf._R() ) {}
	1217	//! constructor function
	1218	void set_parameters ( const mat &A0, const vec &mu0, const ldmat &R0, const ldmat& Lambda0 ) {
	1219	iepdf.set_parameters ( mu0, R0 );// was Lambda, why?
	1220	A = A0;
	1221	mu_const = mu0;
	1222	Re = R0;
	1223	Lambda = Lambda0;
	1224	}
	1225
	1226	void condition ( const vec &cond );
	1227
	1228	void validate() {
	1229	mlnorm<ldmat, enorm>::validate();
	1230	bdm_assert ( A.rows() == mu_const.length(), "mlstudent: A vs. mu mismatch" );
	1231	bdm_assert ( _R.rows() == A.rows(), "mlstudent: A vs. R mismatch" );
	1232
	1233	}
1222	1234	};
1223	1235
…	…
1234	1246	protected:
1235	1247
1236		//! Constant \f$k\f$
1237		double k;
1238
1239		//! cache of iepdf.beta
1240		vec &_beta;
1241
1242		public:
1243		//! Constructor
1244		mgamma() : pdf_internal<egamma>(), k ( 0 ),
1245		_beta ( iepdf._beta() ) {
1246		}
1247
1248		//! Set value of \c k
1249		void set_parameters ( double k, const vec &beta0 );
1250
1251		void condition ( const vec &val ) {
1252		_beta = k / val;
1253		};
1254		/*! Create Gamma density with conditional mean value
1255		\f[ f(rv\|rvc) = \Gamma(k, k/rvc) \f]
1256		from structure
1257		\code
1258		class = 'mgamma';
1259		beta = [...]; // vector of initial alpha
1260		k = 1.1; // multiplicative constant k
1261		rv = RV({'name'}) // description of RV
1262		rvc = RV({'name'}) // description of RV in condition
1263		\endcode
1264		*/
1265		void from_setting ( const Setting &set );
1266		void to_setting (Setting &set) const;
1267		void validate();
	1248	//! Constant \f$k\f$
	1249	double k;
	1250
	1251	//! cache of iepdf.beta
	1252	vec &_beta;
	1253
	1254	public:
	1255	//! Constructor
	1256	mgamma() : pdf_internal<egamma>(), k ( 0 ),
	1257	_beta ( iepdf._beta() ) {
	1258	}
	1259
	1260	//! Set value of \c k
	1261	void set_parameters ( double k, const vec &beta0 );
	1262
	1263	void condition ( const vec &val ) {
	1264	_beta = k / val;
	1265	};
	1266	/*! Create Gamma density with conditional mean value
	1267	\f[ f(rv\|rvc) = \Gamma(k, k/rvc) \f]
	1268	from structure
	1269	\code
	1270	class = 'mgamma';
	1271	beta = [...]; // vector of initial alpha
	1272	k = 1.1; // multiplicative constant k
	1273	rv = RV({'name'}) // description of RV
	1274	rvc = RV({'name'}) // description of RV in condition
	1275	\endcode
	1276	*/
	1277	void from_setting ( const Setting &set );
	1278	void to_setting (Setting &set) const;
	1279	void validate();
1268	1280	};
1269	1281	UIREGISTER ( mgamma );
…	…
1281	1293	class migamma : public pdf_internal<eigamma> {
1282	1294	protected:
1283		//! Constant \f$k\f$
1284		double k;
1285
1286		//! cache of iepdf.alpha
1287		vec &_alpha;
1288
1289		//! cache of iepdf.beta
1290		vec &_beta;
1291
1292		public:
1293		//! \name Constructors
1294		//!@{
1295		migamma() : pdf_internal<eigamma>(),
1296		k ( 0 ),
1297		_alpha ( iepdf._alpha() ),
1298		_beta ( iepdf._beta() ) {
1299		}
1300
1301		migamma ( const migamma &m ) : pdf_internal<eigamma>(),
1302		k ( 0 ),
1303		_alpha ( iepdf._alpha() ),
1304		_beta ( iepdf._beta() ) {
1305		}
1306		//!@}
1307
1308		//! Set value of \c k
1309		void set_parameters ( int len, double k0 ) {
1310		k = k0;
1311		iepdf.set_parameters ( ( 1.0 / ( kk ) + 2.0 ) ones ( len ) /alpha/, ones ( len ) /beta/ );
1312		};
1313
1314		~~void validate ()~~{
1315		pdf_internal<eigamma>::validate();
1316		dimc = dimension();
1317		};
1318
1319		void condition ( const vec &val ) {
1320		_beta = elem_mult ( val, ( _alpha - 1.0 ) );
1321		};
	1295	//! Constant \f$k\f$
	1296	double k;
	1297
	1298	//! cache of iepdf.alpha
	1299	vec &_alpha;
	1300
	1301	//! cache of iepdf.beta
	1302	vec &_beta;
	1303
	1304	public:
	1305	//! \name Constructors
	1306	//!@{
	1307	migamma() : pdf_internal<eigamma>(),
	1308	k ( 0 ),
	1309	_alpha ( iepdf._alpha() ),
	1310	_beta ( iepdf._beta() ) {
	1311	}
	1312
	1313	migamma ( const migamma &m ) : pdf_internal<eigamma>(),
	1314	k ( 0 ),
	1315	_alpha ( iepdf._alpha() ),
	1316	_beta ( iepdf._beta() ) {
	1317	}
	1318	//!@}
	1319
	1320	//! Set value of \c k
	1321	void set_parameters ( int len, double k0 ) {
	1322	k = k0;
	1323	iepdf.set_parameters ( ( 1.0 / ( kk ) + 2.0 ) ones ( len ) /alpha/, ones ( len ) /beta/ );
	1324	};
	1325
	1326	void validate () {
	1327	pdf_internal<eigamma>::validate();
	1328	dimc = dimension();
	1329	};
	1330
	1331	void condition ( const vec &val ) {
	1332	_beta = elem_mult ( val, ( _alpha - 1.0 ) );
	1333	};
1322	1334	};
1323	1335
…	…
1336	1348	class mgamma_fix : public mgamma {
1337	1349	protected:
1338		//! parameter l
1339		double l;
1340		//! reference vector
1341		vec refl;
1342		public:
1343		//! Constructor
1344		mgamma_fix () : mgamma (), refl () {};
1345		//! Set value of \c k
1346		void set_parameters ( double k0 , vec ref0, double l0 ) {
1347		mgamma::set_parameters ( k0, ref0 );
1348		refl = pow ( ref0, 1.0 - l0 );
1349		l = l0;
1350		};
1351
1352		~~void validate ()~~{
1353		mgamma::validate();
1354		dimc = dimension();
1355		};
1356
1357		void condition ( const vec &val ) {
1358		vec mean = elem_mult ( refl, pow ( val, l ) );
1359		_beta = k / mean;
1360		};
	1350	//! parameter l
	1351	double l;
	1352	//! reference vector
	1353	vec refl;
	1354	public:
	1355	//! Constructor
	1356	mgamma_fix () : mgamma (), refl () {};
	1357	//! Set value of \c k
	1358	void set_parameters ( double k0 , vec ref0, double l0 ) {
	1359	mgamma::set_parameters ( k0, ref0 );
	1360	refl = pow ( ref0, 1.0 - l0 );
	1361	l = l0;
	1362	};
	1363
	1364	void validate () {
	1365	mgamma::validate();
	1366	dimc = dimension();
	1367	};
	1368
	1369	void condition ( const vec &val ) {
	1370	vec mean = elem_mult ( refl, pow ( val, l ) );
	1371	_beta = k / mean;
	1372	};
1361	1373	};
1362	1374
…	…
1376	1388	class migamma_ref : public migamma {
1377	1389	protected:
1378		//! parameter l
1379		double l;
1380		//! reference vector
1381		vec refl;
1382		public:
1383		//! Constructor
1384		migamma_ref () : migamma (), refl () {};
1385
1386		//! Set value of \c k
1387		void set_parameters ( double k0 , vec ref0, double l0 ) {
1388		migamma::set_parameters ( ref0.length(), k0 );
1389		refl = pow ( ref0, 1.0 - l0 );
1390		l = l0;
1391		};
1392
1393		~~void validate()~~{
1394		migamma::validate();
1395		dimc = dimension();
1396		};
1397
1398		void condition ( const vec &val ) {
1399		vec mean = elem_mult ( refl, pow ( val, l ) );
1400		migamma::condition ( mean );
1401		};
1402
1403
1404		/*! Create inverse-Gamma density with conditional mean value
1405		\f[ f(rv\|rvc) = i\Gamma(k, k/(rvc^l \circ ref^{(1-l)}) \f]
1406		from structure
1407		\code
1408		class = 'migamma_ref';
1409		ref = [1e-5; 1e-5; 1e-2 1e-3]; // reference vector
1410		l = 0.999; // constant l
1411		k = 0.1; // constant k
1412		rv = RV({'name'}) // description of RV
1413		rvc = RV({'name'}) // description of RV in condition
1414		\endcode
1415		*/
1416		void from_setting ( const Setting &set );
1417
1418		void to_setting (Setting &set) const;
	1390	//! parameter l
	1391	double l;
	1392	//! reference vector
	1393	vec refl;
	1394	public:
	1395	//! Constructor
	1396	migamma_ref () : migamma (), refl () {};
	1397
	1398	//! Set value of \c k
	1399	void set_parameters ( double k0 , vec ref0, double l0 ) {
	1400	migamma::set_parameters ( ref0.length(), k0 );
	1401	refl = pow ( ref0, 1.0 - l0 );
	1402	l = l0;
	1403	};
	1404
	1405	void validate() {
	1406	migamma::validate();
	1407	dimc = dimension();
	1408	};
	1409
	1410	void condition ( const vec &val ) {
	1411	vec mean = elem_mult ( refl, pow ( val, l ) );
	1412	migamma::condition ( mean );
	1413	};
	1414
	1415
	1416	/*! Create inverse-Gamma density with conditional mean value
	1417	\f[ f(rv\|rvc) = i\Gamma(k, k/(rvc^l \circ ref^{(1-l)}) \f]
	1418	from structure
	1419	\code
	1420	class = 'migamma_ref';
	1421	ref = [1e-5; 1e-5; 1e-2 1e-3]; // reference vector
	1422	l = 0.999; // constant l
	1423	k = 0.1; // constant k
	1424	rv = RV({'name'}) // description of RV
	1425	rvc = RV({'name'}) // description of RV in condition
	1426	\endcode
	1427	*/
	1428	void from_setting ( const Setting &set );
	1429
	1430	void to_setting (Setting &set) const;
1419	1431	};
1420	1432
…	…
1435	1447	class elognorm: public enorm<ldmat> {
1436	1448	public:
1437		vec sample() const {
1438		return exp ( enorm<ldmat>::sample() );
1439		};
1440		vec mean() const {
1441		vec var = enorm<ldmat>::variance();
1442		return exp ( mu - 0.5*var );
1443		};
	1449	vec sample() const {
	1450	return exp ( enorm<ldmat>::sample() );
	1451	};
	1452	vec mean() const {
	1453	vec var = enorm<ldmat>::variance();
	1454	return exp ( mu - 0.5*var );
	1455	};
1444	1456
1445	1457	};
…	…
1453	1465	class mlognorm : public pdf_internal<elognorm> {
1454	1466	protected:
1455		//! parameter 1/2*sigma^2
1456		double sig2;
1457
1458		//! access
1459		vec μ
1460		public:
1461		//! Constructor
1462		mlognorm() : pdf_internal<elognorm>(),
1463		sig2 ( 0 ),
1464		mu ( iepdf._mu() ) {
1465		}
1466
1467		//! Set value of \c k
1468		void set_parameters ( int size, double k ) {
1469		sig2 = 0.5 * log ( k * k + 1 );
1470		iepdf.set_parameters ( zeros ( size ), 2sig2eye ( size ) );
1471		};
1472
1473		~~void validate()~~{
1474		pdf_internal<elognorm>::validate();
1475		dimc = iepdf.dimension();
1476		}
1477
1478		void condition ( const vec &val ) {
1479		mu = log ( val ) - sig2;//elem_mult ( refl,pow ( val,l ) );
1480		};
1481
1482		/*! Create logNormal random Walk
1483		\f[ f(rv\|rvc) = log\mathcal{N}( \log(rvc)-0.5\log(k^2+1), k I) \f]
1484		from structure
1485		\code
1486		class = 'mlognorm';
1487		k = 0.1; // "variance" k
1488		mu0 = 0.1; // Initial value of mean
1489		rv = RV({'name'}) // description of RV
1490		rvc = RV({'name'}) // description of RV in condition
1491		\endcode
1492		*/
1493		void from_setting ( const Setting &set );
1494
1495		void to_setting (Setting &set) const;
	1467	//! parameter 1/2*sigma^2
	1468	double sig2;
	1469
	1470	//! access
	1471	vec μ
	1472	public:
	1473	//! Constructor
	1474	mlognorm() : pdf_internal<elognorm>(),
	1475	sig2 ( 0 ),
	1476	mu ( iepdf._mu() ) {
	1477	}
	1478
	1479	//! Set value of \c k
	1480	void set_parameters ( int size, double k ) {
	1481	sig2 = 0.5 * log ( k * k + 1 );
	1482	iepdf.set_parameters ( zeros ( size ), 2sig2eye ( size ) );
	1483	};
	1484
	1485	void validate() {
	1486	pdf_internal<elognorm>::validate();
	1487	dimc = iepdf.dimension();
	1488	}
	1489
	1490	void condition ( const vec &val ) {
	1491	mu = log ( val ) - sig2;//elem_mult ( refl,pow ( val,l ) );
	1492	};
	1493
	1494	/*! Create logNormal random Walk
	1495	\f[ f(rv\|rvc) = log\mathcal{N}( \log(rvc)-0.5\log(k^2+1), k I) \f]
	1496	from structure
	1497	\code
	1498	class = 'mlognorm';
	1499	k = 0.1; // "variance" k
	1500	mu0 = 0.1; // Initial value of mean
	1501	rv = RV({'name'}) // description of RV
	1502	rvc = RV({'name'}) // description of RV in condition
	1503	\endcode
	1504	*/
	1505	void from_setting ( const Setting &set );
	1506
	1507	void to_setting (Setting &set) const;
1496	1508	};
1497	1509
…	…
1504	1516	class eWishartCh : public epdf {
1505	1517	protected:
1506		//! Upper-Triagle of Choleski decomposition of \f$ \Psi \f$
1507		chmat Y;
1508		//! dimension of matrix \f$ \Psi \f$
1509		int p;
1510		//! degrees of freedom \f$ \nu \f$
1511		double delta;
1512		public:
1513		//! Set internal structures
1514		void set_parameters ( const mat &Y0, const double delta0 ) {
1515		Y = chmat ( Y0 );
1516		delta = delta0;
1517		p = Y.rows();
1518		}
1519		//! Set internal structures
1520		void set_parameters ( const chmat &Y0, const double delta0 ) {
1521		Y = Y0;
1522		delta = delta0;
1523		p = Y.rows();
1524		}
1525
1526		~~virtual void validate ()~~{
1527		epdf::validate();
1528		dim = p * p;
1529		}
1530
1531		//! Sample matrix argument
1532		mat sample_mat() const {
1533		mat X = zeros ( p, p );
1534
1535		//sample diagonal
1536		for ( int i = 0; i < p; i++ ) {
1537		GamRNG.setup ( 0.5* ( delta - i ) , 0.5 ); // no +1 !! index if from 0
	1518	//! Upper-Triagle of Choleski decomposition of \f$ \Psi \f$
	1519	chmat Y;
	1520	//! dimension of matrix \f$ \Psi \f$
	1521	int p;
	1522	//! degrees of freedom \f$ \nu \f$
	1523	double delta;
	1524	public:
	1525	//! Set internal structures
	1526	void set_parameters ( const mat &Y0, const double delta0 ) {
	1527	Y = chmat ( Y0 );
	1528	delta = delta0;
	1529	p = Y.rows();
	1530	}
	1531	//! Set internal structures
	1532	void set_parameters ( const chmat &Y0, const double delta0 ) {
	1533	Y = Y0;
	1534	delta = delta0;
	1535	p = Y.rows();
	1536	}
	1537
	1538	virtual void validate () {
	1539	epdf::validate();
	1540	dim = p * p;
	1541	}
	1542
	1543	//! Sample matrix argument
	1544	mat sample_mat() const {
	1545	mat X = zeros ( p, p );
	1546
	1547	//sample diagonal
	1548	for ( int i = 0; i < p; i++ ) {
	1549	GamRNG.setup ( 0.5* ( delta - i ) , 0.5 ); // no +1 !! index if from 0
1538	1550	#pragma omp critical
1539		X ( i, i ) = sqrt ( GamRNG() );
1540		}
1541		//do the rest
1542		for ( int i = 0; i < p; i++ ) {
1543		for ( int j = i + 1; j < p; j++ ) {
	1551	X ( i, i ) = sqrt ( GamRNG() );
	1552	}
	1553	//do the rest
	1554	for ( int i = 0; i < p; i++ ) {
	1555	for ( int j = i + 1; j < p; j++ ) {
1544	1556	#pragma omp critical
1545		X ( i, j ) = NorRNG.sample();
1546		}
1547		}
1548		return X*Y._Ch();// return upper triangular part of the decomposition
1549		}
1550
1551		vec sample () const {
1552		return vec ( sample_mat()._data(), p*p );
1553		}
1554
1555		virtual vec mean() const NOT_IMPLEMENTED(0);
1556
1557		//! return expected variance (not covariance!)
1558		virtual vec variance() const NOT_IMPLEMENTED(0);
1559
1560		virtual double evallog ( const vec &val ) const NOT_IMPLEMENTED(0);
1561
1562		//! fast access function y0 will be copied into Y.Ch.
1563		void setY ( const mat &Ch0 ) {
1564		copy_vector ( dim, Ch0._data(), Y._Ch()._data() );
1565		}
1566
1567		//! fast access function y0 will be copied into Y.Ch.
1568		void _setY ( const vec &ch0 ) {
1569		copy_vector ( dim, ch0._data(), Y._Ch()._data() );
1570		}
1571
1572		//! access function
1573		const chmat& getY() const {
1574		return Y;
1575		}
	1557	X ( i, j ) = NorRNG.sample();
	1558	}
	1559	}
	1560	return X*Y._Ch();// return upper triangular part of the decomposition
	1561	}
	1562
	1563	vec sample () const {
	1564	return vec ( sample_mat()._data(), p*p );
	1565	}
	1566
	1567	virtual vec mean() const NOT_IMPLEMENTED(0);
	1568
	1569	//! return expected variance (not covariance!)
	1570	virtual vec variance() const NOT_IMPLEMENTED(0);
	1571
	1572	virtual double evallog ( const vec &val ) const NOT_IMPLEMENTED(0);
	1573
	1574	//! fast access function y0 will be copied into Y.Ch.
	1575	void setY ( const mat &Ch0 ) {
	1576	copy_vector ( dim, Ch0._data(), Y._Ch()._data() );
	1577	}
	1578
	1579	//! fast access function y0 will be copied into Y.Ch.
	1580	void _setY ( const vec &ch0 ) {
	1581	copy_vector ( dim, ch0._data(), Y._Ch()._data() );
	1582	}
	1583
	1584	//! access function
	1585	const chmat& getY() const {
	1586	return Y;
	1587	}
1576	1588	};
1577	1589
…	…
1581	1593	class eiWishartCh: public epdf {
1582	1594	protected:
1583		//! Internal instance of Wishart density
1584		eWishartCh W;
1585		//! size of Ch
1586		int p;
1587		//! parameter delta
1588		double delta;
1589		public:
1590		//! constructor function
1591		void set_parameters ( const mat &Y0, const double delta0 ) {
1592		delta = delta0;
1593		W.set_parameters ( inv ( Y0 ), delta0 );
1594		p = Y0.rows();
1595		}
1596
1597		~~virtual void validate ()~~{
1598		epdf::validate();
1599		W.validate();
1600		dim = W.dimension();
1601		}
1602
1603
1604		vec sample() const {
1605		mat iCh;
1606		iCh = inv ( W.sample_mat() );
1607		return vec ( iCh._data(), dim );
1608		}
1609		//! access function
1610		void _setY ( const vec &y0 ) {
1611		mat Ch ( p, p );
1612		mat iCh ( p, p );
1613		copy_vector ( dim, y0._data(), Ch._data() );
1614
1615		iCh = inv ( Ch );
1616		W.setY ( iCh );
1617		}
1618
1619		virtual double evallog ( const vec &val ) const {
1620		chmat X ( p );
1621		const chmat& Y = W.getY();
1622
1623		copy_vector ( p*p, val._data(), X._Ch()._data() );
1624		chmat iX ( p );
1625		X.inv ( iX );
1626		// compute
	1595	//! Internal instance of Wishart density
	1596	eWishartCh W;
	1597	//! size of Ch
	1598	int p;
	1599	//! parameter delta
	1600	double delta;
	1601	public:
	1602	//! constructor function
	1603	void set_parameters ( const mat &Y0, const double delta0 ) {
	1604	delta = delta0;
	1605	W.set_parameters ( inv ( Y0 ), delta0 );
	1606	p = Y0.rows();
	1607	}
	1608
	1609	virtual void validate () {
	1610	epdf::validate();
	1611	W.validate();
	1612	dim = W.dimension();
	1613	}
	1614
	1615
	1616	vec sample() const {
	1617	mat iCh;
	1618	iCh = inv ( W.sample_mat() );
	1619	return vec ( iCh._data(), dim );
	1620	}
	1621	//! access function
	1622	void _setY ( const vec &y0 ) {
	1623	mat Ch ( p, p );
	1624	mat iCh ( p, p );
	1625	copy_vector ( dim, y0._data(), Ch._data() );
	1626
	1627	iCh = inv ( Ch );
	1628	W.setY ( iCh );
	1629	}
	1630
	1631	virtual double evallog ( const vec &val ) const {
	1632	chmat X ( p );
	1633	const chmat& Y = W.getY();
	1634
	1635	copy_vector ( p*p, val._data(), X._Ch()._data() );
	1636	chmat iX ( p );
	1637	X.inv ( iX );
	1638	// compute
1627	1639	// \frac{ \|\Psi\|^{m/2}\|X\|^{-(m+p+1)/2}e^{-tr(\Psi X^{-1})/2} }{ 2^{mp/2}\Gamma_p(m/2)},
1628		mat M = Y.to_mat() * iX.to_mat();
1629
1630		double log1 = 0.5 * p * ( 2 * Y.logdet() ) - 0.5 * ( delta + p + 1 ) * ( 2 * X.logdet() ) - 0.5 * trace ( M );
1631		//Fixme! Multivariate gamma omitted!! it is ok for sampling, but not otherwise!!
1632
1633		/* if (0) {
1634		mat XX=X.to_mat();
1635		mat YY=Y.to_mat();
1636
1637		double log2 = 0.5plog(det(YY))-0.5(delta+p+1)log(det(XX))-0.5trace(YYinv(XX));
1638		cout << log1 << "," << log2 << endl;
1639		}*/
1640		return log1;
1641		};
1642
1643		virtual vec mean() const NOT_IMPLEMENTED(0);
1644
1645		//! return expected variance (not covariance!)
1646		virtual vec variance() const NOT_IMPLEMENTED(0);
	1640	mat M = Y.to_mat() * iX.to_mat();
	1641
	1642	double log1 = 0.5 * p * ( 2 * Y.logdet() ) - 0.5 * ( delta + p + 1 ) * ( 2 * X.logdet() ) - 0.5 * trace ( M );
	1643	//Fixme! Multivariate gamma omitted!! it is ok for sampling, but not otherwise!!
	1644
	1645	/* if (0) {
	1646	mat XX=X.to_mat();
	1647	mat YY=Y.to_mat();
	1648
	1649	double log2 = 0.5plog(det(YY))-0.5(delta+p+1)log(det(XX))-0.5trace(YYinv(XX));
	1650	cout << log1 << "," << log2 << endl;
	1651	}*/
	1652	return log1;
	1653	};
	1654
	1655	virtual vec mean() const NOT_IMPLEMENTED(0);
	1656
	1657	//! return expected variance (not covariance!)
	1658	virtual vec variance() const NOT_IMPLEMENTED(0);
1647	1659	};
1648	1660
…	…
1650	1662	class rwiWishartCh : public pdf_internal<eiWishartCh> {
1651	1663	protected:
1652		//!square root of \f$ \nu-p-1 \f$ - needed for computation of \f$ \Psi \f$ from conditions
1653		double sqd;
1654		//!reference point for diagonal
1655		vec refl;
1656		//! power of the reference
1657		double l;
1658		//! dimension
1659		int p;
1660
1661		public:
1662		rwiWishartCh() : sqd ( 0 ), l ( 0 ), p ( 0 ) {}
1663		//! constructor function
1664		void set_parameters ( int p0, double k, vec ref0, double l0 ) {
1665		p = p0;
1666		double delta = 2 / ( k * k ) + p + 3;
1667		sqd = sqrt ( delta - p - 1 );
1668		l = l0;
1669		refl = pow ( ref0, 1 - l );
1670		iepdf.set_parameters ( eye ( p ), delta );
1671		};
1672
1673		~~void validate()~~{
1674		pdf_internal<eiWishartCh>::validate();
1675		dimc = iepdf.dimension();
1676		}
1677
1678		void condition ( const vec &c ) {
1679		vec z = c;
1680		int ri = 0;
1681		for ( int i = 0; i < p*p; i += ( p + 1 ) ) {//trace diagonal element
1682		z ( i ) = pow ( z ( i ), l ) * refl ( ri );
1683		ri++;
1684		}
1685
1686		iepdf._setY ( sqd*z );
1687		}
	1664	//!square root of \f$ \nu-p-1 \f$ - needed for computation of \f$ \Psi \f$ from conditions
	1665	double sqd;
	1666	//!reference point for diagonal
	1667	vec refl;
	1668	//! power of the reference
	1669	double l;
	1670	//! dimension
	1671	int p;
	1672
	1673	public:
	1674	rwiWishartCh() : sqd ( 0 ), l ( 0 ), p ( 0 ) {}
	1675	//! constructor function
	1676	void set_parameters ( int p0, double k, vec ref0, double l0 ) {
	1677	p = p0;
	1678	double delta = 2 / ( k * k ) + p + 3;
	1679	sqd = sqrt ( delta - p - 1 );
	1680	l = l0;
	1681	refl = pow ( ref0, 1 - l );
	1682	iepdf.set_parameters ( eye ( p ), delta );
	1683	};
	1684
	1685	void validate() {
	1686	pdf_internal<eiWishartCh>::validate();
	1687	dimc = iepdf.dimension();
	1688	}
	1689
	1690	void condition ( const vec &c ) {
	1691	vec z = c;
	1692	int ri = 0;
	1693	for ( int i = 0; i < p*p; i += ( p + 1 ) ) {//trace diagonal element
	1694	z ( i ) = pow ( z ( i ), l ) * refl ( ri );
	1695	ri++;
	1696	}
	1697
	1698	iepdf._setY ( sqd*z );
	1699	}
1688	1700	};
1689	1701
…	…
1691	1703	enum RESAMPLING_METHOD { MULTINOMIAL = 0, STRATIFIED = 1, SYSTEMATIC = 3 };
1692	1704
1693		//! Shortcut for multinomial.sample(int n). Various simplifications are allowed see RESAMPLING_METHOD
	1705	//! Shortcut for multinomial.sample(int n). Various simplifications are allowed see RESAMPLING_METHOD
1694	1706	void resample(const vec &w, ivec &ind, RESAMPLING_METHOD=SYSTEMATIC);
1695	1707
…	…
1700	1712	class eEmp: public epdf {
1701	1713	protected :
1702		//! Number of particles
1703		int n;
1704		//! Sample weights \f$w\f$
1705		vec w;
1706		//! Samples \f$x^{(i)}, i=1..n\f$
1707		Array<vec> samples;
1708		public:
1709		//! \name Constructors
1710		//!@{
1711		eEmp () : epdf (), w (), samples () {};
1712		//! copy constructor
1713		eEmp ( const eEmp &e ) : epdf ( e ), w ( e.w ), samples ( e.samples ) {};
1714		//!@}
1715
1716		//! Set samples and weights
1717		void set_statistics ( const vec &w0, const epdf &pdf0 );
1718		//! Set samples and weights
1719		void set_statistics ( const epdf &pdf0 , int n ) {
1720		set_statistics ( ones ( n ) / n, pdf0 );
1721		};
1722		//! Set sample
1723		void set_samples ( const epdf* pdf0 );
1724		//! Set sample
1725		void set_parameters ( int n0, bool copy = true ) {
1726		n = n0;
1727		w.set_size ( n0, copy );
1728		samples.set_size ( n0, copy );
1729		};
1730		//! Set samples
1731		void set_parameters ( const Array<vec> &Av ) {
1732		n = Av.size();
1733		w = 1 / n * ones ( n );
1734		samples = Av;
1735		};
1736		virtual void validate ();
1737		//! Potentially dangerous, use with care.
1738		vec& _w() {
1739		return w;
1740		};
1741		//! Potentially dangerous, use with care.
1742		const vec& _w() const {
1743		return w;
1744		};
1745		//! access function
1746		Array<vec>& _samples() {
1747		return samples;
1748		};
1749		//! access function
1750		const vec& _sample ( int i ) const {
1751		return samples ( i );
1752		};
1753		//! access function
1754		const Array<vec>& _samples() const {
1755		return samples;
1756		};
1757		//! Function performs resampling, i.e. removal of low-weight samples and duplication of high-weight samples such that the new samples represent the same density.
1758		void resample ( RESAMPLING_METHOD method = SYSTEMATIC );
1759
1760		//! inherited operation : NOT implemented
1761		vec sample() const NOT_IMPLEMENTED(0);
1762
1763		//! inherited operation : NOT implemented
1764		double evallog ( const vec &val ) const NOT_IMPLEMENTED(0);
1765
1766		vec mean() const {
1767		vec pom = zeros ( dim );
1768		for ( int i = 0; i < n; i++ ) {
1769		pom += samples ( i ) * w ( i );
1770		}
1771		return pom;
1772		}
1773		vec variance() const {
1774		vec pom = zeros ( dim );
1775		for ( int i = 0; i < n; i++ ) {
1776		pom += pow ( samples ( i ), 2 ) * w ( i );
1777		}
1778		return pom - pow ( mean(), 2 );
1779		}
1780		//! For this class, qbounds are minimum and maximum value of the population!
1781		void qbounds ( vec &lb, vec &ub, double perc = 0.95 ) const;
1782
1783		void to_setting ( Setting &set ) const;
1784
1785		/*! Create object from the following structure
1786
1787		\code
1788		class = 'eEmp';
1789		samples = [...]; % array of samples
1790		w = [...]; % weights of samples stored in vector
1791		--- inherited fields ---
1792		bdm::epdf::from_setting
1793		\endcode
1794		*/
1795		void from_setting ( const Setting &set );
	1714	//! Number of particles
	1715	int n;
	1716	//! Sample weights \f$w\f$
	1717	vec w;
	1718	//! Samples \f$x^{(i)}, i=1..n\f$
	1719	Array<vec> samples;
	1720	public:
	1721	//! \name Constructors
	1722	//!@{
	1723	eEmp () : epdf (), w (), samples () {};
	1724	//! copy constructor
	1725	eEmp ( const eEmp &e ) : epdf ( e ), w ( e.w ), samples ( e.samples ) {};
	1726	//!@}
	1727
	1728	//! Set samples and weights
	1729	void set_statistics ( const vec &w0, const epdf &pdf0 );
	1730	//! Set samples and weights
	1731	void set_statistics ( const epdf &pdf0 , int n ) {
	1732	set_statistics ( ones ( n ) / n, pdf0 );
	1733	};
	1734	//! Set sample
	1735	void set_samples ( const epdf* pdf0 );
	1736	//! Set sample
	1737	void set_parameters ( int n0, bool copy = true ) {
	1738	n = n0;
	1739	w.set_size ( n0, copy );
	1740	samples.set_size ( n0, copy );
	1741	};
	1742	//! Set samples
	1743	void set_parameters ( const Array<vec> &Av ) {
	1744	n = Av.size();
	1745	w = 1 / n * ones ( n );
	1746	samples = Av;
	1747	};
	1748	virtual void validate ();
	1749	//! Potentially dangerous, use with care.
	1750	vec& _w() {
	1751	return w;
	1752	};
	1753	//! Potentially dangerous, use with care.
	1754	const vec& _w() const {
	1755	return w;
	1756	};
	1757	//! access function
	1758	Array<vec>& _samples() {
	1759	return samples;
	1760	};
	1761	//! access function
	1762	const vec& _sample ( int i ) const {
	1763	return samples ( i );
	1764	};
	1765	//! access function
	1766	const Array<vec>& _samples() const {
	1767	return samples;
	1768	};
	1769	//! Function performs resampling, i.e. removal of low-weight samples and duplication of high-weight samples such that the new samples represent the same density.
	1770	void resample ( RESAMPLING_METHOD method = SYSTEMATIC );
	1771
	1772	//! inherited operation : NOT implemented
	1773	vec sample() const NOT_IMPLEMENTED(0);
	1774
	1775	//! inherited operation : NOT implemented
	1776	double evallog ( const vec &val ) const NOT_IMPLEMENTED(0);
	1777
	1778	vec mean() const {
	1779	vec pom = zeros ( dim );
	1780	for ( int i = 0; i < n; i++ ) {
	1781	pom += samples ( i ) * w ( i );
	1782	}
	1783	return pom;
	1784	}
	1785	vec variance() const {
	1786	vec pom = zeros ( dim );
	1787	for ( int i = 0; i < n; i++ ) {
	1788	pom += pow ( samples ( i ), 2 ) * w ( i );
	1789	}
	1790	return pom - pow ( mean(), 2 );
	1791	}
	1792	//! For this class, qbounds are minimum and maximum value of the population!
	1793	void qbounds ( vec &lb, vec &ub, double perc = 0.95 ) const;
	1794
	1795	void to_setting ( Setting &set ) const;
	1796
	1797	/*! Create object from the following structure
	1798
	1799	\code
	1800	class = 'eEmp';
	1801	samples = [...]; % array of samples
	1802	w = [...]; % weights of samples stored in vector
	1803	--- inherited fields ---
	1804	bdm::epdf::from_setting
	1805	\endcode
	1806	*/
	1807	void from_setting ( const Setting &set );
1796	1808	};
1797	1809	UIREGISTER(eEmp);
…	…
1803	1815	void enorm<sq_T>::set_parameters ( const vec &mu0, const sq_T &R0 ) {
1804	1816	//Fixme test dimensions of mu0 and R0;
1805		mu = mu0;
1806		R = R0;
1807		validate();
	1817	mu = mu0;
	1818	R = R0;
	1819	validate();
1808	1820	};
1809	1821
1810	1822	template<class sq_T>
1811	1823	void enorm<sq_T>::from_setting ( const Setting &set ) {
1812		epdf::from_setting ( set ); //reads rv
1813
1814		UI::get ( mu, set, "mu", UI::compulsory );
1815		mat Rtmp;// necessary for conversion
1816		UI::get ( Rtmp, set, "R", UI::compulsory );
1817		R = Rtmp; // conversion
	1824	epdf::from_setting ( set ); //reads rv
	1825
	1826	UI::get ( mu, set, "mu", UI::compulsory );
	1827	mat Rtmp;// necessary for conversion
	1828	UI::get ( Rtmp, set, "R", UI::compulsory );
	1829	R = Rtmp; // conversion
1818	1830	}
1819	1831
1820	1832	template<class sq_T>
1821	1833	void enorm<sq_T>::validate() {
1822		eEF::validate();
1823		bdm_assert ( mu.length() == R.rows(), "mu and R parameters do not match" );
1824		dim = mu.length();
1825		}
	1834	eEF::validate();
	1835	bdm_assert ( mu.length() == R.rows(), "mu and R parameters do not match" );
	1836	dim = mu.length();
	1837	}
1826	1838
1827	1839	template<class sq_T>
1828	1840	void enorm<sq_T>::to_setting ( Setting &set ) const {
1829		epdf::to_setting ( set ); //reads rv
1830		UI::save ( mu, set, "mu");
1831		UI::save ( R.to_mat(), set, "R");
	1841	epdf::to_setting ( set ); //reads rv
	1842	UI::save ( mu, set, "mu");
	1843	UI::save ( R.to_mat(), set, "R");
1832	1844	}
1833	1845
…	…
1836	1848	template<class sq_T>
1837	1849	void enorm<sq_T>::dupdate ( mat &v, double nu ) {
1838		//
	1850	//
1839	1851	};
1840	1852
…	…
1846	1858	template<class sq_T>
1847	1859	vec enorm<sq_T>::sample() const {
1848		vec x ( dim );
	1860	vec x ( dim );
1849	1861	#pragma omp critical
1850		NorRNG.sample_vector ( dim, x );
1851		vec smp = R.sqrt_mult ( x );
1852
1853		smp += mu;
1854		return smp;
	1862	NorRNG.sample_vector ( dim, x );
	1863	vec smp = R.sqrt_mult ( x );
	1864
	1865	smp += mu;
	1866	return smp;
1855	1867	};
1856	1868
…	…
1865	1877	template<class sq_T>
1866	1878	double enorm<sq_T>::evallog_nn ( const vec &val ) const {
1867		// 1.83787706640935 = log(2pi)
1868		double tmp = -0.5 * ( R.invqform ( mu - val ) );// - lognc();
1869		return tmp;
	1879	// 1.83787706640935 = log(2pi)
	1880	double tmp = -0.5 * ( R.invqform ( mu - val ) );// - lognc();
	1881	return tmp;
1870	1882	};
1871	1883
1872	1884	template<class sq_T>
1873	1885	inline double enorm<sq_T>::lognc () const {
1874		// 1.83787706640935 = log(2pi)
1875		double tmp = 0.5 * ( R.cols() * 1.83787706640935 + R.logdet() );
1876		return tmp;
	1886	// 1.83787706640935 = log(2pi)
	1887	double tmp = 0.5 * ( R.cols() * 1.83787706640935 + R.logdet() );
	1888	return tmp;
1877	1889	};
1878	1890
…	…
1906	1918	template<class sq_T>
1907	1919	shared_ptr<epdf> enorm<sq_T>::marginal ( const RV &rvn ) const {
1908		enorm<sq_T> *tmp = new enorm<sq_T> ();
1909		shared_ptr<epdf> narrow ( tmp );
1910		marginal ( rvn, *tmp );
1911		return narrow;
	1920	enorm<sq_T> *tmp = new enorm<sq_T> ();
	1921	shared_ptr<epdf> narrow ( tmp );
	1922	marginal ( rvn, *tmp );
	1923	return narrow;
1912	1924	}
1913	1925
1914	1926	template<class sq_T>
1915	1927	void enorm<sq_T>::marginal ( const RV &rvn, enorm<sq_T> &target ) const {
1916		bdm_assert ( isnamed(), "rv description is not assigned" );
1917		ivec irvn = rvn.dataind ( rv );
1918
1919		sq_T Rn ( R, irvn ); // select rows and columns of R
1920
1921		target.set_rv ( rvn );
1922		target.set_parameters ( mu ( irvn ), Rn );
	1928	bdm_assert ( isnamed(), "rv description is not assigned" );
	1929	ivec irvn = rvn.dataind ( rv );
	1930
	1931	sq_T Rn ( R, irvn ); // select rows and columns of R
	1932
	1933	target.set_rv ( rvn );
	1934	target.set_parameters ( mu ( irvn ), Rn );
1923	1935	}
1924	1936
1925	1937	template<class sq_T>
1926	1938	shared_ptr<pdf> enorm<sq_T>::condition ( const RV &rvn ) const {
1927		mlnorm<sq_T> *tmp = new mlnorm<sq_T> ();
1928		shared_ptr<pdf> narrow ( tmp );
1929		condition ( rvn, *tmp );
1930		return narrow;
	1939	mlnorm<sq_T> *tmp = new mlnorm<sq_T> ();
	1940	shared_ptr<pdf> narrow ( tmp );
	1941	condition ( rvn, *tmp );
	1942	return narrow;
1931	1943	}
1932	1944
1933	1945	template<class sq_T>
1934	1946	void enorm<sq_T>::condition ( const RV &rvn, pdf &target ) const {
1935		typedef mlnorm<sq_T> TMlnorm;
1936
1937		bdm_assert ( isnamed(), "rvs are not assigned" );
1938		TMlnorm &uptarget = dynamic_cast<TMlnorm &> ( target );
1939
1940		RV rvc = rv.subt ( rvn );
1941		bdm_assert ( ( rvc._dsize() + rvn._dsize() == rv._dsize() ), "wrong rvn" );
1942		//Permutation vector of the new R
1943		ivec irvn = rvn.dataind ( rv );
1944		ivec irvc = rvc.dataind ( rv );
1945		ivec perm = concat ( irvn , irvc );
1946		sq_T Rn ( R, perm );
1947
1948		//fixme - could this be done in general for all sq_T?
1949		mat S = Rn.to_mat();
1950		//fixme
1951		int n = rvn._dsize() - 1;
1952		int end = R.rows() - 1;
1953		mat S11 = S.get ( 0, n, 0, n );
1954		mat S12 = S.get ( 0, n , rvn._dsize(), end );
1955		mat S22 = S.get ( rvn._dsize(), end, rvn._dsize(), end );
1956
1957		vec mu1 = mu ( irvn );
1958		vec mu2 = mu ( irvc );
1959		mat A = S12 * inv ( S22 );
1960		sq_T R_n ( S11 - A *S12.T() );
1961
1962		uptarget.set_rv ( rvn );
1963		uptarget.set_rvc ( rvc );
1964		uptarget.set_parameters ( A, mu1 - A*mu2, R_n );
1965		uptarget.validate();
	1947	typedef mlnorm<sq_T> TMlnorm;
	1948
	1949	bdm_assert ( isnamed(), "rvs are not assigned" );
	1950	TMlnorm &uptarget = dynamic_cast<TMlnorm &> ( target );
	1951
	1952	RV rvc = rv.subt ( rvn );
	1953	bdm_assert ( ( rvc._dsize() + rvn._dsize() == rv._dsize() ), "wrong rvn" );
	1954	//Permutation vector of the new R
	1955	ivec irvn = rvn.dataind ( rv );
	1956	ivec irvc = rvc.dataind ( rv );
	1957	ivec perm = concat ( irvn , irvc );
	1958	sq_T Rn ( R, perm );
	1959
	1960	//fixme - could this be done in general for all sq_T?
	1961	mat S = Rn.to_mat();
	1962	//fixme
	1963	int n = rvn._dsize() - 1;
	1964	int end = R.rows() - 1;
	1965	mat S11 = S.get ( 0, n, 0, n );
	1966	mat S12 = S.get ( 0, n , rvn._dsize(), end );
	1967	mat S22 = S.get ( rvn._dsize(), end, rvn._dsize(), end );
	1968
	1969	vec mu1 = mu ( irvn );
	1970	vec mu2 = mu ( irvc );
	1971	mat A = S12 * inv ( S22 );
	1972	sq_T R_n ( S11 - A *S12.T() );
	1973
	1974	uptarget.set_rv ( rvn );
	1975	uptarget.set_rvc ( rvc );
	1976	uptarget.set_parameters ( A, mu1 - A*mu2, R_n );
	1977	uptarget.validate();
1966	1978	}
1967	1979
1968	1980	/! \brief Dirac delta function distribution /
1969		class dirac: public epdf{
1970		public:
1971		vec point;
1972		public:
1973		double evallog (const vec &dt) const {return -inf;}
1974		vec mean () const {return point;}
1975		vec variance () const {return zeros(point.length());}
1976		void qbounds ( vec &lb, vec &ub, double percentage = 0.95 ) const { lb = point; ub = point;}
1977		//! access
1978		const vec& _point() {return point;}
1979		void set_point(const vec& p){point =p; dim=p.length();}
1980		vec sample() const {return point;}
1981		};
	1981	class dirac: public epdf {
	1982	public:
	1983	vec point;
	1984	public:
	1985	double evallog (const vec &dt) const {
	1986	return -inf;
	1987	}
	1988	vec mean () const {
	1989	return point;
	1990	}
	1991	vec variance () const {
	1992	return zeros(point.length());
	1993	}
	1994	void qbounds ( vec &lb, vec &ub, double percentage = 0.95 ) const {
	1995	lb = point;
	1996	ub = point;
	1997	}
	1998	//! access
	1999	const vec& _point() {
	2000	return point;
	2001	}
	2002	void set_point(const vec& p) {
	2003	point =p;
	2004	dim=p.length();
	2005	}
	2006	vec sample() const {
	2007	return point;
	2008	}
	2009	};
1982	2010
1983	2011
…	…
1986	2014	template<class sq_T>
1987	2015	void mgnorm<sq_T >::set_parameters ( const shared_ptr<fnc> &g0, const sq_T &R0 ) {
1988		g = g0;
1989		this->iepdf.set_parameters ( zeros ( g->dimension() ), R0 );
	2016	g = g0;
	2017	this->iepdf.set_parameters ( zeros ( g->dimension() ), R0 );
1990	2018	}
1991	2019
1992	2020	template<class sq_T>
1993	2021	void mgnorm<sq_T >::condition ( const vec &cond ) {
1994		this->iepdf._mu() = g->eval ( cond );
	2022	this->iepdf._mu() = g->eval ( cond );
1995	2023	};
1996	2024
…	…
1998	2026	template<class sq_T>
1999	2027	std::ostream &operator<< ( std::ostream &os, mlnorm<sq_T> &ml ) {
2000		os << "A:" << ml.A << endl;
2001		os << "mu:" << ml.mu_const << endl;
2002		os << "R:" << ml._R() << endl;
2003		return os;
	2028	os << "A:" << ml.A << endl;
	2029	os << "mu:" << ml.mu_const << endl;
	2030	os << "R:" << ml._R() << endl;
	2031	return os;
2004	2032	};
2005	2033

Context Navigation

Changeset 1064 for library/bdm/stat/exp_family.h

Legend:

library/bdm/stat/exp_family.h

Download in other formats: