Context Navigation

← Previous Change
Next Change →

exp_family.h

Timestamp:

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

Author:

mido

Message:

ASTYLER RUN OVER THE WHOLE LIBRARY, JUPEE

Files:

: 1 modified

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

Legend:

: Unmodified
: Added
: Removed

library/bdm/stat/exp_family.h

r735	r737
19	19	#include "../math/chmat.h"
20	20
21		namespace bdm
22		{
	21	namespace bdm {
23	22
24	23
…	…
36	35	*/
37	36
38		class eEF : public epdf
39		{
40		public:
	37	class eEF : public epdf {
	38	public:
41	39	// eEF() :epdf() {};
42		//! default constructor
43		eEF () : epdf () {};
44		//! logarithm of the normalizing constant, \f$\mathcal{I}\f$
45		virtual double lognc() const = 0;
46
47		//!Evaluate normalized log-probability
48		virtual double evallog_nn (const vec &val) const {
49		bdm_error ("Not implemented");
50		return 0.0;
51		}
52
53		//!Evaluate normalized log-probability
54		virtual double evallog (const vec &val) const {
55		double tmp;
56		tmp = evallog_nn (val) - lognc();
57		return tmp;
58		}
59		//!Evaluate normalized log-probability for many samples
60		virtual vec evallog_mat (const mat &Val) const {
61		vec x (Val.cols());
62		for (int i = 0;i < Val.cols();i++) {x (i) = evallog_nn (Val.get_col (i)) ;}
63		return x -lognc();
64		}
65		//!Evaluate normalized log-probability for many samples
66		virtual vec evallog_mat (const Array<vec> &Val) const {
67		vec x (Val.length());
68		for (int i = 0;i < Val.length();i++) {x (i) = evallog_nn (Val (i)) ;}
69		return x -lognc();
70		}
71
72		//!Power of the density, used e.g. to flatten the density
73		virtual void pow (double p) {
74		bdm_error ("Not implemented");
75		}
	40	//! default constructor
	41	eEF () : epdf () {};
	42	//! logarithm of the normalizing constant, \f$\mathcal{I}\f$
	43	virtual double lognc() const = 0;
	44
	45	//!Evaluate normalized log-probability
	46	virtual double evallog_nn ( const vec &val ) const {
	47	bdm_error ( "Not implemented" );
	48	return 0.0;
	49	}
	50
	51	//!Evaluate normalized log-probability
	52	virtual double evallog ( const vec &val ) const {
	53	double tmp;
	54	tmp = evallog_nn ( val ) - lognc();
	55	return tmp;
	56	}
	57	//!Evaluate normalized log-probability for many samples
	58	virtual vec evallog_mat ( const mat &Val ) const {
	59	vec x ( Val.cols() );
	60	for ( int i = 0; i < Val.cols(); i++ ) {
	61	x ( i ) = evallog_nn ( Val.get_col ( i ) ) ;
	62	}
	63	return x - lognc();
	64	}
	65	//!Evaluate normalized log-probability for many samples
	66	virtual vec evallog_mat ( const Array<vec> &Val ) const {
	67	vec x ( Val.length() );
	68	for ( int i = 0; i < Val.length(); i++ ) {
	69	x ( i ) = evallog_nn ( Val ( i ) ) ;
	70	}
	71	return x - lognc();
	72	}
	73
	74	//!Power of the density, used e.g. to flatten the density
	75	virtual void pow ( double p ) {
	76	bdm_error ( "Not implemented" );
	77	}
76	78	};
77	79
78	80
79	81	//! Estimator for Exponential family
80		class BMEF : public BM
81		{
82		protected:
83		//! forgetting factor
84		double frg;
85		//! cached value of lognc() in the previous step (used in evaluation of \c ll )
86		double last_lognc;
87		public:
88		//! Default constructor (=empty constructor)
89		BMEF (double frg0 = 1.0) : BM (), frg (frg0) {}
90		//! Copy constructor
91		BMEF (const BMEF &B) : BM (B), frg (B.frg), last_lognc (B.last_lognc) {}
92		//!get statistics from another model
93		virtual void set_statistics (const BMEF* BM0) {
94		bdm_error ("Not implemented");
95		}
96
97		//! Weighted update of sufficient statistics (Bayes rule)
98		virtual void bayes_weighted (const vec &data, const vec &cond=empty_vec, const double w=1.0) {};
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) {
104		bdm_error ("Not implemented");
105		}
106
107		BMEF* _copy_ () const {
108		bdm_error ("function _copy_ not implemented for this BM");
109		return NULL;
110		}
	82	class BMEF : public BM {
	83	protected:
	84	//! forgetting factor
	85	double frg;
	86	//! cached value of lognc() in the previous step (used in evaluation of \c ll )
	87	double last_lognc;
	88	public:
	89	//! Default constructor (=empty constructor)
	90	BMEF ( double frg0 = 1.0 ) : BM (), frg ( frg0 ) {}
	91	//! Copy constructor
	92	BMEF ( const BMEF &B ) : BM ( B ), frg ( B.frg ), last_lognc ( B.last_lognc ) {}
	93	//!get statistics from another model
	94	virtual void set_statistics ( const BMEF* BM0 ) {
	95	bdm_error ( "Not implemented" );
	96	}
	97
	98	//! Weighted update of sufficient statistics (Bayes rule)
	99	virtual void bayes_weighted ( const vec &data, const vec &cond = empty_vec, const double w = 1.0 ) {};
	100	//original Bayes
	101	void bayes ( const vec &yt, const vec &cond = empty_vec );
	102
	103	//!Flatten the posterior according to the given BMEF (of the same type!)
	104	virtual void flatten ( const BMEF * B ) {
	105	bdm_error ( "Not implemented" );
	106	}
	107
	108	BMEF* _copy_ () const {
	109	bdm_error ( "function _copy_ not implemented for this BM" );
	110	return NULL;
	111	}
111	112	};
112	113
…	…
120	121	*/
121	122	template<class sq_T>
122		class enorm : public eEF
123		{
124		protected:
125		//! mean value
126		vec mu;
127		//! Covariance matrix in decomposed form
128		sq_T R;
129		public:
130		//!\name Constructors
131		//!@{
132
133		enorm () : eEF (), mu (), R () {};
134		enorm (const vec &mu, const sq_T &R) {set_parameters (mu, R);}
135		void set_parameters (const vec &mu, const sq_T &R);
136		/*! Create Normal density
137		\f[ f(rv) = N(\mu, R) \f]
138		from structure
139		\code
140		class = 'enorm<ldmat>', (OR) 'enorm<chmat>', (OR) 'enorm<fsqmat>';
141		mu = []; // mean value
142		R = []; // variance, square matrix of appropriate dimension
143		\endcode
144		*/
145		void from_setting (const Setting &root);
146		void validate() {
147		bdm_assert (mu.length() == R.rows(), "mu and R parameters do not match");
148		dim = mu.length();
149		}
150		//!@}
151
152		//! \name Mathematical operations
153		//!@{
154
155		//! dupdate in exponential form (not really handy)
156		void dupdate (mat &v, double nu = 1.0);
157
158		vec sample() const;
159
160		double evallog_nn (const vec &val) const;
161		double lognc () const;
162		vec mean() const {return mu;}
163		vec variance() const {return diag (R.to_mat());}
	123	class enorm : public eEF {
	124	protected:
	125	//! mean value
	126	vec mu;
	127	//! Covariance matrix in decomposed form
	128	sq_T R;
	129	public:
	130	//!\name Constructors
	131	//!@{
	132
	133	enorm () : eEF (), mu (), R () {};
	134	enorm ( const vec &mu, const sq_T &R ) {
	135	set_parameters ( mu, R );
	136	}
	137	void set_parameters ( const vec &mu, const sq_T &R );
	138	/*! Create Normal density
	139	\f[ f(rv) = N(\mu, R) \f]
	140	from structure
	141	\code
	142	class = 'enorm<ldmat>', (OR) 'enorm<chmat>', (OR) 'enorm<fsqmat>';
	143	mu = []; // mean value
	144	R = []; // variance, square matrix of appropriate dimension
	145	\endcode
	146	*/
	147	void from_setting ( const Setting &root );
	148	void validate() {
	149	bdm_assert ( mu.length() == R.rows(), "mu and R parameters do not match" );
	150	dim = mu.length();
	151	}
	152	//!@}
	153
	154	//! \name Mathematical operations
	155	//!@{
	156
	157	//! dupdate in exponential form (not really handy)
	158	void dupdate ( mat &v, double nu = 1.0 );
	159
	160	vec sample() const;
	161
	162	double evallog_nn ( const vec &val ) const;
	163	double lognc () const;
	164	vec mean() const {
	165	return mu;
	166	}
	167	vec variance() const {
	168	return diag ( R.to_mat() );
	169	}
164	170	// mlnorm<sq_T>* condition ( const RV &rvn ) const ; <=========== fails to cmpile. Why?
165		shared_ptr<pdf> condition ( const RV &rvn ) const;
166
167		// target not typed to mlnorm<sq_T, enorm<sq_T> > &
168		// because that doesn't compile (perhaps because we
169		// haven't finished defining enorm yet), but the type
170		// is required
171		void condition ( const RV &rvn, pdf &target ) const;
172
173		shared_ptr<epdf> marginal (const RV &rvn ) const;
174		void marginal ( const RV &rvn, enorm<sq_T> &target ) const;
175		//!@}
176
177		//! \name Access to attributes
178		//!@{
179
180		vec& _mu() {return mu;}
181		const vec& _mu() const {return mu;}
182		void set_mu (const vec mu0) { mu = mu0;}
183		sq_T& _R() {return R;}
184		const sq_T& _R() const {return R;}
185		//!@}
186
187		};
188		UIREGISTER2 (enorm, chmat);
	171	shared_ptr<pdf> condition ( const RV &rvn ) const;
	172
	173	// target not typed to mlnorm<sq_T, enorm<sq_T> > &
	174	// because that doesn't compile (perhaps because we
	175	// haven't finished defining enorm yet), but the type
	176	// is required
	177	void condition ( const RV &rvn, pdf &target ) const;
	178
	179	shared_ptr<epdf> marginal ( const RV &rvn ) const;
	180	void marginal ( const RV &rvn, enorm<sq_T> &target ) const;
	181	//!@}
	182
	183	//! \name Access to attributes
	184	//!@{
	185
	186	vec& _mu() {
	187	return mu;
	188	}
	189	const vec& _mu() const {
	190	return mu;
	191	}
	192	void set_mu ( const vec mu0 ) {
	193	mu = mu0;
	194	}
	195	sq_T& _R() {
	196	return R;
	197	}
	198	const sq_T& _R() const {
	199	return R;
	200	}
	201	//!@}
	202
	203	};
	204	UIREGISTER2 ( enorm, chmat );
189	205	SHAREDPTR2 ( enorm, chmat );
190		UIREGISTER2 (~~enorm, ldmat~~);
	206	UIREGISTER2 ( enorm, ldmat );
191	207	SHAREDPTR2 ( enorm, ldmat );
192		UIREGISTER2 (~~enorm, fsqmat~~);
	208	UIREGISTER2 ( enorm, fsqmat );
193	209	SHAREDPTR2 ( enorm, fsqmat );
194	210
…	…
200	216	*
201	217	*/
202		class egiw : public eEF
203		{
204		protected:
205		//! Extended information matrix of sufficient statistics
206		ldmat V;
207		//! Number of data records (degrees of freedom) of sufficient statistics
208		double nu;
209		//! Dimension of the output
210		int dimx;
211		//! Dimension of the regressor
212		int nPsi;
213		public:
214		//!\name Constructors
215		//!@{
216		egiw() : eEF() {};
217		egiw (int dimx0, ldmat V0, double nu0 = -1.0) : eEF() {set_parameters (dimx0, V0, nu0);};
218
219		void set_parameters (int dimx0, ldmat V0, double nu0 = -1.0);
220		//!@}
221
222		vec sample() const;
223		mat sample_mat(int n) const;
224		vec mean() const;
225		vec variance() const;
226		void sample_mat(mat &Mi, chmat &Ri)const;
227
228		void factorize(mat &M, ldmat &Vz, ldmat &Lam) const;
229		//! LS estimate of \f$\theta\f$
230		vec est_theta() const;
231
232		//! Covariance of the LS estimate
233		ldmat est_theta_cov() const;
234
235		//! expected values of the linear coefficient and the covariance matrix are written to \c M and \c R , respectively
236		void mean_mat (mat &M, mat&R) const;
237		//! In this instance, val= [theta, r]. For multivariate instances, it is stored columnwise val = [theta_1 theta_2 ... r_1 r_2 ]
238		double evallog_nn (const vec &val) const;
239		double lognc () const;
240		void pow (double p) {V = p;nu = p;};
241
242		//! \name Access attributes
243		//!@{
244
245		ldmat& _V() {return V;}
246		const ldmat& _V() const {return V;}
247		double& _nu() {return nu;}
248		const double& _nu() const {return nu;}
249		const int & _dimx() const {return dimx;}
250
251		/*! Create Gauss-inverse-Wishart density
252		\f[ f(rv) = GiW(V,\nu) \f]
253		from structure
254		\code
255		class = 'egiw';
256		V = []; // square matrix
257		dV = []; // vector of diagonal of V (when V not given)
258		nu = []; // scalar \nu ((almost) degrees of freedom)
259		// when missing, it will be computed to obtain proper pdf
260		dimx = []; // dimension of the wishart part
261		rv = RV({'name'}) // description of RV
262		rvc = RV({'name'}) // description of RV in condition
263		\endcode
264		*/
265
266		void from_setting (const Setting &set) {
267		epdf::from_setting(set);
268		UI::get (dimx, set, "dimx", UI::compulsory);
269		if (!UI::get (nu, set, "nu", UI::optional)) {nu=-1;}
270		mat V;
271		if (!UI::get (V, set, "V", UI::optional)){
272		vec dV;
273		UI::get (dV, set, "dV", UI::compulsory);
274		set_parameters (dimx, ldmat(dV), nu);
275
276		} else {
277		set_parameters (dimx, V, nu);
278		}
279		}
280
281		void to_setting ( Setting& set ) const{
282		epdf::to_setting(set);
283		string s("egiw");
284		UI::save(s,set,"class");
285		UI::save(dimx,set,"dimx");
286		UI::save(V.to_mat(),set,"V");
287		UI::save(nu,set,"nu");
288		};
289
290		void validate(){
291		// check sizes, rvs etc.
292		}
293		void log_register( bdm::logger& L, const string& prefix ){
294		if (log_level==3){
295		root::log_register(L,prefix);
296		logrec->ids.set_length(2);
297		int th_dim=dimension()-dimx*(dimx+1)/2;
298		logrec->ids(0)=L.add_vector(RV("",th_dim), prefix + logrec->L.prefix_sep() +"mean");
299		logrec->ids(1)=L.add_vector(RV("",th_dim*th_dim),prefix + logrec->L.prefix_sep() + "variance");
300		} else {
301		epdf::log_register(L,prefix);
302		}
303		}
304		void log_write() const {
305		if (log_level==3){
306		mat M;
307		ldmat Lam;
308		ldmat Vz;
309		factorize(M,Vz,Lam);
310		logrec->L.log_vector(logrec->ids(0), est_theta() );
311		logrec->L.log_vector(logrec->ids(1), cvectorize(est_theta_cov().to_mat()));
312		} else {
313		epdf::log_write();
314		}
315
316		}
317		//!@}
	218	class egiw : public eEF {
	219	protected:
	220	//! Extended information matrix of sufficient statistics
	221	ldmat V;
	222	//! Number of data records (degrees of freedom) of sufficient statistics
	223	double nu;
	224	//! Dimension of the output
	225	int dimx;
	226	//! Dimension of the regressor
	227	int nPsi;
	228	public:
	229	//!\name Constructors
	230	//!@{
	231	egiw() : eEF() {};
	232	egiw ( int dimx0, ldmat V0, double nu0 = -1.0 ) : eEF() {
	233	set_parameters ( dimx0, V0, nu0 );
	234	};
	235
	236	void set_parameters ( int dimx0, ldmat V0, double nu0 = -1.0 );
	237	//!@}
	238
	239	vec sample() const;
	240	mat sample_mat ( int n ) const;
	241	vec mean() const;
	242	vec variance() const;
	243	void sample_mat ( mat &Mi, chmat &Ri ) const;
	244
	245	void factorize ( mat &M, ldmat &Vz, ldmat &Lam ) const;
	246	//! LS estimate of \f$\theta\f$
	247	vec est_theta() const;
	248
	249	//! Covariance of the LS estimate
	250	ldmat est_theta_cov() const;
	251
	252	//! expected values of the linear coefficient and the covariance matrix are written to \c M and \c R , respectively
	253	void mean_mat ( mat &M, mat&R ) const;
	254	//! In this instance, val= [theta, r]. For multivariate instances, it is stored columnwise val = [theta_1 theta_2 ... r_1 r_2 ]
	255	double evallog_nn ( const vec &val ) const;
	256	double lognc () const;
	257	void pow ( double p ) {
	258	V *= p;
	259	nu *= p;
	260	};
	261
	262	//! \name Access attributes
	263	//!@{
	264
	265	ldmat& _V() {
	266	return V;
	267	}
	268	const ldmat& _V() const {
	269	return V;
	270	}
	271	double& _nu() {
	272	return nu;
	273	}
	274	const double& _nu() const {
	275	return nu;
	276	}
	277	const int & _dimx() const {
	278	return dimx;
	279	}
	280
	281	/*! Create Gauss-inverse-Wishart density
	282	\f[ f(rv) = GiW(V,\nu) \f]
	283	from structure
	284	\code
	285	class = 'egiw';
	286	V = []; // square matrix
	287	dV = []; // vector of diagonal of V (when V not given)
	288	nu = []; // scalar \nu ((almost) degrees of freedom)
	289	// when missing, it will be computed to obtain proper pdf
	290	dimx = []; // dimension of the wishart part
	291	rv = RV({'name'}) // description of RV
	292	rvc = RV({'name'}) // description of RV in condition
	293	\endcode
	294	*/
	295
	296	void from_setting ( const Setting &set ) {
	297	epdf::from_setting ( set );
	298	UI::get ( dimx, set, "dimx", UI::compulsory );
	299	if ( !UI::get ( nu, set, "nu", UI::optional ) ) {
	300	nu = -1;
	301	}
	302	mat V;
	303	if ( !UI::get ( V, set, "V", UI::optional ) ) {
	304	vec dV;
	305	UI::get ( dV, set, "dV", UI::compulsory );
	306	set_parameters ( dimx, ldmat ( dV ), nu );
	307
	308	} else {
	309	set_parameters ( dimx, V, nu );
	310	}
	311	}
	312
	313	void to_setting ( Setting& set ) const {
	314	epdf::to_setting ( set );
	315	string s ( "egiw" );
	316	UI::save ( s, set, "class" );
	317	UI::save ( dimx, set, "dimx" );
	318	UI::save ( V.to_mat(), set, "V" );
	319	UI::save ( nu, set, "nu" );
	320	};
	321
	322	void validate() {
	323	// check sizes, rvs etc.
	324	}
	325	void log_register ( bdm::logger& L, const string& prefix ) {
	326	if ( log_level == 3 ) {
	327	root::log_register ( L, prefix );
	328	logrec->ids.set_length ( 2 );
	329	int th_dim = dimension() - dimx * ( dimx + 1 ) / 2;
	330	logrec->ids ( 0 ) = L.add_vector ( RV ( "", th_dim ), prefix + logrec->L.prefix_sep() + "mean" );
	331	logrec->ids ( 1 ) = L.add_vector ( RV ( "", th_dim * th_dim ), prefix + logrec->L.prefix_sep() + "variance" );
	332	} else {
	333	epdf::log_register ( L, prefix );
	334	}
	335	}
	336	void log_write() const {
	337	if ( log_level == 3 ) {
	338	mat M;
	339	ldmat Lam;
	340	ldmat Vz;
	341	factorize ( M, Vz, Lam );
	342	logrec->L.log_vector ( logrec->ids ( 0 ), est_theta() );
	343	logrec->L.log_vector ( logrec->ids ( 1 ), cvectorize ( est_theta_cov().to_mat() ) );
	344	} else {
	345	epdf::log_write();
	346	}
	347
	348	}
	349	//!@}
318	350	};
319	351	UIREGISTER ( egiw );
…	…
328	360	where \f$\gamma=\sum_i \beta_i\f$.
329	361	*/
330		class eDirich: public eEF
331		{
332		protected:
333		//!sufficient statistics
334		vec beta;
335		public:
336		//!\name Constructors
337		//!@{
338
339		eDirich () : eEF () {};
340		eDirich (const eDirich &D0) : eEF () {set_parameters (D0.beta);};
341		eDirich (const vec &beta0) {set_parameters (beta0);};
342		void set_parameters (const vec &beta0) {
343		beta = beta0;
344		dim = beta.length();
345		}
346		//!@}
347
348		//! using sampling procedure from wikipedia
349		vec sample() const {
350		vec y(beta.length());
351		for (int i=0; i<beta.length(); i++){
352		GamRNG.setup(beta(i),1);
353		#pragma omp critical
354		y(i)=GamRNG();
355		}
356		return y/sum(y);
357		}
358
359		vec mean() const {return beta / sum (beta);};
360		vec variance() const {double gamma = sum (beta); return elem_mult (beta, (gamma-beta)) / (gammagamma (gamma + 1));}
361		//! In this instance, val is ...
362		double evallog_nn (const vec &val) const {
363		double tmp; tmp = (beta - 1) * log (val);
364		return tmp;
365		}
366
367		double lognc () const {
368		double tmp;
369		double gam = sum (beta);
370		double lgb = 0.0;
371		for (int i = 0;i < beta.length();i++) {lgb += lgamma (beta (i));}
372		tmp = lgb - lgamma (gam);
373		return tmp;
374		}
375
376		//!access function
377		vec& _beta() {return beta;}
378		/*! configuration structure
379		\code
380		class = 'eDirich';
381		beta = []; //parametr beta
382		\endcode
383		*/
384		void from_setting(const Setting &set){
385		epdf::from_setting(set);
386		UI::get(beta,set, "beta", UI::compulsory);
387		validate();
388		}
389		void validate() {
390		//check rv
391		dim = beta.length();
392		}
393		};
394		UIREGISTER(eDirich);
	362	class eDirich: public eEF {
	363	protected:
	364	//!sufficient statistics
	365	vec beta;
	366	public:
	367	//!\name Constructors
	368	//!@{
	369
	370	eDirich () : eEF () {};
	371	eDirich ( const eDirich &D0 ) : eEF () {
	372	set_parameters ( D0.beta );
	373	};
	374	eDirich ( const vec &beta0 ) {
	375	set_parameters ( beta0 );
	376	};
	377	void set_parameters ( const vec &beta0 ) {
	378	beta = beta0;
	379	dim = beta.length();
	380	}
	381	//!@}
	382
	383	//! using sampling procedure from wikipedia
	384	vec sample() const {
	385	vec y ( beta.length() );
	386	for ( int i = 0; i < beta.length(); i++ ) {
	387	GamRNG.setup ( beta ( i ), 1 );
	388	#pragma omp critical
	389	y ( i ) = GamRNG();
	390	}
	391	return y / sum ( y );
	392	}
	393
	394	vec mean() const {
	395	return beta / sum ( beta );
	396	};
	397	vec variance() const {
	398	double gamma = sum ( beta );
	399	return elem_mult ( beta, ( gamma - beta ) ) / ( gammagamma ( gamma + 1 ) );
	400	}
	401	//! In this instance, val is ...
	402	double evallog_nn ( const vec &val ) const {
	403	double tmp;
	404	tmp = ( beta - 1 ) * log ( val );
	405	return tmp;
	406	}
	407
	408	double lognc () const {
	409	double tmp;
	410	double gam = sum ( beta );
	411	double lgb = 0.0;
	412	for ( int i = 0; i < beta.length(); i++ ) {
	413	lgb += lgamma ( beta ( i ) );
	414	}
	415	tmp = lgb - lgamma ( gam );
	416	return tmp;
	417	}
	418
	419	//!access function
	420	vec& _beta() {
	421	return beta;
	422	}
	423	/*! configuration structure
	424	\code
	425	class = 'eDirich';
	426	beta = []; //parametr beta
	427	\endcode
	428	*/
	429	void from_setting ( const Setting &set ) {
	430	epdf::from_setting ( set );
	431	UI::get ( beta, set, "beta", UI::compulsory );
	432	validate();
	433	}
	434	void validate() {
	435	//check rv
	436	dim = beta.length();
	437	}
	438	};
	439	UIREGISTER ( eDirich );
395	440
396	441	/*! Random Walk on Dirichlet
397		Using simple assignment
	442	Using simple assignment
398	443	\f[ \beta = rvc / k + \beta_c \f]
399	444	hence, mean value = rvc, variance = (k+1)meanmean;
400
	445
401	446	The greater k is, the greater is the variance of the random walk;
402
	447
403	448	\f$ \beta_c \f$ is used as regularizing element to avoid corner cases, i.e. when one element of rvc is zero.
404	449	By default is it set to 0.1;
…	…
406	451
407	452	class mDirich: public pdf_internal<eDirich> {
408		protected:
409		//! constant \f$ k \f$ of the random walk
410		double k;
411		//! cache of beta_i
412		vec &_beta;
413		//! stabilizing coefficient \f$ \beta_c \f$
414		vec betac;
415		public:
416		mDirich(): pdf_internal<eDirich>(), _beta(iepdf._beta()){};
417		void condition (const vec &val) {_beta = val/k+betac; };
418		/*! Create Dirichlet random walk
419		\f[ f(rv\|rvc) = Di(rvc*k) \f]
420		from structure
421		\code
422		class = 'mDirich';
423		k = 1; // multiplicative constant k
424		--- optional ---
425		rv = RV({'name'},size) // description of RV
426		beta0 = []; // initial value of beta
427		betac = []; // initial value of beta
428		\endcode
429		*/
430		void from_setting (const Setting &set) {
431		pdf::from_setting (set); // reads rv and rvc
432		if (_rv()._dsize()>0){
433		rvc = _rv().copy_t(-1);
434		}
435		vec beta0;
436		if (!UI::get (beta0, set, "beta0", UI::optional)){
437		beta0 = ones(_rv()._dsize());
438		}
439		if (!UI::get (betac, set, "betac", UI::optional)){
440		betac = 0.1*ones(_rv()._dsize());
441		}
442		_beta = beta0;
443
444		UI::get (k, set, "k", UI::compulsory);
445		validate();
446		}
447		void validate() {
448		pdf_internal<eDirich>::validate();
449		bdm_assert(_beta.length()==betac.length(),"beta0 and betac are not compatible");
450		if (_rv()._dsize()>0){
451		bdm_assert( (_rv()._dsize()==dimension()) , "Size of rv does not match with beta");
452		}
453		dimc = _beta.length();
454		};
455		};
456		UIREGISTER(mDirich);
	453	protected:
	454	//! constant \f$ k \f$ of the random walk
	455	double k;
	456	//! cache of beta_i
	457	vec &_beta;
	458	//! stabilizing coefficient \f$ \beta_c \f$
	459	vec betac;
	460	public:
	461	mDirich() : pdf_internal<eDirich>(), _beta ( iepdf._beta() ) {};
	462	void condition ( const vec &val ) {
	463	_beta = val / k + betac;
	464	};
	465	/*! Create Dirichlet random walk
	466	\f[ f(rv\|rvc) = Di(rvc*k) \f]
	467	from structure
	468	\code
	469	class = 'mDirich';
	470	k = 1; // multiplicative constant k
	471	--- optional ---
	472	rv = RV({'name'},size) // description of RV
	473	beta0 = []; // initial value of beta
	474	betac = []; // initial value of beta
	475	\endcode
	476	*/
	477	void from_setting ( const Setting &set ) {
	478	pdf::from_setting ( set ); // reads rv and rvc
	479	if ( _rv()._dsize() > 0 ) {
	480	rvc = _rv().copy_t ( -1 );
	481	}
	482	vec beta0;
	483	if ( !UI::get ( beta0, set, "beta0", UI::optional ) ) {
	484	beta0 = ones ( _rv()._dsize() );
	485	}
	486	if ( !UI::get ( betac, set, "betac", UI::optional ) ) {
	487	betac = 0.1 * ones ( _rv()._dsize() );
	488	}
	489	_beta = beta0;
	490
	491	UI::get ( k, set, "k", UI::compulsory );
	492	validate();
	493	}
	494	void validate() {
	495	pdf_internal<eDirich>::validate();
	496	bdm_assert ( _beta.length() == betac.length(), "beta0 and betac are not compatible" );
	497	if ( _rv()._dsize() > 0 ) {
	498	bdm_assert ( ( _rv()._dsize() == dimension() ) , "Size of rv does not match with beta" );
	499	}
	500	dimc = _beta.length();
	501	};
	502	};
	503	UIREGISTER ( mDirich );
457	504
458	505	//! \brief Estimator for Multinomial density
459		class multiBM : public BMEF
460		{
461		protected:
462		//! Conjugate prior and posterior
463		eDirich est;
464		//! Pointer inside est to sufficient statistics
465		vec β
466		public:
467		//!Default constructor
468		multiBM () : BMEF (), est (), beta (est._beta()) {
469		if (beta.length() > 0) {last_lognc = est.lognc();}
470		else{last_lognc = 0.0;}
471		}
472		//!Copy constructor
473		multiBM (const multiBM &B) : BMEF (B), est (B.est), beta (est._beta()) {}
474		//! Sets sufficient statistics to match that of givefrom mB0
475		void set_statistics (const BM* mB0) {const multiBM* mB = dynamic_cast<const multiBM*> (mB0); beta = mB->beta;}
476		void bayes (const vec &yt, const vec &cond=empty_vec) {
477		if (frg < 1.0) {beta *= frg;last_lognc = est.lognc();}
478		beta += yt;
479		if (evalll) {ll = est.lognc() - last_lognc;}
480		}
481		double logpred (const vec &yt) const {
482		eDirich pred (est);
483		vec &beta = pred._beta();
484
485		double lll;
486		if (frg < 1.0)
487		{beta *= frg;lll = pred.lognc();}
488		else
489		if (evalll) {lll = last_lognc;}
490		else{lll = pred.lognc();}
491
492		beta += yt;
493		return pred.lognc() - lll;
494		}
495		void flatten (const BMEF* B) {
496		const multiBM* E = dynamic_cast<const multiBM*> (B);
497		// sum(beta) should be equal to sum(B.beta)
498		const vec &Eb = E->beta;//const_cast<multiBM*> ( E )->_beta();
499		beta *= (sum (Eb) / sum (beta));
500		if (evalll) {last_lognc = est.lognc();}
501		}
502		//! return correctly typed posterior (covariant return)
503		const eDirich& posterior() const {return est;};
504		//! constructor function
505		void set_parameters (const vec &beta0) {
506		est.set_parameters (beta0);
507		if (evalll) {last_lognc = est.lognc();}
508		}
509		void to_setting(Setting &set) const{
510		BMEF::to_setting(set);
511		Setting& prior= set.add("prior", Setting::TypeGroup);
512		est.to_setting(prior);
513		}
	506	class multiBM : public BMEF {
	507	protected:
	508	//! Conjugate prior and posterior
	509	eDirich est;
	510	//! Pointer inside est to sufficient statistics
	511	vec β
	512	public:
	513	//!Default constructor
	514	multiBM () : BMEF (), est (), beta ( est._beta() ) {
	515	if ( beta.length() > 0 ) {
	516	last_lognc = est.lognc();
	517	} else {
	518	last_lognc = 0.0;
	519	}
	520	}
	521	//!Copy constructor
	522	multiBM ( const multiBM &B ) : BMEF ( B ), est ( B.est ), beta ( est._beta() ) {}
	523	//! Sets sufficient statistics to match that of givefrom mB0
	524	void set_statistics ( const BM* mB0 ) {
	525	const multiBM* mB = dynamic_cast<const multiBM*> ( mB0 );
	526	beta = mB->beta;
	527	}
	528	void bayes ( const vec &yt, const vec &cond = empty_vec ) {
	529	if ( frg < 1.0 ) {
	530	beta *= frg;
	531	last_lognc = est.lognc();
	532	}
	533	beta += yt;
	534	if ( evalll ) {
	535	ll = est.lognc() - last_lognc;
	536	}
	537	}
	538	double logpred ( const vec &yt ) const {
	539	eDirich pred ( est );
	540	vec &beta = pred._beta();
	541
	542	double lll;
	543	if ( frg < 1.0 ) {
	544	beta *= frg;
	545	lll = pred.lognc();
	546	} else if ( evalll ) {
	547	lll = last_lognc;
	548	} else {
	549	lll = pred.lognc();
	550	}
	551
	552	beta += yt;
	553	return pred.lognc() - lll;
	554	}
	555	void flatten ( const BMEF* B ) {
	556	const multiBM* E = dynamic_cast<const multiBM*> ( B );
	557	// sum(beta) should be equal to sum(B.beta)
	558	const vec &Eb = E->beta;//const_cast<multiBM*> ( E )->_beta();
	559	beta *= ( sum ( Eb ) / sum ( beta ) );
	560	if ( evalll ) {
	561	last_lognc = est.lognc();
	562	}
	563	}
	564	//! return correctly typed posterior (covariant return)
	565	const eDirich& posterior() const {
	566	return est;
	567	};
	568	//! constructor function
	569	void set_parameters ( const vec &beta0 ) {
	570	est.set_parameters ( beta0 );
	571	if ( evalll ) {
	572	last_lognc = est.lognc();
	573	}
	574	}
	575	void to_setting ( Setting &set ) const {
	576	BMEF::to_setting ( set );
	577	Setting& prior = set.add ( "prior", Setting::TypeGroup );
	578	est.to_setting ( prior );
	579	}
514	580	};
515	581
…	…
523	589	*/
524	590
525		class egamma : public eEF
526		{
527		protected:
528		//! Vector \f$\alpha\f$
529		vec alpha;
530		//! Vector \f$\beta\f$
531		vec beta;
532		public :
533		//! \name Constructors
534		//!@{
535		egamma () : eEF (), alpha (0), beta (0) {};
536		egamma (const vec &a, const vec &b) {set_parameters (a, b);};
537		void set_parameters (const vec &a, const vec &b) {alpha = a, beta = b;dim = alpha.length();};
538		//!@}
539
540		vec sample() const;
541		double evallog (const vec &val) const;
542		double lognc () const;
543		//! Returns pointer to internal alpha. Potentially dengerous: use with care!
544		vec& _alpha() {return alpha;}
545		//! Returns pointer to internal beta. Potentially dengerous: use with care!
546		vec& _beta() {return beta;}
547		vec mean() const {return elem_div (alpha, beta);}
548		vec variance() const {return elem_div (alpha, elem_mult (beta, beta)); }
549
550		/*! Create Gamma density
551		\f[ f(rv\|rvc) = \Gamma(\alpha, \beta) \f]
552		from structure
553		\code
554		class = 'egamma';
555		alpha = [...]; // vector of alpha
556		beta = [...]; // vector of beta
557		rv = RV({'name'}) // description of RV
558		\endcode
559		*/
560		void from_setting (const Setting &set) {
561		epdf::from_setting (set); // reads rv
562		UI::get (alpha, set, "alpha", UI::compulsory);
563		UI::get (beta, set, "beta", UI::compulsory);
564		validate();
565		}
566		void validate() {
567		bdm_assert (alpha.length() == beta.length(), "parameters do not match");
568		dim = alpha.length();
569		}
570		};
571		UIREGISTER (egamma);
	591	class egamma : public eEF {
	592	protected:
	593	//! Vector \f$\alpha\f$
	594	vec alpha;
	595	//! Vector \f$\beta\f$
	596	vec beta;
	597	public :
	598	//! \name Constructors
	599	//!@{
	600	egamma () : eEF (), alpha ( 0 ), beta ( 0 ) {};
	601	egamma ( const vec &a, const vec &b ) {
	602	set_parameters ( a, b );
	603	};
	604	void set_parameters ( const vec &a, const vec &b ) {
	605	alpha = a, beta = b;
	606	dim = alpha.length();
	607	};
	608	//!@}
	609
	610	vec sample() const;
	611	double evallog ( const vec &val ) const;
	612	double lognc () const;
	613	//! Returns pointer to internal alpha. Potentially dengerous: use with care!
	614	vec& _alpha() {
	615	return alpha;
	616	}
	617	//! Returns pointer to internal beta. Potentially dengerous: use with care!
	618	vec& _beta() {
	619	return beta;
	620	}
	621	vec mean() const {
	622	return elem_div ( alpha, beta );
	623	}
	624	vec variance() const {
	625	return elem_div ( alpha, elem_mult ( beta, beta ) );
	626	}
	627
	628	/*! Create Gamma density
	629	\f[ f(rv\|rvc) = \Gamma(\alpha, \beta) \f]
	630	from structure
	631	\code
	632	class = 'egamma';
	633	alpha = [...]; // vector of alpha
	634	beta = [...]; // vector of beta
	635	rv = RV({'name'}) // description of RV
	636	\endcode
	637	*/
	638	void from_setting ( const Setting &set ) {
	639	epdf::from_setting ( set ); // reads rv
	640	UI::get ( alpha, set, "alpha", UI::compulsory );
	641	UI::get ( beta, set, "beta", UI::compulsory );
	642	validate();
	643	}
	644	void validate() {
	645	bdm_assert ( alpha.length() == beta.length(), "parameters do not match" );
	646	dim = alpha.length();
	647	}
	648	};
	649	UIREGISTER ( egamma );
572	650	SHAREDPTR ( egamma );
573	651
…	…
588	666	*/
589	667
590		class eigamma : public egamma
591		{
592		protected:
593		public :
594		//! \name Constructors
595		//! All constructors are inherited
596		//!@{
597		//!@}
598
599		vec sample() const {return 1.0 / egamma::sample();};
600		//! Returns poiter to alpha and beta. Potentially dangerous: use with care!
601		vec mean() const {return elem_div (beta, alpha - 1);}
602		vec variance() const {vec mea = mean(); return elem_div (elem_mult (mea, mea), alpha - 2);}
	668	class eigamma : public egamma {
	669	protected:
	670	public :
	671	//! \name Constructors
	672	//! All constructors are inherited
	673	//!@{
	674	//!@}
	675
	676	vec sample() const {
	677	return 1.0 / egamma::sample();
	678	};
	679	//! Returns poiter to alpha and beta. Potentially dangerous: use with care!
	680	vec mean() const {
	681	return elem_div ( beta, alpha - 1 );
	682	}
	683	vec variance() const {
	684	vec mea = mean();
	685	return elem_div ( elem_mult ( mea, mea ), alpha - 2 );
	686	}
603	687	};
604	688	/*
…	…
619	703	//! Uniform distributed density on a rectangular support
620	704
621		class euni: public epdf
622		{
623		protected:
	705	class euni: public epdf {
	706	protected:
624	707	//! lower bound on support
625		vec low;
	708	vec low;
626	709	//! upper bound on support
627		vec high;
	710	vec high;
628	711	//! internal
629		vec distance;
	712	vec distance;
630	713	//! normalizing coefficients
631		double nk;
	714	double nk;
632	715	//! cache of log( \c nk )
633		double lnk;
634		public:
635		//! \name Constructors
636		//!@{
637		euni () : epdf () {}
638		euni (const vec &low0, const vec &high0) {set_parameters (low0, high0);}
639		void set_parameters (const vec &low0, const vec &high0) {
640		distance = high0 - low0;
641		low = low0;
642		high = high0;
643		nk = prod (1.0 / distance);
644		lnk = log (nk);
645		dim = low.length();
646		}
647		//!@}
648
649		double evallog (const vec &val) const {
650		if (any (val < low) && any (val > high)) {return -inf;}
651		else return lnk;
652		}
653		vec sample() const {
654		vec smp (dim);
	716	double lnk;
	717	public:
	718	//! \name Constructors
	719	//!@{
	720	euni () : epdf () {}
	721	euni ( const vec &low0, const vec &high0 ) {
	722	set_parameters ( low0, high0 );
	723	}
	724	void set_parameters ( const vec &low0, const vec &high0 ) {
	725	distance = high0 - low0;
	726	low = low0;
	727	high = high0;
	728	nk = prod ( 1.0 / distance );
	729	lnk = log ( nk );
	730	dim = low.length();
	731	}
	732	//!@}
	733
	734	double evallog ( const vec &val ) const {
	735	if ( any ( val < low ) && any ( val > high ) ) {
	736	return -inf;
	737	} else return lnk;
	738	}
	739	vec sample() const {
	740	vec smp ( dim );
655	741	#pragma omp critical
656		UniRNG.sample_vector (dim , smp);
657		return low + elem_mult (distance, smp);
658		}
659		//! set values of \c low and \c high
660		vec mean() const {return (high -low) / 2.0;}
661		vec variance() const {return (pow (high, 2) + pow (low, 2) + elem_mult (high, low)) / 3.0;}
662		/*! Create Uniform density
663		\f[ f(rv) = U(low,high) \f]
664		from structure
665		\code
666		class = 'euni'
667		high = [...]; // vector of upper bounds
668		low = [...]; // vector of lower bounds
669		rv = RV({'name'}); // description of RV
670		\endcode
671		*/
672		void from_setting (const Setting &set) {
673		epdf::from_setting (set); // reads rv and rvc
674
675		UI::get (high, set, "high", UI::compulsory);
676		UI::get (low, set, "low", UI::compulsory);
677		set_parameters(low,high);
678		validate();
679		}
680		void validate() {
681		bdm_assert(high.length()==low.length(), "Incompatible high and low vectors");
682		dim = high.length();
683		bdm_assert (min (distance) > 0.0, "bad support");
684		}
685		};
686		UIREGISTER(euni);
	742	UniRNG.sample_vector ( dim , smp );
	743	return low + elem_mult ( distance, smp );
	744	}
	745	//! set values of \c low and \c high
	746	vec mean() const {
	747	return ( high - low ) / 2.0;
	748	}
	749	vec variance() const {
	750	return ( pow ( high, 2 ) + pow ( low, 2 ) + elem_mult ( high, low ) ) / 3.0;
	751	}
	752	/*! Create Uniform density
	753	\f[ f(rv) = U(low,high) \f]
	754	from structure
	755	\code
	756	class = 'euni'
	757	high = [...]; // vector of upper bounds
	758	low = [...]; // vector of lower bounds
	759	rv = RV({'name'}); // description of RV
	760	\endcode
	761	*/
	762	void from_setting ( const Setting &set ) {
	763	epdf::from_setting ( set ); // reads rv and rvc
	764
	765	UI::get ( high, set, "high", UI::compulsory );
	766	UI::get ( low, set, "low", UI::compulsory );
	767	set_parameters ( low, high );
	768	validate();
	769	}
	770	void validate() {
	771	bdm_assert ( high.length() == low.length(), "Incompatible high and low vectors" );
	772	dim = high.length();
	773	bdm_assert ( min ( distance ) > 0.0, "bad support" );
	774	}
	775	};
	776	UIREGISTER ( euni );
687	777
688	778	//! Uniform density with conditional mean value
689		class mguni : public pdf_internal<euni>{
	779	class mguni : public pdf_internal<euni> {
690	780	//! function of the mean value
691	781	shared_ptr<fnc> mean;
692	782	//! distance from mean to both sides
693	783	vec delta;
694		public:
695		~~void condition(const vec &cond)~~{
696		~~vec mea=mean->eval(cond~~);
697		~~iepdf.set_parameters(mea-delta,mea+delta~~);
698		}
699		//! load from
700		~~void from_setting(const Setting &set)~~{
701		~~pdf::from_setting(set~~); //reads rv and rvc
702		~~UI::get(delta,set,"delta",UI::compulsory~~);
703		~~mean = UI::build<fnc>(set,"mean",UI::compulsory~~);
704
705		~~iepdf.set_parameters(-delta,delta~~);
706		dimc = mean->dimensionc();
707		validate();
708		}
709		};
710		UIREGISTER~~(mguni~~);
	784	public:
	785	void condition ( const vec &cond ) {
	786	vec mea = mean->eval ( cond );
	787	iepdf.set_parameters ( mea - delta, mea + delta );
	788	}
	789	//! load from
	790	void from_setting ( const Setting &set ) {
	791	pdf::from_setting ( set ); //reads rv and rvc
	792	UI::get ( delta, set, "delta", UI::compulsory );
	793	mean = UI::build<fnc> ( set, "mean", UI::compulsory );
	794
	795	iepdf.set_parameters ( -delta, delta );
	796	dimc = mean->dimensionc();
	797	validate();
	798	}
	799	};
	800	UIREGISTER ( mguni );
711	801	/*!
712	802	\brief Normal distributed linear function with linear function of mean value;
…	…
715	805	*/
716	806	template < class sq_T, template <typename> class TEpdf = enorm >
717		class mlnorm : public pdf_internal< TEpdf<sq_T> >
718		{
719		protected:
720		//! Internal epdf that arise by conditioning on \c rvc
721		mat A;
722		//! Constant additive term
723		vec mu_const;
	807	class mlnorm : public pdf_internal< TEpdf<sq_T> > {
	808	protected:
	809	//! Internal epdf that arise by conditioning on \c rvc
	810	mat A;
	811	//! Constant additive term
	812	vec mu_const;
724	813	// vec& _mu; //cached epdf.mu; !!!!!! WHY NOT?
725		public:
726		//! \name Constructors
727		//!@{
728		mlnorm() : pdf_internal< TEpdf<sq_T> >() {};
729		~~mlnorm (const mat &A, const vec &mu0, const sq_T &R~~) : pdf_internal< TEpdf<sq_T> >() {
730		~~set_parameters (A, mu0, R~~);
731		}
732
733		//! Set \c A and \c R
734		~~void set_parameters (const mat &A0, const vec &mu0, const sq_T &R0) {~~
735		~~this->iepdf.set_parameters (zeros (A0.rows()), R0~~);
736		A = A0;
737		mu_const = mu0;
738		this->dimc = A0.cols();
739		}
740		//!@}
741		//! Set value of \c rvc . Result of this operation is stored in \c epdf use function \c _ep to access it.
742		~~void condition (const vec &cond~~) {
743		this->iepdf._mu() = A * cond + mu_const;
	814	public:
	815	//! \name Constructors
	816	//!@{
	817	mlnorm() : pdf_internal< TEpdf<sq_T> >() {};
	818	mlnorm ( const mat &A, const vec &mu0, const sq_T &R ) : pdf_internal< TEpdf<sq_T> >() {
	819	set_parameters ( A, mu0, R );
	820	}
	821
	822	//! Set \c A and \c R
	823	void set_parameters ( const mat &A0, const vec &mu0, const sq_T &R0 ) {
	824	this->iepdf.set_parameters ( zeros ( A0.rows() ), R0 );
	825	A = A0;
	826	mu_const = mu0;
	827	this->dimc = A0.cols();
	828	}
	829	//!@}
	830	//! Set value of \c rvc . Result of this operation is stored in \c epdf use function \c _ep to access it.
	831	void condition ( const vec &cond ) {
	832	this->iepdf._mu() = A * cond + mu_const;
744	833	//R is already assigned;
745		}
746
747		//!access function
748		const vec& _mu_const() const {return mu_const;}
749		//!access function
750		const mat& _A() const {return A;}
751		//!access function
752		mat _R() const { return this->iepdf._R().to_mat(); }
753		//!access function
754		sq_T __R() const { return this->iepdf._R(); }
755
756		//! Debug stream
757		template<typename sq_M>
758		friend std::ostream &operator<< (std::ostream &os, mlnorm<sq_M, enorm> &ml);
759
760		/*! Create Normal density with linear function of mean value
761		\f[ f(rv\|rvc) = N(A*rvc+const, R) \f]
762		from structure
763		\code
764		class = 'mlnorm<ldmat>', (OR) 'mlnorm<chmat>', (OR) 'mlnorm<fsqmat>';
765		A = []; // matrix or vector of appropriate dimension
766		const = []; // vector of constant term
767		R = []; // square matrix of appropriate dimension
768		\endcode
769		*/
770		void from_setting (const Setting &set) {
771		pdf::from_setting (set);
772
773		UI::get (A, set, "A", UI::compulsory);
774		UI::get (mu_const, set, "const", UI::compulsory);
775		mat R0;
776		UI::get (R0, set, "R", UI::compulsory);
777		set_parameters (A, mu_const, R0);
778		validate();
779		};
780		void validate() {
781		pdf_internal<TEpdf<sq_T> >::validate();
782		bdm_assert (A.rows() == mu_const.length(), "mlnorm: A vs. mu mismatch");
783		bdm_assert (A.rows() == _R().rows(), "mlnorm: A vs. R mismatch");
784
785		}
786		};
787		UIREGISTER2 (mlnorm,ldmat);
	834	}
	835
	836	//!access function
	837	const vec& _mu_const() const {
	838	return mu_const;
	839	}
	840	//!access function
	841	const mat& _A() const {
	842	return A;
	843	}
	844	//!access function
	845	mat _R() const {
	846	return this->iepdf._R().to_mat();
	847	}
	848	//!access function
	849	sq_T __R() const {
	850	return this->iepdf._R();
	851	}
	852
	853	//! Debug stream
	854	template<typename sq_M>
	855	friend std::ostream &operator<< ( std::ostream &os, mlnorm<sq_M, enorm> &ml );
	856
	857	/*! Create Normal density with linear function of mean value
	858	\f[ f(rv\|rvc) = N(A*rvc+const, R) \f]
	859	from structure
	860	\code
	861	class = 'mlnorm<ldmat>', (OR) 'mlnorm<chmat>', (OR) 'mlnorm<fsqmat>';
	862	A = []; // matrix or vector of appropriate dimension
	863	const = []; // vector of constant term
	864	R = []; // square matrix of appropriate dimension
	865	\endcode
	866	*/
	867	void from_setting ( const Setting &set ) {
	868	pdf::from_setting ( set );
	869
	870	UI::get ( A, set, "A", UI::compulsory );
	871	UI::get ( mu_const, set, "const", UI::compulsory );
	872	mat R0;
	873	UI::get ( R0, set, "R", UI::compulsory );
	874	set_parameters ( A, mu_const, R0 );
	875	validate();
	876	};
	877	void validate() {
	878	pdf_internal<TEpdf<sq_T> >::validate();
	879	bdm_assert ( A.rows() == mu_const.length(), "mlnorm: A vs. mu mismatch" );
	880	bdm_assert ( A.rows() == _R().rows(), "mlnorm: A vs. R mismatch" );
	881
	882	}
	883	};
	884	UIREGISTER2 ( mlnorm, ldmat );
788	885	SHAREDPTR2 ( mlnorm, ldmat );
789		UIREGISTER2 (~~mlnorm,fsqmat~~);
	886	UIREGISTER2 ( mlnorm, fsqmat );
790	887	SHAREDPTR2 ( mlnorm, fsqmat );
791		UIREGISTER2 (~~mlnorm, chmat~~);
	888	UIREGISTER2 ( mlnorm, chmat );
792	889	SHAREDPTR2 ( mlnorm, chmat );
793	890
794	891	//! pdf with general function for mean value
795	892	template<class sq_T>
796		class mgnorm : public pdf_internal< enorm< sq_T > >
797		{
798		private:
	893	class mgnorm : public pdf_internal< enorm< sq_T > > {
	894	private:
799	895	// vec μ WHY NOT?
800		shared_ptr<fnc> g;
801
802		public:
803		//!default constructor
804		mgnorm() : pdf_internal<enorm<sq_T> >() { }
805		//!set mean function
806		~~inline void set_parameters (const shared_ptr<fnc> &g0, const sq_T &R0~~);
807		~~inline void condition (const vec &cond~~);
808
809
810		~~/*! Create Normal density with given function of mean value~~
811		\f[ f(rv\|rvc) = N( g(rvc), R) \f]
812		from structure
813		\code
814		class = 'mgnorm';
815		g.class = 'fnc'; // function for mean value evolution
816		g._fields_of_fnc = ...;
817
818		R = [1, 0; // covariance matrix
819		0, 1];
820		--OR --
821		dR = [1, 1]; // diagonal of cavariance matrix
822
823		rv = RV({'name'}) // description of RV
824		rvc = RV({'name'}) // description of RV in condition
825		\endcode
826		*/
827
828		~~void from_setting (const Setting &set~~) {
829		~~pdf::from_setting(set~~);
830		~~shared_ptr<fnc> g = UI::build<fnc> (set, "g", UI::compulsory~~);
831
832		mat R;
833		vec dR;
834		~~if (UI::get (dR, set, "dR")~~)
835		~~R = diag (dR~~);
836		else
837		~~UI::get (R, set, "R", UI::compulsory~~);
838
839		~~set_parameters (g, R~~);
840		validate();
841		}
842		void validate() {
843		~~bdm_assert(g->dimension()==this->dimension(),"incompatible function"~~);
844		}
845		};
846
847		UIREGISTER2 (~~mgnorm, chmat~~);
	896	shared_ptr<fnc> g;
	897
	898	public:
	899	//!default constructor
	900	mgnorm() : pdf_internal<enorm<sq_T> >() { }
	901	//!set mean function
	902	inline void set_parameters ( const shared_ptr<fnc> &g0, const sq_T &R0 );
	903	inline void condition ( const vec &cond );
	904
	905
	906	/*! Create Normal density with given function of mean value
	907	\f[ f(rv\|rvc) = N( g(rvc), R) \f]
	908	from structure
	909	\code
	910	class = 'mgnorm';
	911	g.class = 'fnc'; // function for mean value evolution
	912	g._fields_of_fnc = ...;
	913
	914	R = [1, 0; // covariance matrix
	915	0, 1];
	916	--OR --
	917	dR = [1, 1]; // diagonal of cavariance matrix
	918
	919	rv = RV({'name'}) // description of RV
	920	rvc = RV({'name'}) // description of RV in condition
	921	\endcode
	922	*/
	923
	924	void from_setting ( const Setting &set ) {
	925	pdf::from_setting ( set );
	926	shared_ptr<fnc> g = UI::build<fnc> ( set, "g", UI::compulsory );
	927
	928	mat R;
	929	vec dR;
	930	if ( UI::get ( dR, set, "dR" ) )
	931	R = diag ( dR );
	932	else
	933	UI::get ( R, set, "R", UI::compulsory );
	934
	935	set_parameters ( g, R );
	936	validate();
	937	}
	938	void validate() {
	939	bdm_assert ( g->dimension() == this->dimension(), "incompatible function" );
	940	}
	941	};
	942
	943	UIREGISTER2 ( mgnorm, chmat );
848	944	SHAREDPTR2 ( mgnorm, chmat );
849	945
…	…
856	952	Perhaps a moment-matching technique?
857	953	*/
858		class mlstudent : public mlnorm<ldmat, enorm>
859		{
860		protected:
861		//! Variable \f$ \Lambda \f$ from theory
862		ldmat Lambda;
863		//! Reference to variable \f$ R \f$
864		ldmat &_R;
865		//! Variable \f$ R_e \f$
866		ldmat Re;
867		public:
868		mlstudent () : mlnorm<ldmat, enorm> (),
869		Lambda (), _R (iepdf._R()) {}
870		//! constructor function
871		void set_parameters (const mat &A0, const vec &mu0, const ldmat &R0, const ldmat& Lambda0) {
872		iepdf.set_parameters (mu0, R0);// was Lambda, why?
873		A = A0;
874		mu_const = mu0;
875		Re = R0;
876		Lambda = Lambda0;
877		}
878		void condition (const vec &cond) {
879		if (cond.length()>0) {
880		iepdf._mu() = A * cond + mu_const;
881		} else {
882		iepdf._mu() = mu_const;
883		}
884		double zeta;
885		//ugly hack!
886		if ( (cond.length() + 1) == Lambda.rows()) {
887		zeta = Lambda.invqform (concat (cond, vec_1 (1.0)));
888		} else {
889		zeta = Lambda.invqform (cond);
890		}
891		_R = Re;
892		_R *= (1 + zeta);// / ( nu ); << nu is in Re!!!!!!
893		};
894
895		void validate() {
896		bdm_assert (A.rows() == mu_const.length(), "mlstudent: A vs. mu mismatch");
897		bdm_assert (_R.rows() == A.rows(), "mlstudent: A vs. R mismatch");
898
899		}
	954	class mlstudent : public mlnorm<ldmat, enorm> {
	955	protected:
	956	//! Variable \f$ \Lambda \f$ from theory
	957	ldmat Lambda;
	958	//! Reference to variable \f$ R \f$
	959	ldmat &_R;
	960	//! Variable \f$ R_e \f$
	961	ldmat Re;
	962	public:
	963	mlstudent () : mlnorm<ldmat, enorm> (),
	964	Lambda (), _R ( iepdf._R() ) {}
	965	//! constructor function
	966	void set_parameters ( const mat &A0, const vec &mu0, const ldmat &R0, const ldmat& Lambda0 ) {
	967	iepdf.set_parameters ( mu0, R0 );// was Lambda, why?
	968	A = A0;
	969	mu_const = mu0;
	970	Re = R0;
	971	Lambda = Lambda0;
	972	}
	973	void condition ( const vec &cond ) {
	974	if ( cond.length() > 0 ) {
	975	iepdf._mu() = A * cond + mu_const;
	976	} else {
	977	iepdf._mu() = mu_const;
	978	}
	979	double zeta;
	980	//ugly hack!
	981	if ( ( cond.length() + 1 ) == Lambda.rows() ) {
	982	zeta = Lambda.invqform ( concat ( cond, vec_1 ( 1.0 ) ) );
	983	} else {
	984	zeta = Lambda.invqform ( cond );
	985	}
	986	_R = Re;
	987	_R *= ( 1 + zeta );// / ( nu ); << nu is in Re!!!!!!
	988	};
	989
	990	void validate() {
	991	bdm_assert ( A.rows() == mu_const.length(), "mlstudent: A vs. mu mismatch" );
	992	bdm_assert ( _R.rows() == A.rows(), "mlstudent: A vs. R mismatch" );
	993
	994	}
900	995	};
901	996	/*!
…	…
908	1003	The standard deviation of the walk is then: \f$\mu/\sqrt(k)\f$.
909	1004	*/
910		class mgamma : public pdf_internal<egamma>
911		{
912		protected:
913
914		//! Constant \f$k\f$
915		double k;
916
917		//! cache of iepdf.beta
918		vec &_beta;
919
920		public:
921		//! Constructor
922		mgamma() : pdf_internal<egamma>(), k (0),
923		_beta (iepdf._beta()) {
924		}
925
926		//! Set value of \c k
927		void set_parameters (double k, const vec &beta0);
928
929		void condition (const vec &val) {_beta = k / val;};
930		/*! Create Gamma density with conditional mean value
931		\f[ f(rv\|rvc) = \Gamma(k, k/rvc) \f]
932		from structure
933		\code
934		class = 'mgamma';
935		beta = [...]; // vector of initial alpha
936		k = 1.1; // multiplicative constant k
937		rv = RV({'name'}) // description of RV
938		rvc = RV({'name'}) // description of RV in condition
939		\endcode
940		*/
941		void from_setting (const Setting &set) {
942		pdf::from_setting (set); // reads rv and rvc
943		vec betatmp; // ugly but necessary
944		UI::get (betatmp, set, "beta", UI::compulsory);
945		UI::get (k, set, "k", UI::compulsory);
946		set_parameters (k, betatmp);
947		validate();
948		}
949		void validate() {
950		pdf_internal<egamma>::validate();
951
952		dim = _beta.length();
953		dimc = _beta.length();
954		}
955		};
956		UIREGISTER (mgamma);
957		SHAREDPTR (mgamma);
	1005	class mgamma : public pdf_internal<egamma> {
	1006	protected:
	1007
	1008	//! Constant \f$k\f$
	1009	double k;
	1010
	1011	//! cache of iepdf.beta
	1012	vec &_beta;
	1013
	1014	public:
	1015	//! Constructor
	1016	mgamma() : pdf_internal<egamma>(), k ( 0 ),
	1017	_beta ( iepdf._beta() ) {
	1018	}
	1019
	1020	//! Set value of \c k
	1021	void set_parameters ( double k, const vec &beta0 );
	1022
	1023	void condition ( const vec &val ) {
	1024	_beta = k / val;
	1025	};
	1026	/*! Create Gamma density with conditional mean value
	1027	\f[ f(rv\|rvc) = \Gamma(k, k/rvc) \f]
	1028	from structure
	1029	\code
	1030	class = 'mgamma';
	1031	beta = [...]; // vector of initial alpha
	1032	k = 1.1; // multiplicative constant k
	1033	rv = RV({'name'}) // description of RV
	1034	rvc = RV({'name'}) // description of RV in condition
	1035	\endcode
	1036	*/
	1037	void from_setting ( const Setting &set ) {
	1038	pdf::from_setting ( set ); // reads rv and rvc
	1039	vec betatmp; // ugly but necessary
	1040	UI::get ( betatmp, set, "beta", UI::compulsory );
	1041	UI::get ( k, set, "k", UI::compulsory );
	1042	set_parameters ( k, betatmp );
	1043	validate();
	1044	}
	1045	void validate() {
	1046	pdf_internal<egamma>::validate();
	1047
	1048	dim = _beta.length();
	1049	dimc = _beta.length();
	1050	}
	1051	};
	1052	UIREGISTER ( mgamma );
	1053	SHAREDPTR ( mgamma );
958	1054
959	1055	/*!
…	…
966	1062	The standard deviation of the walk is then: \f$ \mu/\sqrt(k)\f$.
967	1063	*/
968		class migamma : public pdf_internal<eigamma>
969		{
970		protected:
971		//! Constant \f$k\f$
972		double k;
973
974		//! cache of iepdf.alpha
975		vec &_alpha;
976
977		//! cache of iepdf.beta
978		vec &_beta;
979
980		public:
981		//! \name Constructors
982		//!@{
983		migamma() : pdf_internal<eigamma>(),
984		k (0),
985		_alpha (iepdf._alpha()),
986		_beta (iepdf._beta()) {
987		}
988
989		migamma (const migamma &m) : pdf_internal<eigamma>(),
990		k (0),
991		_alpha (iepdf._alpha()),
992		_beta (iepdf._beta()) {
993		}
994		//!@}
995
996		//! Set value of \c k
997		void set_parameters (int len, double k0) {
998		k = k0;
999		iepdf.set_parameters ( (1.0 / (kk) + 2.0) ones (len) /alpha/, ones (len) /beta/);
1000		dimc = dimension();
1001		};
1002		void condition (const vec &val) {
1003		_beta = elem_mult (val, (_alpha - 1.0));
1004		};
	1064	class migamma : public pdf_internal<eigamma> {
	1065	protected:
	1066	//! Constant \f$k\f$
	1067	double k;
	1068
	1069	//! cache of iepdf.alpha
	1070	vec &_alpha;
	1071
	1072	//! cache of iepdf.beta
	1073	vec &_beta;
	1074
	1075	public:
	1076	//! \name Constructors
	1077	//!@{
	1078	migamma() : pdf_internal<eigamma>(),
	1079	k ( 0 ),
	1080	_alpha ( iepdf._alpha() ),
	1081	_beta ( iepdf._beta() ) {
	1082	}
	1083
	1084	migamma ( const migamma &m ) : pdf_internal<eigamma>(),
	1085	k ( 0 ),
	1086	_alpha ( iepdf._alpha() ),
	1087	_beta ( iepdf._beta() ) {
	1088	}
	1089	//!@}
	1090
	1091	//! Set value of \c k
	1092	void set_parameters ( int len, double k0 ) {
	1093	k = k0;
	1094	iepdf.set_parameters ( ( 1.0 / ( kk ) + 2.0 ) ones ( len ) /alpha/, ones ( len ) /beta/ );
	1095	dimc = dimension();
	1096	};
	1097	void condition ( const vec &val ) {
	1098	_beta = elem_mult ( val, ( _alpha - 1.0 ) );
	1099	};
1005	1100	};
1006	1101
…	…
1017	1112	The standard deviation of the walk is then: \f$\mu/\sqrt(k)\f$.
1018	1113	*/
1019		class mgamma_fix : public mgamma
1020		{
1021		protected:
1022		//! parameter l
1023		double l;
1024		//! reference vector
1025		vec refl;
1026		public:
1027		//! Constructor
1028		mgamma_fix () : mgamma (), refl () {};
1029		//! Set value of \c k
1030		void set_parameters (double k0 , vec ref0, double l0) {
1031		mgamma::set_parameters (k0, ref0);
1032		refl = pow (ref0, 1.0 - l0);l = l0;
1033		dimc = dimension();
1034		};
1035
1036		void condition (const vec &val) {vec mean = elem_mult (refl, pow (val, l)); _beta = k / mean;};
	1114	class mgamma_fix : public mgamma {
	1115	protected:
	1116	//! parameter l
	1117	double l;
	1118	//! reference vector
	1119	vec refl;
	1120	public:
	1121	//! Constructor
	1122	mgamma_fix () : mgamma (), refl () {};
	1123	//! Set value of \c k
	1124	void set_parameters ( double k0 , vec ref0, double l0 ) {
	1125	mgamma::set_parameters ( k0, ref0 );
	1126	refl = pow ( ref0, 1.0 - l0 );
	1127	l = l0;
	1128	dimc = dimension();
	1129	};
	1130
	1131	void condition ( const vec &val ) {
	1132	vec mean = elem_mult ( refl, pow ( val, l ) );
	1133	_beta = k / mean;
	1134	};
1037	1135	};
1038	1136
…	…
1050	1148	The standard deviation of the walk is then: \f$\mu/\sqrt(k)\f$.
1051	1149	*/
1052		class migamma_ref : public migamma
1053		{
1054		protected:
1055		//! parameter l
1056		double l;
1057		//! reference vector
1058		vec refl;
1059		public:
1060		//! Constructor
1061		migamma_ref () : migamma (), refl () {};
1062		//! Set value of \c k
1063		void set_parameters (double k0 , vec ref0, double l0) {
1064		migamma::set_parameters (ref0.length(), k0);
1065		refl = pow (ref0, 1.0 - l0);
1066		l = l0;
1067		dimc = dimension();
1068		};
1069
1070		void condition (const vec &val) {
1071		vec mean = elem_mult (refl, pow (val, l));
1072		migamma::condition (mean);
1073		};
1074
1075
1076		/*! Create inverse-Gamma density with conditional mean value
1077		\f[ f(rv\|rvc) = i\Gamma(k, k/(rvc^l \circ ref^{(1-l)}) \f]
1078		from structure
1079		\code
1080		class = 'migamma_ref';
1081		ref = [1e-5; 1e-5; 1e-2 1e-3]; // reference vector
1082		l = 0.999; // constant l
1083		k = 0.1; // constant k
1084		rv = RV({'name'}) // description of RV
1085		rvc = RV({'name'}) // description of RV in condition
1086		\endcode
1087		*/
1088		void from_setting (const Setting &set);
1089
1090		// TODO dodelat void to_setting( Setting &set ) const;
1091		};
1092
1093
1094		UIREGISTER (migamma_ref);
1095		SHAREDPTR (migamma_ref);
	1150	class migamma_ref : public migamma {
	1151	protected:
	1152	//! parameter l
	1153	double l;
	1154	//! reference vector
	1155	vec refl;
	1156	public:
	1157	//! Constructor
	1158	migamma_ref () : migamma (), refl () {};
	1159	//! Set value of \c k
	1160	void set_parameters ( double k0 , vec ref0, double l0 ) {
	1161	migamma::set_parameters ( ref0.length(), k0 );
	1162	refl = pow ( ref0, 1.0 - l0 );
	1163	l = l0;
	1164	dimc = dimension();
	1165	};
	1166
	1167	void condition ( const vec &val ) {
	1168	vec mean = elem_mult ( refl, pow ( val, l ) );
	1169	migamma::condition ( mean );
	1170	};
	1171
	1172
	1173	/*! Create inverse-Gamma density with conditional mean value
	1174	\f[ f(rv\|rvc) = i\Gamma(k, k/(rvc^l \circ ref^{(1-l)}) \f]
	1175	from structure
	1176	\code
	1177	class = 'migamma_ref';
	1178	ref = [1e-5; 1e-5; 1e-2 1e-3]; // reference vector
	1179	l = 0.999; // constant l
	1180	k = 0.1; // constant k
	1181	rv = RV({'name'}) // description of RV
	1182	rvc = RV({'name'}) // description of RV in condition
	1183	\endcode
	1184	*/
	1185	void from_setting ( const Setting &set );
	1186
	1187	// TODO dodelat void to_setting( Setting &set ) const;
	1188	};
	1189
	1190
	1191	UIREGISTER ( migamma_ref );
	1192	SHAREDPTR ( migamma_ref );
1096	1193
1097	1194	/*! Log-Normal probability density
…	…
1105	1202	Function from_setting loads mu and R in the same way as it does for enorm<>!
1106	1203	*/
1107		class elognorm: public enorm<ldmat>
1108		{
1109		public:
1110		vec sample() const {return exp (enorm<ldmat>::sample());};
1111		vec mean() const {vec var = enorm<ldmat>::variance();return exp (mu - 0.5*var);};
	1204	class elognorm: public enorm<ldmat> {
	1205	public:
	1206	vec sample() const {
	1207	return exp ( enorm<ldmat>::sample() );
	1208	};
	1209	vec mean() const {
	1210	vec var = enorm<ldmat>::variance();
	1211	return exp ( mu - 0.5*var );
	1212	};
1112	1213
1113	1214	};
…	…
1119	1220
1120	1221	*/
1121		class mlognorm : public pdf_internal<elognorm>
1122		{
1123		protected:
1124		//! parameter 1/2*sigma^2
1125		double sig2;
1126
1127		//! access
1128		vec μ
1129		public:
1130		//! Constructor
1131		mlognorm() : pdf_internal<elognorm>(),
1132		sig2 (0),
1133		mu (iepdf._mu()) {
1134		}
1135
1136		//! Set value of \c k
1137		void set_parameters (int size, double k) {
1138		sig2 = 0.5 * log (k * k + 1);
1139		iepdf.set_parameters (zeros (size), 2sig2eye (size));
1140
1141		dimc = size;
1142		};
1143
1144		void condition (const vec &val) {
1145		mu = log (val) - sig2;//elem_mult ( refl,pow ( val,l ) );
1146		};
1147
1148		/*! Create logNormal random Walk
1149		\f[ f(rv\|rvc) = log\mathcal{N}( \log(rvc)-0.5\log(k^2+1), k I) \f]
1150		from structure
1151		\code
1152		class = 'mlognorm';
1153		k = 0.1; // "variance" k
1154		mu0 = 0.1; // Initial value of mean
1155		rv = RV({'name'}) // description of RV
1156		rvc = RV({'name'}) // description of RV in condition
1157		\endcode
1158		*/
1159		void from_setting (const Setting &set);
1160
1161		// TODO dodelat void to_setting( Setting &set ) const;
1162
1163		};
1164
1165		UIREGISTER (mlognorm);
1166		SHAREDPTR (mlognorm);
	1222	class mlognorm : public pdf_internal<elognorm> {
	1223	protected:
	1224	//! parameter 1/2*sigma^2
	1225	double sig2;
	1226
	1227	//! access
	1228	vec μ
	1229	public:
	1230	//! Constructor
	1231	mlognorm() : pdf_internal<elognorm>(),
	1232	sig2 ( 0 ),
	1233	mu ( iepdf._mu() ) {
	1234	}
	1235
	1236	//! Set value of \c k
	1237	void set_parameters ( int size, double k ) {
	1238	sig2 = 0.5 * log ( k * k + 1 );
	1239	iepdf.set_parameters ( zeros ( size ), 2sig2eye ( size ) );
	1240
	1241	dimc = size;
	1242	};
	1243
	1244	void condition ( const vec &val ) {
	1245	mu = log ( val ) - sig2;//elem_mult ( refl,pow ( val,l ) );
	1246	};
	1247
	1248	/*! Create logNormal random Walk
	1249	\f[ f(rv\|rvc) = log\mathcal{N}( \log(rvc)-0.5\log(k^2+1), k I) \f]
	1250	from structure
	1251	\code
	1252	class = 'mlognorm';
	1253	k = 0.1; // "variance" k
	1254	mu0 = 0.1; // Initial value of mean
	1255	rv = RV({'name'}) // description of RV
	1256	rvc = RV({'name'}) // description of RV in condition
	1257	\endcode
	1258	*/
	1259	void from_setting ( const Setting &set );
	1260
	1261	// TODO dodelat void to_setting( Setting &set ) const;
	1262
	1263	};
	1264
	1265	UIREGISTER ( mlognorm );
	1266	SHAREDPTR ( mlognorm );
1167	1267
1168	1268	/*! inverse Wishart density defined on Choleski decomposition
1169	1269
1170	1270	*/
1171		class eWishartCh : public epdf
1172		{
1173		protected:
1174		//! Upper-Triagle of Choleski decomposition of \f$ \Psi \f$
1175		chmat Y;
1176		//! dimension of matrix \f$ \Psi \f$
1177		int p;
1178		//! degrees of freedom \f$ \nu \f$
1179		double delta;
1180		public:
1181		//! Set internal structures
1182		void set_parameters (const mat &Y0, const double delta0) {Y = chmat (Y0);delta = delta0; p = Y.rows(); dim = p * p; }
1183		//! Set internal structures
1184		void set_parameters (const chmat &Y0, const double delta0) {Y = Y0;delta = delta0; p = Y.rows(); dim = p * p; }
1185		//! Sample matrix argument
1186		mat sample_mat() const {
1187		mat X = zeros (p, p);
1188
1189		//sample diagonal
1190		for (int i = 0;i < p;i++) {
1191		GamRNG.setup (0.5* (delta - i) , 0.5); // no +1 !! index if from 0
	1271	class eWishartCh : public epdf {
	1272	protected:
	1273	//! Upper-Triagle of Choleski decomposition of \f$ \Psi \f$
	1274	chmat Y;
	1275	//! dimension of matrix \f$ \Psi \f$
	1276	int p;
	1277	//! degrees of freedom \f$ \nu \f$
	1278	double delta;
	1279	public:
	1280	//! Set internal structures
	1281	void set_parameters ( const mat &Y0, const double delta0 ) {
	1282	Y = chmat ( Y0 );
	1283	delta = delta0;
	1284	p = Y.rows();
	1285	dim = p * p;
	1286	}
	1287	//! Set internal structures
	1288	void set_parameters ( const chmat &Y0, const double delta0 ) {
	1289	Y = Y0;
	1290	delta = delta0;
	1291	p = Y.rows();
	1292	dim = p * p;
	1293	}
	1294	//! Sample matrix argument
	1295	mat sample_mat() const {
	1296	mat X = zeros ( p, p );
	1297
	1298	//sample diagonal
	1299	for ( int i = 0; i < p; i++ ) {
	1300	GamRNG.setup ( 0.5* ( delta - i ) , 0.5 ); // no +1 !! index if from 0
1192	1301	#pragma omp critical
1193		X (i, i) = sqrt (GamRNG());
	1302	X ( i, i ) = sqrt ( GamRNG() );
	1303	}
	1304	//do the rest
	1305	for ( int i = 0; i < p; i++ ) {
	1306	for ( int j = i + 1; j < p; j++ ) {
	1307	#pragma omp critical
	1308	X ( i, j ) = NorRNG.sample();
1194	1309	}
1195		~~//do the rest~~
1196		~~for (int i = 0;i < p;i++) {~~
1197		~~for (int j = i + 1;j < p;j++) {~~
1198		#pragma omp critical
1199		~~X (i, j) = NorRNG.sample(~~);
1200		}
1201		}
1202		~~return X*Y._Ch();// return upper triangular part of the decomposition~~
1203		}
1204		~~vec sample () const {~~
1205		~~return vec (sample_mat()._data(), p*p);~~
1206		}
1207		~~//! fast access function y0 will be copied into Y.Ch.~~
1208		~~void setY (const mat &Ch0) {copy_vector (dim, Ch0._data(), Y._Ch()._data());~~}
1209		~~//! fast access function y0 will be copied into Y.Ch.~~
1210		~~void _setY (const vec &ch0) {copy_vector (dim, ch0._data(), Y._Ch()._data()); }~~
1211		~~//! access function~~
1212		~~const chmat& getY() const {return Y;~~}
	1310	}
	1311	return X*Y._Ch();// return upper triangular part of the decomposition
	1312	}
	1313	vec sample () const {
	1314	return vec ( sample_mat()._data(), p*p );
	1315	}
	1316	//! fast access function y0 will be copied into Y.Ch.
	1317	void setY ( const mat &Ch0 ) {
	1318	copy_vector ( dim, Ch0._data(), Y._Ch()._data() );
	1319	}
	1320	//! fast access function y0 will be copied into Y.Ch.
	1321	void _setY ( const vec &ch0 ) {
	1322	copy_vector ( dim, ch0._data(), Y._Ch()._data() );
	1323	}
	1324	//! access function
	1325	const chmat& getY() const {
	1326	return Y;
	1327	}
1213	1328	};
1214	1329
…	…
1216	1331	/*! Being computed by conversion from `standard' Wishart
1217	1332	*/
1218		class eiWishartCh: public epdf
1219		{
1220		protected:
1221		//! Internal instance of Wishart density
1222		eWishartCh W;
1223		//! size of Ch
1224		int p;
1225		//! parameter delta
1226		double delta;
1227		public:
1228		//! constructor function
1229		void set_parameters (const mat &Y0, const double delta0) {
1230		delta = delta0;
1231		W.set_parameters (inv (Y0), delta0);
1232		dim = W.dimension(); p = Y0.rows();
1233		}
1234		vec sample() const {mat iCh; iCh = inv (W.sample_mat()); return vec (iCh._data(), dim);}
1235		//! access function
1236		void _setY (const vec &y0) {
1237		mat Ch (p, p);
1238		mat iCh (p, p);
1239		copy_vector (dim, y0._data(), Ch._data());
1240
1241		iCh = inv (Ch);
1242		W.setY (iCh);
1243		}
1244		virtual double evallog (const vec &val) const {
1245		chmat X (p);
1246		const chmat& Y = W.getY();
1247
1248		copy_vector (p*p, val._data(), X._Ch()._data());
1249		chmat iX (p);X.inv (iX);
1250		// compute
	1333	class eiWishartCh: public epdf {
	1334	protected:
	1335	//! Internal instance of Wishart density
	1336	eWishartCh W;
	1337	//! size of Ch
	1338	int p;
	1339	//! parameter delta
	1340	double delta;
	1341	public:
	1342	//! constructor function
	1343	void set_parameters ( const mat &Y0, const double delta0 ) {
	1344	delta = delta0;
	1345	W.set_parameters ( inv ( Y0 ), delta0 );
	1346	dim = W.dimension();
	1347	p = Y0.rows();
	1348	}
	1349	vec sample() const {
	1350	mat iCh;
	1351	iCh = inv ( W.sample_mat() );
	1352	return vec ( iCh._data(), dim );
	1353	}
	1354	//! access function
	1355	void _setY ( const vec &y0 ) {
	1356	mat Ch ( p, p );
	1357	mat iCh ( p, p );
	1358	copy_vector ( dim, y0._data(), Ch._data() );
	1359
	1360	iCh = inv ( Ch );
	1361	W.setY ( iCh );
	1362	}
	1363	virtual double evallog ( const vec &val ) const {
	1364	chmat X ( p );
	1365	const chmat& Y = W.getY();
	1366
	1367	copy_vector ( p*p, val._data(), X._Ch()._data() );
	1368	chmat iX ( p );
	1369	X.inv ( iX );
	1370	// compute
1251	1371	// \frac{ \|\Psi\|^{m/2}\|X\|^{-(m+p+1)/2}e^{-tr(\Psi X^{-1})/2} }{ 2^{mp/2}\Gamma_p(m/2)},
1252		mat M = Y.to_mat() * iX.to_mat();
1253
1254		~~double log1 = 0.5 * p * (2 * Y.logdet()) - 0.5 * (delta + p + 1) * (2 * X.logdet()) - 0.5 * trace (M~~);
1255		//Fixme! Multivariate gamma omitted!! it is ok for sampling, but not otherwise!!
1256
1257		/* if (0) {
1258		mat XX=X.to_mat();
1259		mat YY=Y.to_mat();
1260
1261		double log2 = 0.5plog(det(YY))-0.5(delta+p+1)log(det(XX))-0.5trace(YYinv(XX));
1262		cout << log1 << "," << log2 << endl;
1263		}*/
1264		return log1;
1265		};
	1372	mat M = Y.to_mat() * iX.to_mat();
	1373
	1374	double log1 = 0.5 * p * ( 2 * Y.logdet() ) - 0.5 * ( delta + p + 1 ) * ( 2 * X.logdet() ) - 0.5 * trace ( M );
	1375	//Fixme! Multivariate gamma omitted!! it is ok for sampling, but not otherwise!!
	1376
	1377	/* if (0) {
	1378	mat XX=X.to_mat();
	1379	mat YY=Y.to_mat();
	1380
	1381	double log2 = 0.5plog(det(YY))-0.5(delta+p+1)log(det(XX))-0.5trace(YYinv(XX));
	1382	cout << log1 << "," << log2 << endl;
	1383	}*/
	1384	return log1;
	1385	};
1266	1386
1267	1387	};
1268	1388
1269	1389	//! Random Walk on inverse Wishart
1270		class rwiWishartCh : public pdf_internal<eiWishartCh>
1271		{
1272		protected:
1273		//!square root of \f$ \nu-p-1 \f$ - needed for computation of \f$ \Psi \f$ from conditions
1274		double sqd;
1275		//!reference point for diagonal
1276		vec refl;
1277		//! power of the reference
1278		double l;
1279		//! dimension
1280		int p;
1281
1282		public:
1283		rwiWishartCh() : sqd (0), l (0), p (0) {}
1284		//! constructor function
1285		void set_parameters (int p0, double k, vec ref0, double l0) {
1286		p = p0;
1287		double delta = 2 / (k * k) + p + 3;
1288		sqd = sqrt (delta - p - 1);
1289		l = l0;
1290		refl = pow (ref0, 1 - l);
1291
1292		iepdf.set_parameters (eye (p), delta);
1293		dimc = iepdf.dimension();
1294		}
1295		void condition (const vec &c) {
1296		vec z = c;
1297		int ri = 0;
1298		for (int i = 0;i < p*p;i += (p + 1)) {//trace diagonal element
1299		z (i) = pow (z (i), l) * refl (ri);
1300		ri++;
1301		}
1302
1303		iepdf._setY (sqd*z);
1304		}
	1390	class rwiWishartCh : public pdf_internal<eiWishartCh> {
	1391	protected:
	1392	//!square root of \f$ \nu-p-1 \f$ - needed for computation of \f$ \Psi \f$ from conditions
	1393	double sqd;
	1394	//!reference point for diagonal
	1395	vec refl;
	1396	//! power of the reference
	1397	double l;
	1398	//! dimension
	1399	int p;
	1400
	1401	public:
	1402	rwiWishartCh() : sqd ( 0 ), l ( 0 ), p ( 0 ) {}
	1403	//! constructor function
	1404	void set_parameters ( int p0, double k, vec ref0, double l0 ) {
	1405	p = p0;
	1406	double delta = 2 / ( k * k ) + p + 3;
	1407	sqd = sqrt ( delta - p - 1 );
	1408	l = l0;
	1409	refl = pow ( ref0, 1 - l );
	1410
	1411	iepdf.set_parameters ( eye ( p ), delta );
	1412	dimc = iepdf.dimension();
	1413	}
	1414	void condition ( const vec &c ) {
	1415	vec z = c;
	1416	int ri = 0;
	1417	for ( int i = 0; i < p*p; i += ( p + 1 ) ) {//trace diagonal element
	1418	z ( i ) = pow ( z ( i ), l ) * refl ( ri );
	1419	ri++;
	1420	}
	1421
	1422	iepdf._setY ( sqd*z );
	1423	}
1305	1424	};
1306	1425
…	…
1312	1431	Used e.g. in particle filters.
1313	1432	*/
1314		class eEmp: public epdf
1315		{
1316		protected :
1317		//! Number of particles
1318		int n;
1319		//! Sample weights \f$w\f$
1320		vec w;
1321		//! Samples \f$x^{(i)}, i=1..n\f$
1322		Array<vec> samples;
1323		public:
1324		//! \name Constructors
1325		//!@{
1326		eEmp () : epdf (), w (), samples () {};
1327		//! copy constructor
1328		eEmp (const eEmp &e) : epdf (e), w (e.w), samples (e.samples) {};
1329		//!@}
1330
1331		//! Set samples and weights
1332		void set_statistics (const vec &w0, const epdf &pdf0);
1333		//! Set samples and weights
1334		void set_statistics (const epdf &pdf0 , int n) {set_statistics (ones (n) / n, pdf0);};
1335		//! Set sample
1336		void set_samples (const epdf* pdf0);
1337		//! Set sample
1338		void set_parameters (int n0, bool copy = true) {n = n0; w.set_size (n0, copy);samples.set_size (n0, copy);};
1339		//! Set samples
1340		void set_parameters (const Array<vec> &Av) {
1341		bdm_assert(Av.size()>0,"Empty samples");
1342		n = Av.size();
1343		epdf::set_parameters(Av(0).length());
1344		w=1/n*ones(n);
1345		samples=Av;
1346		};
1347		//! Potentially dangerous, use with care.
1348		vec& _w() {return w;};
1349		//! Potentially dangerous, use with care.
1350		const vec& _w() const {return w;};
1351		//! access function
1352		Array<vec>& _samples() {return samples;};
1353		//! access function
1354		const vec& _sample(int i) const {return samples(i);};
1355		//! access function
1356		const Array<vec>& _samples() const {return samples;};
1357		//! 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.
1358		//! The vector with indeces of new samples is returned in variable \c index.
1359		void resample ( ivec &index, RESAMPLING_METHOD method = SYSTEMATIC);
1360
1361		//! Resampling without returning index of new particles.
1362		void resample (RESAMPLING_METHOD method = SYSTEMATIC){ivec ind; resample(ind,method);};
1363
1364		//! inherited operation : NOT implemented
1365		vec sample() const {
1366		bdm_error ("Not implemented");
1367		return vec();
1368		}
1369
1370		//! inherited operation : NOT implemented
1371		double evallog (const vec &val) const {
1372		bdm_error ("Not implemented");
1373		return 0.0;
1374		}
1375
1376		vec mean() const {
1377		vec pom = zeros (dim);
1378		for (int i = 0;i < n;i++) {pom += samples (i) * w (i);}
1379		return pom;
1380		}
1381		vec variance() const {
1382		vec pom = zeros (dim);
1383		for (int i = 0;i < n;i++) {pom += pow (samples (i), 2) * w (i);}
1384		return pom -pow (mean(), 2);
1385		}
1386		//! For this class, qbounds are minimum and maximum value of the population!
1387		void qbounds (vec &lb, vec &ub, double perc = 0.95) const {
1388		// lb in inf so than it will be pushed below;
1389		lb.set_size (dim);
1390		ub.set_size (dim);
1391		lb = std::numeric_limits<double>::infinity();
1392		ub = -std::numeric_limits<double>::infinity();
1393		int j;
1394		for (int i = 0;i < n;i++) {
1395		for (j = 0;j < dim; j++) {
1396		if (samples (i) (j) < lb (j)) {lb (j) = samples (i) (j);}
1397		if (samples (i) (j) > ub (j)) {ub (j) = samples (i) (j);}
	1433	class eEmp: public epdf {
	1434	protected :
	1435	//! Number of particles
	1436	int n;
	1437	//! Sample weights \f$w\f$
	1438	vec w;
	1439	//! Samples \f$x^{(i)}, i=1..n\f$
	1440	Array<vec> samples;
	1441	public:
	1442	//! \name Constructors
	1443	//!@{
	1444	eEmp () : epdf (), w (), samples () {};
	1445	//! copy constructor
	1446	eEmp ( const eEmp &e ) : epdf ( e ), w ( e.w ), samples ( e.samples ) {};
	1447	//!@}
	1448
	1449	//! Set samples and weights
	1450	void set_statistics ( const vec &w0, const epdf &pdf0 );
	1451	//! Set samples and weights
	1452	void set_statistics ( const epdf &pdf0 , int n ) {
	1453	set_statistics ( ones ( n ) / n, pdf0 );
	1454	};
	1455	//! Set sample
	1456	void set_samples ( const epdf* pdf0 );
	1457	//! Set sample
	1458	void set_parameters ( int n0, bool copy = true ) {
	1459	n = n0;
	1460	w.set_size ( n0, copy );
	1461	samples.set_size ( n0, copy );
	1462	};
	1463	//! Set samples
	1464	void set_parameters ( const Array<vec> &Av ) {
	1465	bdm_assert ( Av.size() > 0, "Empty samples" );
	1466	n = Av.size();
	1467	epdf::set_parameters ( Av ( 0 ).length() );
	1468	w = 1 / n * ones ( n );
	1469	samples = Av;
	1470	};
	1471	//! Potentially dangerous, use with care.
	1472	vec& _w() {
	1473	return w;
	1474	};
	1475	//! Potentially dangerous, use with care.
	1476	const vec& _w() const {
	1477	return w;
	1478	};
	1479	//! access function
	1480	Array<vec>& _samples() {
	1481	return samples;
	1482	};
	1483	//! access function
	1484	const vec& _sample ( int i ) const {
	1485	return samples ( i );
	1486	};
	1487	//! access function
	1488	const Array<vec>& _samples() const {
	1489	return samples;
	1490	};
	1491	//! 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.
	1492	//! The vector with indeces of new samples is returned in variable \c index.
	1493	void resample ( ivec &index, RESAMPLING_METHOD method = SYSTEMATIC );
	1494
	1495	//! Resampling without returning index of new particles.
	1496	void resample ( RESAMPLING_METHOD method = SYSTEMATIC ) {
	1497	ivec ind;
	1498	resample ( ind, method );
	1499	};
	1500
	1501	//! inherited operation : NOT implemented
	1502	vec sample() const {
	1503	bdm_error ( "Not implemented" );
	1504	return vec();
	1505	}
	1506
	1507	//! inherited operation : NOT implemented
	1508	double evallog ( const vec &val ) const {
	1509	bdm_error ( "Not implemented" );
	1510	return 0.0;
	1511	}
	1512
	1513	vec mean() const {
	1514	vec pom = zeros ( dim );
	1515	for ( int i = 0; i < n; i++ ) {
	1516	pom += samples ( i ) * w ( i );
	1517	}
	1518	return pom;
	1519	}
	1520	vec variance() const {
	1521	vec pom = zeros ( dim );
	1522	for ( int i = 0; i < n; i++ ) {
	1523	pom += pow ( samples ( i ), 2 ) * w ( i );
	1524	}
	1525	return pom - pow ( mean(), 2 );
	1526	}
	1527	//! For this class, qbounds are minimum and maximum value of the population!
	1528	void qbounds ( vec &lb, vec &ub, double perc = 0.95 ) const {
	1529	// lb in inf so than it will be pushed below;
	1530	lb.set_size ( dim );
	1531	ub.set_size ( dim );
	1532	lb = std::numeric_limits<double>::infinity();
	1533	ub = -std::numeric_limits<double>::infinity();
	1534	int j;
	1535	for ( int i = 0; i < n; i++ ) {
	1536	for ( j = 0; j < dim; j++ ) {
	1537	if ( samples ( i ) ( j ) < lb ( j ) ) {
	1538	lb ( j ) = samples ( i ) ( j );
	1539	}
	1540	if ( samples ( i ) ( j ) > ub ( j ) ) {
	1541	ub ( j ) = samples ( i ) ( j );
1398	1542	}
1399	1543	}
1400	1544	}
	1545	}
1401	1546	};
1402	1547
…	…
1405	1550
1406	1551	template<class sq_T>
1407		void enorm<sq_T>::set_parameters (const vec &mu0, const sq_T &R0)
1408		{
	1552	void enorm<sq_T>::set_parameters ( const vec &mu0, const sq_T &R0 ) {
1409	1553	//Fixme test dimensions of mu0 and R0;
1410	1554	mu = mu0;
…	…
1414	1558
1415	1559	template<class sq_T>
1416		void enorm<sq_T>::from_setting (const Setting &set)
1417		{
1418		epdf::from_setting (set); //reads rv
1419
1420		UI::get (mu, set, "mu", UI::compulsory);
	1560	void enorm<sq_T>::from_setting ( const Setting &set ) {
	1561	epdf::from_setting ( set ); //reads rv
	1562
	1563	UI::get ( mu, set, "mu", UI::compulsory );
1421	1564	mat Rtmp;// necessary for conversion
1422		UI::get (~~Rtmp, set, "R", UI::compulsory~~);
	1565	UI::get ( Rtmp, set, "R", UI::compulsory );
1423	1566	R = Rtmp; // conversion
1424	1567	validate();
…	…
1426	1569
1427	1570	template<class sq_T>
1428		void enorm<sq_T>::dupdate (mat &v, double nu)
1429		{
	1571	void enorm<sq_T>::dupdate ( mat &v, double nu ) {
1430	1572	//
1431	1573	};
…	…
1437	1579
1438	1580	template<class sq_T>
1439		vec enorm<sq_T>::sample() const
1440		{
1441		vec x (dim);
	1581	vec enorm<sq_T>::sample() const {
	1582	vec x ( dim );
1442	1583	#pragma omp critical
1443		NorRNG.sample_vector (~~dim, x~~);
1444		vec smp = R.sqrt_mult (x);
	1584	NorRNG.sample_vector ( dim, x );
	1585	vec smp = R.sqrt_mult ( x );
1445	1586
1446	1587	smp += mu;
…	…
1457	1598
1458	1599	template<class sq_T>
1459		double enorm<sq_T>::evallog_nn (const vec &val) const
1460		{
	1600	double enorm<sq_T>::evallog_nn ( const vec &val ) const {
1461	1601	// 1.83787706640935 = log(2pi)
1462		double tmp = -0.5 * (~~R.invqform (mu - val)~~);// - lognc();
	1602	double tmp = -0.5 * ( R.invqform ( mu - val ) );// - lognc();
1463	1603	return tmp;
1464	1604	};
1465	1605
1466	1606	template<class sq_T>
1467		inline double enorm<sq_T>::lognc () const
1468		{
	1607	inline double enorm<sq_T>::lognc () const {
1469	1608	// 1.83787706640935 = log(2pi)
1470		double tmp = 0.5 * (~~R.cols() * 1.83787706640935 + R.logdet()~~);
	1609	double tmp = 0.5 * ( R.cols() * 1.83787706640935 + R.logdet() );
1471	1610	return tmp;
1472	1611	};
…	…
1500	1639
1501	1640	template<class sq_T>
1502		shared_ptr<epdf> enorm<sq_T>::marginal ( const RV &rvn ) const
1503		{
	1641	shared_ptr<epdf> enorm<sq_T>::marginal ( const RV &rvn ) const {
1504	1642	enorm<sq_T> *tmp = new enorm<sq_T> ();
1505		shared_ptr<epdf> narrow~~(tmp~~);
	1643	shared_ptr<epdf> narrow ( tmp );
1506	1644	marginal ( rvn, *tmp );
1507	1645	return narrow;
…	…
1509	1647
1510	1648	template<class sq_T>
1511		void enorm<sq_T>::marginal ( const RV &rvn, enorm<sq_T> &target ) const
1512		{
1513		bdm_assert (isnamed(), "rv description is not assigned");
1514		ivec irvn = rvn.dataind (rv);
1515
1516		sq_T Rn (R, irvn); // select rows and columns of R
	1649	void enorm<sq_T>::marginal ( const RV &rvn, enorm<sq_T> &target ) const {
	1650	bdm_assert ( isnamed(), "rv description is not assigned" );
	1651	ivec irvn = rvn.dataind ( rv );
	1652
	1653	sq_T Rn ( R, irvn ); // select rows and columns of R
1517	1654
1518	1655	target.set_rv ( rvn );
1519		target.set_parameters (~~mu (irvn), Rn~~);
	1656	target.set_parameters ( mu ( irvn ), Rn );
1520	1657	}
1521	1658
1522	1659	template<class sq_T>
1523		shared_ptr<pdf> enorm<sq_T>::condition ( const RV &rvn ) const
1524		{
	1660	shared_ptr<pdf> enorm<sq_T>::condition ( const RV &rvn ) const {
1525	1661	mlnorm<sq_T> *tmp = new mlnorm<sq_T> ();
1526		shared_ptr<pdf> narrow~~(tmp~~);
	1662	shared_ptr<pdf> narrow ( tmp );
1527	1663	condition ( rvn, *tmp );
1528	1664	return narrow;
…	…
1530	1666
1531	1667	template<class sq_T>
1532		void enorm<sq_T>::condition ( const RV &rvn, pdf &target ) const
1533		{
	1668	void enorm<sq_T>::condition ( const RV &rvn, pdf &target ) const {
1534	1669	typedef mlnorm<sq_T> TMlnorm;
1535	1670
1536		bdm_assert (~~isnamed(), "rvs are not assigned"~~);
1537		TMlnorm &uptarget = dynamic_cast<TMlnorm &>~~(target~~);
1538
1539		RV rvc = rv.subt (~~rvn~~);
1540		bdm_assert ( (~~rvc._dsize() + rvn._dsize() == rv._dsize()), "wrong rvn"~~);
	1671	bdm_assert ( isnamed(), "rvs are not assigned" );
	1672	TMlnorm &uptarget = dynamic_cast<TMlnorm &> ( target );
	1673
	1674	RV rvc = rv.subt ( rvn );
	1675	bdm_assert ( ( rvc._dsize() + rvn._dsize() == rv._dsize() ), "wrong rvn" );
1541	1676	//Permutation vector of the new R
1542		ivec irvn = rvn.dataind (rv);
1543		ivec irvc = rvc.dataind (rv);
1544		ivec perm = concat (~~irvn , irvc~~);
1545		sq_T Rn (~~R, perm~~);
	1677	ivec irvn = rvn.dataind ( rv );
	1678	ivec irvc = rvc.dataind ( rv );
	1679	ivec perm = concat ( irvn , irvc );
	1680	sq_T Rn ( R, perm );
1546	1681
1547	1682	//fixme - could this be done in general for all sq_T?
…	…
1550	1685	int n = rvn._dsize() - 1;
1551	1686	int end = R.rows() - 1;
1552		mat S11 = S.get (~~0, n, 0, n~~);
1553		mat S12 = S.get (~~0, n , rvn._dsize(), end~~);
1554		mat S22 = S.get (~~rvn._dsize(), end, rvn._dsize(), end~~);
1555
1556		vec mu1 = mu (~~irvn~~);
1557		vec mu2 = mu (~~irvc~~);
1558		mat A = S12 * inv (~~S22~~);
1559		sq_T R_n (~~S11 - A *S12.T()~~);
1560
1561		uptarget.set_rv (~~rvn~~);
1562		uptarget.set_rvc (~~rvc~~);
1563		uptarget.set_parameters (~~A, mu1 - A*mu2, R_n~~);
	1687	mat S11 = S.get ( 0, n, 0, n );
	1688	mat S12 = S.get ( 0, n , rvn._dsize(), end );
	1689	mat S22 = S.get ( rvn._dsize(), end, rvn._dsize(), end );
	1690
	1691	vec mu1 = mu ( irvn );
	1692	vec mu2 = mu ( irvc );
	1693	mat A = S12 * inv ( S22 );
	1694	sq_T R_n ( S11 - A *S12.T() );
	1695
	1696	uptarget.set_rv ( rvn );
	1697	uptarget.set_rvc ( rvc );
	1698	uptarget.set_parameters ( A, mu1 - A*mu2, R_n );
1564	1699	}
1565	1700
…	…
1567	1702	///////
1568	1703	template<class sq_T>
1569		void mgnorm<sq_T >::set_parameters (~~const shared_ptr<fnc> &g0, const sq_T &R0~~) {
	1704	void mgnorm<sq_T >::set_parameters ( const shared_ptr<fnc> &g0, const sq_T &R0 ) {
1570	1705	g = g0;
1571		this->iepdf.set_parameters (~~zeros (g->dimension()), R0~~);
	1706	this->iepdf.set_parameters ( zeros ( g->dimension() ), R0 );
1572	1707	}
1573	1708
1574	1709	template<class sq_T>
1575		void mgnorm<sq_T >::condition (const vec &cond) {this->iepdf._mu() = g->eval (cond);};
	1710	void mgnorm<sq_T >::condition ( const vec &cond ) {
	1711	this->iepdf._mu() = g->eval ( cond );
	1712	};
1576	1713
1577	1714	//! \todo unify this stuff with to_string()
1578	1715	template<class sq_T>
1579		std::ostream &operator<< (std::ostream &os, mlnorm<sq_T> &ml)
1580		{
	1716	std::ostream &operator<< ( std::ostream &os, mlnorm<sq_T> &ml ) {
1581	1717	os << "A:" << ml.A << endl;
1582	1718	os << "mu:" << ml.mu_const << endl;

Context Navigation

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

Legend:

library/bdm/stat/exp_family.h

Download in other formats: