Context Navigation

← Previous Change
Next Change →

exp_family.h

Timestamp:

08/08/09 13:43:13 (15 years ago)

Author:

smidl

Message:

changes in mpdf -> compile OK, broken tests!

Files:

: 1 modified

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

Legend:

: Unmodified
: Added
: Removed

library/bdm/stat/exp_family.h

r487	r488
24	24
25	25	//! Global Uniform_RNG
26		extern Uniform_RNG UniRNG;
	26	extern Uniform_RNG UniRNG;
27	27	//! Global Normal_RNG
28		extern Normal_RNG NorRNG;
	28	extern Normal_RNG NorRNG;
29	29	//! Global Gamma_RNG
30		extern Gamma_RNG GamRNG;
31
32		/*!
33		* \brief General conjugate exponential family posterior density.
34
35		* More?...
36		*/
37
38		class eEF : public epdf
39		{
40		public:
	30	extern Gamma_RNG GamRNG;
	31
	32	/*!
	33	* \brief General conjugate exponential family posterior density.
	34
	35	* More?...
	36	*/
	37
	38	class eEF : public epdf
	39	{
	40	public:
41	41	// 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		//!TODO decide if it is really needed
47		virtual void dupdate ( mat &v ) {it_error ( "Not implemented" );};
48		//!Evaluate normalized log-probability
49		virtual double evallog_nn ( const vec &val ) const{it_error ( "Not implemented" );return 0.0;};
50		//!Evaluate normalized log-probability
51		virtual double evallog ( const vec &val ) const {
52		double tmp;
53		tmp= evallog_nn ( val )-lognc();
54		// it_assert_debug ( std::isfinite ( tmp ),"Infinite value" );
55		return tmp;}
56		//!Evaluate normalized log-probability for many samples
57		virtual vec evallog_m ( const mat &Val ) const
58		{
59		vec x ( Val.cols() );
60		for ( int i=0;i<Val.cols();i++ ) {x ( i ) =evallog_nn ( Val.get_col ( i ) ) ;}
61		return x-lognc();
	42	//! default constructor
	43	eEF () : epdf () {};
	44	//! logarithm of the normalizing constant, \f$\mathcal{I}\f$
	45	virtual double lognc() const = 0;
	46	//!TODO decide if it is really needed
	47	virtual void dupdate (mat &v) {it_error ("Not implemented");};
	48	//!Evaluate normalized log-probability
	49	virtual double evallog_nn (const vec &val) const{it_error ("Not implemented");return 0.0;};
	50	//!Evaluate normalized log-probability
	51	virtual double evallog (const vec &val) const {
	52	double tmp;
	53	tmp = evallog_nn (val) - lognc();
	54	// it_assert_debug ( std::isfinite ( tmp ),"Infinite value" );
	55	return tmp;
	56	}
	57	//!Evaluate normalized log-probability for many samples
	58	virtual vec evallog_m (const mat &Val) const {
	59	vec x (Val.cols());
	60	for (int i = 0;i < Val.cols();i++) {x (i) = evallog_nn (Val.get_col (i)) ;}
	61	return x -lognc();
	62	}
	63	//!Evaluate normalized log-probability for many samples
	64	virtual vec evallog_m (const Array<vec> &Val) const {
	65	vec x (Val.length());
	66	for (int i = 0;i < Val.length();i++) {x (i) = evallog_nn (Val (i)) ;}
	67	return x -lognc();
	68	}
	69	//!Power of the density, used e.g. to flatten the density
	70	virtual void pow (double p) {it_error ("Not implemented");};
	71	};
	72
	73
	74	//! Estimator for Exponential family
	75	class BMEF : public BM
	76	{
	77	protected:
	78	//! forgetting factor
	79	double frg;
	80	//! cached value of lognc() in the previous step (used in evaluation of \c ll )
	81	double last_lognc;
	82	public:
	83	//! Default constructor (=empty constructor)
	84	BMEF (double frg0 = 1.0) : BM (), frg (frg0) {}
	85	//! Copy constructor
	86	BMEF (const BMEF &B) : BM (B), frg (B.frg), last_lognc (B.last_lognc) {}
	87	//!get statistics from another model
	88	virtual void set_statistics (const BMEF* BM0) {it_error ("Not implemented");};
	89	//! Weighted update of sufficient statistics (Bayes rule)
	90	virtual void bayes (const vec &data, const double w) {};
	91	//original Bayes
	92	void bayes (const vec &dt);
	93	//!Flatten the posterior according to the given BMEF (of the same type!)
	94	virtual void flatten (const BMEF * B) {it_error ("Not implemented");}
	95	//!Flatten the posterior as if to keep nu0 data
	96	// virtual void flatten ( double nu0 ) {it_error ( "Not implemented" );}
	97
	98	BMEF* _copy_ () const {it_error ("function _copy_ not implemented for this BM"); return NULL;};
	99	};
	100
	101	template<class sq_T, template <typename> class TEpdf >
	102	class mlnorm;
	103
	104	/*!
	105	* \brief Gaussian density with positive definite (decomposed) covariance matrix.
	106
	107	* More?...
	108	*/
	109	template<class sq_T>
	110	class enorm : public eEF
	111	{
	112	protected:
	113	//! mean value
	114	vec mu;
	115	//! Covariance matrix in decomposed form
	116	sq_T R;
	117	public:
	118	//!\name Constructors
	119	//!@{
	120
	121	enorm () : eEF (), mu (), R () {};
	122	enorm (const vec &mu, const sq_T &R) {set_parameters (mu, R);}
	123	void set_parameters (const vec &mu, const sq_T &R);
	124	void from_setting (const Setting &root);
	125	void validate() {
	126	it_assert (mu.length() == R.rows(), "parameters mismatch");
	127	dim = mu.length();
	128	}
	129	//!@}
	130
	131	//! \name Mathematical operations
	132	//!@{
	133
	134	//! dupdate in exponential form (not really handy)
	135	void dupdate (mat &v, double nu = 1.0);
	136
	137	vec sample() const;
	138
	139	double evallog_nn (const vec &val) const;
	140	double lognc () const;
	141	vec mean() const {return mu;}
	142	vec variance() const {return diag (R.to_mat());}
	143	// mlnorm<sq_T>* condition ( const RV &rvn ) const ; <=========== fails to cmpile. Why?
	144	mpdf* condition (const RV &rvn) const ;
	145	enorm<sq_T>* marginal (const RV &rv) const;
	146	// epdf* marginal ( const RV &rv ) const;
	147	//!@}
	148
	149	//! \name Access to attributes
	150	//!@{
	151
	152	vec& _mu() {return mu;}
	153	void set_mu (const vec mu0) { mu = mu0;}
	154	sq_T& _R() {return R;}
	155	const sq_T& _R() const {return R;}
	156	//!@}
	157
	158	};
	159	UIREGISTER (enorm<chmat>);
	160	UIREGISTER (enorm<ldmat>);
	161	UIREGISTER (enorm<fsqmat>);
	162
	163
	164	/*!
	165	* \brief Gauss-inverse-Wishart density stored in LD form
	166
	167	* For \f$p\f$-variate densities, given rv.count() should be \f$p\times\f$ V.rows().
	168	*
	169	*/
	170	class egiw : public eEF
	171	{
	172	protected:
	173	//! Extended information matrix of sufficient statistics
	174	ldmat V;
	175	//! Number of data records (degrees of freedom) of sufficient statistics
	176	double nu;
	177	//! Dimension of the output
	178	int dimx;
	179	//! Dimension of the regressor
	180	int nPsi;
	181	public:
	182	//!\name Constructors
	183	//!@{
	184	egiw() : eEF() {};
	185	egiw (int dimx0, ldmat V0, double nu0 = -1.0) : eEF() {set_parameters (dimx0, V0, nu0);};
	186
	187	void set_parameters (int dimx0, ldmat V0, double nu0 = -1.0) {
	188	dimx = dimx0;
	189	nPsi = V0.rows() - dimx;
	190	dim = dimx * (dimx + nPsi); // size(R) + size(Theta)
	191
	192	V = V0;
	193	if (nu0 < 0) {
	194	nu = 0.1 + nPsi + 2 * dimx + 2; // +2 assures finite expected value of R
	195	// terms before that are sufficient for finite normalization
	196	} else {
	197	nu = nu0;
62	198	}
63		//!Power of the density, used e.g. to flatten the density
64		virtual void pow ( double p ) {it_error ( "Not implemented" );};
65		};
66
67
68		//! Estimator for Exponential family
69		class BMEF : public BM
70		{
71		protected:
72		//! forgetting factor
73		double frg;
74		//! cached value of lognc() in the previous step (used in evaluation of \c ll )
75		double last_lognc;
76		public:
77		//! Default constructor (=empty constructor)
78		BMEF ( double frg0=1.0 ) :BM ( ), frg ( frg0 ) {}
79		//! Copy constructor
80		BMEF ( const BMEF &B ) :BM ( B ), frg ( B.frg ), last_lognc ( B.last_lognc ) {}
81		//!get statistics from another model
82		virtual void set_statistics ( const BMEF* BM0 ) {it_error ( "Not implemented" );};
83		//! Weighted update of sufficient statistics (Bayes rule)
84		virtual void bayes ( const vec &data, const double w ) {};
85		//original Bayes
86		void bayes ( const vec &dt );
87		//!Flatten the posterior according to the given BMEF (of the same type!)
88		virtual void flatten ( const BMEF * B ) {it_error ( "Not implemented" );}
89		//!Flatten the posterior as if to keep nu0 data
90		// virtual void flatten ( double nu0 ) {it_error ( "Not implemented" );}
91
92		BMEF* _copy_ () const {it_error ( "function _copy_ not implemented for this BM" ); return NULL;};
93		};
94
95		template<class sq_T, template <typename> class TEpdf >
96		class mlnorm;
97
98		/*!
99		* \brief Gaussian density with positive definite (decomposed) covariance matrix.
100
101		* More?...
102		*/
103		template<class sq_T>
104		class enorm : public eEF
105		{
106		protected:
107		//! mean value
108		vec mu;
109		//! Covariance matrix in decomposed form
110		sq_T R;
111		public:
112		//!\name Constructors
113		//!@{
114
115		enorm ( ) :eEF ( ), mu ( ),R ( ) {};
116		enorm ( const vec &mu,const sq_T &R ) {set_parameters ( mu,R );}
117		void set_parameters ( const vec &mu,const sq_T &R );
118		void from_setting(const Setting &root);
119		void validate() {
120		it_assert(mu.length()==R.rows(),"parameters mismatch");
121		dim = mu.length();
	199	}
	200	//!@}
	201
	202	vec sample() const;
	203	vec mean() const;
	204	vec variance() const;
	205
	206	//! LS estimate of \f$\theta\f$
	207	vec est_theta() const;
	208
	209	//! Covariance of the LS estimate
	210	ldmat est_theta_cov() const;
	211
	212	void mean_mat (mat &M, mat&R) const;
	213	//! In this instance, val= [theta, r]. For multivariate instances, it is stored columnwise val = [theta_1 theta_2 ... r_1 r_2 ]
	214	double evallog_nn (const vec &val) const;
	215	double lognc () const;
	216	void pow (double p) {V = p;nu = p;};
	217
	218	//! \name Access attributes
	219	//!@{
	220
	221	ldmat& _V() {return V;}
	222	const ldmat& _V() const {return V;}
	223	double& _nu() {return nu;}
	224	const double& _nu() const {return nu;}
	225	void from_setting (const Setting &set) {
	226	UI::get (nu, set, "nu", UI::compulsory);
	227	UI::get (dimx, set, "dimx", UI::compulsory);
	228	mat V;
	229	UI::get (V, set, "V", UI::compulsory);
	230	set_parameters (dimx, V, nu);
	231	RV* rv = UI::build<RV> (set, "rv", UI::compulsory);
	232	set_rv (*rv);
	233	delete rv;
	234	}
	235	//!@}
	236	};
	237	UIREGISTER (egiw);
	238
	239	/*! \brief Dirichlet posterior density
	240
	241	Continuous Dirichlet density of \f$n\f$-dimensional variable \f$x\f$
	242	\f[
	243	f(x\|\beta) = \frac{\Gamma[\gamma]}{\prod_{i=1}^{n}\Gamma(\beta_i)} \prod_{i=1}^{n}x_i^{\beta_i-1}
	244	\f]
	245	where \f$\gamma=\sum_i \beta_i\f$.
	246	*/
	247	class eDirich: public eEF
	248	{
	249	protected:
	250	//!sufficient statistics
	251	vec beta;
	252	public:
	253	//!\name Constructors
	254	//!@{
	255
	256	eDirich () : eEF () {};
	257	eDirich (const eDirich &D0) : eEF () {set_parameters (D0.beta);};
	258	eDirich (const vec &beta0) {set_parameters (beta0);};
	259	void set_parameters (const vec &beta0) {
	260	beta = beta0;
	261	dim = beta.length();
	262	}
	263	//!@}
	264
	265	vec sample() const {it_error ("Not implemented");return vec_1 (0.0);};
	266	vec mean() const {return beta / sum (beta);};
	267	vec variance() const {double gamma = sum (beta); return elem_mult (beta, (beta + 1)) / (gamma* (gamma + 1));}
	268	//! In this instance, val is ...
	269	double evallog_nn (const vec &val) const {
	270	double tmp; tmp = (beta - 1) * log (val);
	271	// it_assert_debug ( std::isfinite ( tmp ),"Infinite value" );
	272	return tmp;
	273	};
	274	double lognc () const {
	275	double tmp;
	276	double gam = sum (beta);
	277	double lgb = 0.0;
	278	for (int i = 0;i < beta.length();i++) {lgb += lgamma (beta (i));}
	279	tmp = lgb - lgamma (gam);
	280	// it_assert_debug ( std::isfinite ( tmp ),"Infinite value" );
	281	return tmp;
	282	};
	283	//!access function
	284	vec& _beta() {return beta;}
	285	//!Set internal parameters
	286	};
	287
	288	//! \brief Estimator for Multinomial density
	289	class multiBM : public BMEF
	290	{
	291	protected:
	292	//! Conjugate prior and posterior
	293	eDirich est;
	294	//! Pointer inside est to sufficient statistics
	295	vec β
	296	public:
	297	//!Default constructor
	298	multiBM () : BMEF (), est (), beta (est._beta()) {
	299	if (beta.length() > 0) {last_lognc = est.lognc();}
	300	else{last_lognc = 0.0;}
	301	}
	302	//!Copy constructor
	303	multiBM (const multiBM &B) : BMEF (B), est (B.est), beta (est._beta()) {}
	304	//! Sets sufficient statistics to match that of givefrom mB0
	305	void set_statistics (const BM* mB0) {const multiBM* mB = dynamic_cast<const multiBM*> (mB0); beta = mB->beta;}
	306	void bayes (const vec &dt) {
	307	if (frg < 1.0) {beta *= frg;last_lognc = est.lognc();}
	308	beta += dt;
	309	if (evalll) {ll = est.lognc() - last_lognc;}
	310	}
	311	double logpred (const vec &dt) const {
	312	eDirich pred (est);
	313	vec &beta = pred._beta();
	314
	315	double lll;
	316	if (frg < 1.0)
	317	{beta *= frg;lll = pred.lognc();}
	318	else
	319	if (evalll) {lll = last_lognc;}
	320	else{lll = pred.lognc();}
	321
	322	beta += dt;
	323	return pred.lognc() - lll;
	324	}
	325	void flatten (const BMEF* B) {
	326	const multiBM* E = dynamic_cast<const multiBM*> (B);
	327	// sum(beta) should be equal to sum(B.beta)
	328	const vec &Eb = E->beta;//const_cast<multiBM*> ( E )->_beta();
	329	beta *= (sum (Eb) / sum (beta));
	330	if (evalll) {last_lognc = est.lognc();}
	331	}
	332	const epdf& posterior() const {return est;};
	333	const eDirich* _e() const {return &est;};
	334	void set_parameters (const vec &beta0) {
	335	est.set_parameters (beta0);
	336	if (evalll) {last_lognc = est.lognc();}
	337	}
	338	};
	339
	340	/*!
	341	\brief Gamma posterior density
	342
	343	Multivariate Gamma density as product of independent univariate densities.
	344	\f[
	345	f(x\|\alpha,\beta) = \prod f(x_i\|\alpha_i,\beta_i)
	346	\f]
	347	*/
	348
	349	class egamma : public eEF
	350	{
	351	protected:
	352	//! Vector \f$\alpha\f$
	353	vec alpha;
	354	//! Vector \f$\beta\f$
	355	vec beta;
	356	public :
	357	//! \name Constructors
	358	//!@{
	359	egamma () : eEF (), alpha (0), beta (0) {};
	360	egamma (const vec &a, const vec &b) {set_parameters (a, b);};
	361	void set_parameters (const vec &a, const vec &b) {alpha = a, beta = b;dim = alpha.length();};
	362	//!@}
	363
	364	vec sample() const;
	365	//! TODO: is it used anywhere?
	366	// mat sample ( int N ) const;
	367	double evallog (const vec &val) const;
	368	double lognc () const;
	369	//! Returns poiter to alpha and beta. Potentially dengerous: use with care!
	370	vec& _alpha() {return alpha;}
	371	vec& _beta() {return beta;}
	372	vec mean() const {return elem_div (alpha, beta);}
	373	vec variance() const {return elem_div (alpha, elem_mult (beta, beta)); }
	374
	375	//! Load from structure with elements:
	376	//! \code
	377	//! { alpha = [...]; // vector of alpha
	378	//! beta = [...]; // vector of beta
	379	//! rv = {class="RV",...} // description
	380	//! }
	381	//! \endcode
	382	//!@}
	383	void from_setting (const Setting &set) {
	384	epdf::from_setting (set); // reads rv
	385	UI::get (alpha, set, "alpha", UI::compulsory);
	386	UI::get (beta, set, "beta", UI::compulsory);
	387	validate();
	388	}
	389	void validate() {
	390	it_assert (alpha.length() == beta.length(), "parameters do not match");
	391	dim = alpha.length();
	392	}
	393	};
	394	UIREGISTER (egamma);
	395	/*!
	396	\brief Inverse-Gamma posterior density
	397
	398	Multivariate inverse-Gamma density as product of independent univariate densities.
	399	\f[
	400	f(x\|\alpha,\beta) = \prod f(x_i\|\alpha_i,\beta_i)
	401	\f]
	402
	403	Vector \f$\beta\f$ has different meaning (in fact it is 1/beta as used in definition of iG)
	404
	405	Inverse Gamma can be converted to Gamma using \f[
	406	x\sim iG(a,b) => 1/x\sim G(a,1/b)
	407	\f]
	408	This relation is used in sampling.
	409	*/
	410
	411	class eigamma : public egamma
	412	{
	413	protected:
	414	public :
	415	//! \name Constructors
	416	//! All constructors are inherited
	417	//!@{
	418	//!@}
	419
	420	vec sample() const {return 1.0 / egamma::sample();};
	421	//! Returns poiter to alpha and beta. Potentially dangerous: use with care!
	422	vec mean() const {return elem_div (beta, alpha - 1);}
	423	vec variance() const {vec mea = mean(); return elem_div (elem_mult (mea, mea), alpha - 2);}
	424	};
	425	/*
	426	//! Weighted mixture of epdfs with external owned components.
	427	class emix : public epdf {
	428	protected:
	429	int n;
	430	vec &w;
	431	Array<epdf*> Coms;
	432	public:
	433	//! Default constructor
	434	emix ( const RV &rv, vec &w0): epdf(rv), n(w0.length()), w(w0), Coms(n) {};
	435	void set_parameters( int &i, double wi, epdf* ep){w(i)=wi;Coms(i)=ep;}
	436	vec mean(){vec pom; for(int i=0;i<n;i++){pom+=Coms(i)->mean()*w(i);} return pom;};
	437	vec sample() {it_error ( "Not implemented" );return 0;}
	438	};
	439	*/
	440
	441	//! Uniform distributed density on a rectangular support
	442
	443	class euni: public epdf
	444	{
	445	protected:
	446	//! lower bound on support
	447	vec low;
	448	//! upper bound on support
	449	vec high;
	450	//! internal
	451	vec distance;
	452	//! normalizing coefficients
	453	double nk;
	454	//! cache of log( \c nk )
	455	double lnk;
	456	public:
	457	//! \name Constructors
	458	//!@{
	459	euni () : epdf () {}
	460	euni (const vec &low0, const vec &high0) {set_parameters (low0, high0);}
	461	void set_parameters (const vec &low0, const vec &high0) {
	462	distance = high0 - low0;
	463	it_assert_debug (min (distance) > 0.0, "bad support");
	464	low = low0;
	465	high = high0;
	466	nk = prod (1.0 / distance);
	467	lnk = log (nk);
	468	dim = low.length();
	469	}
	470	//!@}
	471
	472	double eval (const vec &val) const {return nk;}
	473	double evallog (const vec &val) const {
	474	if (any (val < low) && any (val > high)) {return inf;}
	475	else return lnk;
	476	}
	477	vec sample() const {
	478	vec smp (dim);
	479	#pragma omp critical
	480	UniRNG.sample_vector (dim , smp);
	481	return low + elem_mult (distance, smp);
	482	}
	483	//! set values of \c low and \c high
	484	vec mean() const {return (high -low) / 2.0;}
	485	vec variance() const {return (pow (high, 2) + pow (low, 2) + elem_mult (high, low)) / 3.0;}
	486	//! Load from structure with elements:
	487	//! \code
	488	//! { high = [...]; // vector of upper bounds
	489	//! low = [...]; // vector of lower bounds
	490	//! rv = {class="RV",...} // description of RV
	491	//! }
	492	//! \endcode
	493	//!@}
	494	void from_setting (const Setting &set) {
	495	epdf::from_setting (set); // reads rv and rvc
	496
	497	UI::get (high, set, "high", UI::compulsory);
	498	UI::get (low, set, "low", UI::compulsory);
	499	}
	500	};
	501
	502
	503	/*!
	504	\brief Normal distributed linear function with linear function of mean value;
	505
	506	Mean value \f$mu=A*rvc+mu_0\f$.
	507	*/
	508	template < class sq_T, template <typename> class TEpdf = enorm >
	509	class mlnorm : public mpdf_internal< TEpdf<sq_T> >
	510	{
	511	protected:
	512	//! Internal epdf that arise by conditioning on \c rvc
	513	mat A;
	514	vec mu_const;
	515	// vec& _mu; //cached epdf.mu; !!!!!! WHY NOT?
	516	public:
	517	//! \name Constructors
	518	//!@{
	519	mlnorm() : mpdf_internal< TEpdf<sq_T> >() {};
	520	mlnorm (const mat &A, const vec &mu0, const sq_T &R) : mpdf_internal< TEpdf<sq_T> >() {
	521	set_parameters (A, mu0, R);
	522	}
	523
	524	//! Set \c A and \c R
	525	void set_parameters (const mat &A0, const vec &mu0, const sq_T &R0) {
	526	it_assert_debug (A0.rows() == mu0.length(), "");
	527	it_assert_debug (A0.rows() == R0.rows(), "");
	528
	529	this->iepdf.set_parameters (zeros (A0.rows()), R0);
	530	A = A0;
	531	mu_const = mu0;
	532	this->dimc = A0.cols();
	533	}
	534	//!@}
	535	//! Set value of \c rvc . Result of this operation is stored in \c epdf use function \c _ep to access it.
	536	void condition (const vec &cond) {
	537	this->iepdf._mu() = A * cond + mu_const;
	538	//R is already assigned;
	539	}
	540
	541	//!access function
	542	vec& _mu_const() {return mu_const;}
	543	//!access function
	544	mat& _A() {return A;}
	545	//!access function
	546	mat _R() { return this->iepdf._R().to_mat(); }
	547
	548	template<typename sq_M>
	549	friend std::ostream &operator<< (std::ostream &os, mlnorm<sq_M, enorm> &ml);
	550
	551	void from_setting (const Setting &set) {
	552	mpdf::from_setting (set);
	553
	554	UI::get (A, set, "A", UI::compulsory);
	555	UI::get (mu_const, set, "const", UI::compulsory);
	556	mat R0;
	557	UI::get (R0, set, "R", UI::compulsory);
	558	set_parameters (A, mu_const, R0);
	559	};
	560	};
	561	UIREGISTER (mlnorm<ldmat>);
	562	UIREGISTER (mlnorm<fsqmat>);
	563	UIREGISTER (mlnorm<chmat>);
	564
	565	//! Mpdf with general function for mean value
	566	template<class sq_T>
	567	class mgnorm : public mpdf_internal< enorm< sq_T > >
	568	{
	569	protected:
	570	// vec μ WHY NOT?
	571	fnc* g;
	572	public:
	573	//!default constructor
	574	mgnorm() : mpdf_internal<enorm<sq_T> >() { }
	575	//!set mean function
	576	inline void set_parameters (fnc* g0, const sq_T &R0);
	577	inline void condition (const vec &cond);
	578
	579
	580	/*! UI for mgnorm
	581
	582	The mgnorm is constructed from a structure with fields:
	583	\code
	584	system = {
	585	type = "mgnorm";
	586	// function for mean value evolution
	587	g = {type="fnc"; ... }
	588
	589	// variance
	590	R = [1, 0,
	591	0, 1];
	592	// --OR --
	593	dR = [1, 1];
	594
	595	// == OPTIONAL ==
	596
	597	// description of y variables
	598	y = {type="rv"; names=["y", "u"];};
	599	// description of u variable
	600	u = {type="rv"; names=[];}
	601	};
	602	\endcode
	603
	604	Result if
	605	*/
	606
	607	void from_setting (const Setting &set) {
	608	fnc* g = UI::build<fnc> (set, "g", UI::compulsory);
	609
	610	mat R;
	611	vec dR;
	612	if (UI::get (dR, set, "dR"))
	613	R = diag (dR);
	614	else
	615	UI::get (R, set, "R", UI::compulsory);
	616
	617	set_parameters (g, R);
	618	}
	619	};
	620
	621	UIREGISTER (mgnorm<chmat>);
	622
	623
	624	/*! (Approximate) Student t density with linear function of mean value
	625
	626	The internal epdf of this class is of the type of a Gaussian (enorm).
	627	However, each conditioning is trying to assure the best possible approximation by taking into account the zeta function. See [] for reference.
	628
	629	Perhaps a moment-matching technique?
	630	*/
	631	class mlstudent : public mlnorm<ldmat, enorm>
	632	{
	633	protected:
	634	ldmat Lambda;
	635	ldmat &_R;
	636	ldmat Re;
	637	public:
	638	mlstudent () : mlnorm<ldmat, enorm> (),
	639	Lambda (), _R (iepdf._R()) {}
	640	void set_parameters (const mat &A0, const vec &mu0, const ldmat &R0, const ldmat& Lambda0) {
	641	it_assert_debug (A0.rows() == mu0.length(), "");
	642	it_assert_debug (R0.rows() == A0.rows(), "");
	643
	644	iepdf.set_parameters (mu0, Lambda); //
	645	A = A0;
	646	mu_const = mu0;
	647	Re = R0;
	648	Lambda = Lambda0;
	649	}
	650	void condition (const vec &cond) {
	651	iepdf._mu() = A * cond + mu_const;
	652	double zeta;
	653	//ugly hack!
	654	if ( (cond.length() + 1) == Lambda.rows()) {
	655	zeta = Lambda.invqform (concat (cond, vec_1 (1.0)));
	656	} else {
	657	zeta = Lambda.invqform (cond);
122	658	}
123		//!@}
124
125		//! \name Mathematical operations
126		//!@{
127
128		//! dupdate in exponential form (not really handy)
129		void dupdate ( mat &v,double nu=1.0 );
130
131		vec sample() const;
132
133		double evallog_nn ( const vec &val ) const;
134		double lognc () const;
135		vec mean() const {return mu;}
136		vec variance() const {return diag ( R.to_mat() );}
137		// mlnorm<sq_T>* condition ( const RV &rvn ) const ; <=========== fails to cmpile. Why?
138		mpdf* condition ( const RV &rvn ) const ;
139		enorm<sq_T>* marginal ( const RV &rv ) const;
140		// epdf* marginal ( const RV &rv ) const;
141		//!@}
142
143		//! \name Access to attributes
144		//!@{
145
146		vec& _mu() {return mu;}
147		void set_mu ( const vec mu0 ) { mu=mu0;}
148		sq_T& _R() {return R;}
149		const sq_T& _R() const {return R;}
150		//!@}
151
152		};
153		UIREGISTER(enorm<chmat>);
154		UIREGISTER(enorm<ldmat>);
155		UIREGISTER(enorm<fsqmat>);
156
157
158		/*!
159		* \brief Gauss-inverse-Wishart density stored in LD form
160
161		* For \f$p\f$-variate densities, given rv.count() should be \f$p\times\f$ V.rows().
162		*
163		*/
164		class egiw : public eEF
165		{
166		protected:
167		//! Extended information matrix of sufficient statistics
168		ldmat V;
169		//! Number of data records (degrees of freedom) of sufficient statistics
170		double nu;
171		//! Dimension of the output
172		int dimx;
173		//! Dimension of the regressor
174		int nPsi;
175		public:
176		//!\name Constructors
177		//!@{
178		egiw() :eEF() {};
179		egiw ( int dimx0, ldmat V0, double nu0=-1.0 ) :eEF() {set_parameters ( dimx0,V0, nu0 );};
180
181		void set_parameters ( int dimx0, ldmat V0, double nu0=-1.0 )
182		{
183		dimx=dimx0;
184		nPsi = V0.rows()-dimx;
185		dim = dimx* ( dimx+nPsi ); // size(R) + size(Theta)
186
187		V=V0;
188		if ( nu0<0 )
189		{
190		nu = 0.1 +nPsi +2*dimx +2; // +2 assures finite expected value of R
191		// terms before that are sufficient for finite normalization
192		}
193		else
194		{
195		nu=nu0;
	659	_R = Re;
	660	_R *= (1 + zeta);// / ( nu ); << nu is in Re!!!!!!
	661	};
	662
	663	};
	664	/*!
	665	\brief Gamma random walk
	666
	667	Mean value, \f$\mu\f$, of this density is given by \c rvc .
	668	Standard deviation of the random walk is proportional to one \f$k\f$-th the mean.
	669	This is achieved by setting \f$\alpha=k\f$ and \f$\beta=k/\mu\f$.
	670
	671	The standard deviation of the walk is then: \f$\mu/\sqrt(k)\f$.
	672	*/
	673	class mgamma : public mpdf_internal<egamma>
	674	{
	675	protected:
	676
	677	//! Constant \f$k\f$
	678	double k;
	679
	680	//! cache of iepdf.beta
	681	vec &_beta;
	682
	683	public:
	684	//! Constructor
	685	mgamma() : mpdf_internal<egamma>(), k (0),
	686	_beta (iepdf._beta()) {
	687	}
	688
	689	//! Set value of \c k
	690	void set_parameters (double k, const vec &beta0);
	691
	692	void condition (const vec &val) {_beta = k / val;};
	693	//! Load from structure with elements:
	694	//! \code
	695	//! { alpha = [...]; // vector of alpha
	696	//! k = 1.1; // multiplicative constant k
	697	//! rv = {class="RV",...} // description of RV
	698	//! rvc = {class="RV",...} // description of RV in condition
	699	//! }
	700	//! \endcode
	701	//!@}
	702	void from_setting (const Setting &set) {
	703	mpdf::from_setting (set); // reads rv and rvc
	704	vec betatmp; // ugly but necessary
	705	UI::get (betatmp, set, "beta", UI::compulsory);
	706	UI::get (k, set, "k", UI::compulsory);
	707	set_parameters (k, betatmp);
	708	}
	709	};
	710	UIREGISTER (mgamma);
	711
	712	/*!
	713	\brief Inverse-Gamma random walk
	714
	715	Mean value, \f$ \mu \f$, of this density is given by \c rvc .
	716	Standard deviation of the random walk is proportional to one \f$ k \f$-th the mean.
	717	This is achieved by setting \f$ \alpha=\mu/k^2+2 \f$ and \f$ \beta=\mu(\alpha-1)\f$.
	718
	719	The standard deviation of the walk is then: \f$ \mu/\sqrt(k)\f$.
	720	*/
	721	class migamma : public mpdf_internal<eigamma>
	722	{
	723	protected:
	724	//! Constant \f$k\f$
	725	double k;
	726
	727	//! cache of iepdf.alpha
	728	vec &_alpha;
	729
	730	//! cache of iepdf.beta
	731	vec &_beta;
	732
	733	public:
	734	//! \name Constructors
	735	//!@{
	736	migamma() : mpdf_internal<eigamma>(),
	737	k (0),
	738	_alpha (iepdf._alpha()),
	739	_beta (iepdf._beta()) {
	740	}
	741
	742	migamma (const migamma &m) : mpdf_internal<eigamma>(),
	743	k (0),
	744	_alpha (iepdf._alpha()),
	745	_beta (iepdf._beta()) {
	746	}
	747	//!@}
	748
	749	//! Set value of \c k
	750	void set_parameters (int len, double k0) {
	751	k = k0;
	752	iepdf.set_parameters ( (1.0 / (kk) + 2.0) ones (len) /alpha/, ones (len) /beta/);
	753	dimc = dimension();
	754	};
	755	void condition (const vec &val) {
	756	_beta = elem_mult (val, (_alpha - 1.0));
	757	};
	758	};
	759
	760
	761	/*!
	762	\brief Gamma random walk around a fixed point
	763
	764	Mean value, \f$\mu\f$, of this density is given by a geometric combination of \c rvc and given fixed point, \f$p\f$. \f$l\f$ is the coefficient of the geometric combimation
	765	\f[ \mu = \mu_{t-1} ^{l} p^{1-l}\f]
	766
	767	Standard deviation of the random walk is proportional to one \f$k\f$-th the mean.
	768	This is achieved by setting \f$\alpha=k\f$ and \f$\beta=k/\mu\f$.
	769
	770	The standard deviation of the walk is then: \f$\mu/\sqrt(k)\f$.
	771	*/
	772	class mgamma_fix : public mgamma
	773	{
	774	protected:
	775	//! parameter l
	776	double l;
	777	//! reference vector
	778	vec refl;
	779	public:
	780	//! Constructor
	781	mgamma_fix () : mgamma (), refl () {};
	782	//! Set value of \c k
	783	void set_parameters (double k0 , vec ref0, double l0) {
	784	mgamma::set_parameters (k0, ref0);
	785	refl = pow (ref0, 1.0 - l0);l = l0;
	786	dimc = dimension();
	787	};
	788
	789	void condition (const vec &val) {vec mean = elem_mult (refl, pow (val, l)); _beta = k / mean;};
	790	};
	791
	792
	793	/*!
	794	\brief Inverse-Gamma random walk around a fixed point
	795
	796	Mean value, \f$\mu\f$, of this density is given by a geometric combination of \c rvc and given fixed point, \f$p\f$. \f$l\f$ is the coefficient of the geometric combimation
	797	\f[ \mu = \mu_{t-1} ^{l} p^{1-l}\f]
	798
	799	==== Check == vv =
	800	Standard deviation of the random walk is proportional to one \f$k\f$-th the mean.
	801	This is achieved by setting \f$\alpha=k\f$ and \f$\beta=k/\mu\f$.
	802
	803	The standard deviation of the walk is then: \f$\mu/\sqrt(k)\f$.
	804	*/
	805	class migamma_ref : public migamma
	806	{
	807	protected:
	808	//! parameter l
	809	double l;
	810	//! reference vector
	811	vec refl;
	812	public:
	813	//! Constructor
	814	migamma_ref () : migamma (), refl () {};
	815	//! Set value of \c k
	816	void set_parameters (double k0 , vec ref0, double l0) {
	817	migamma::set_parameters (ref0.length(), k0);
	818	refl = pow (ref0, 1.0 - l0);
	819	l = l0;
	820	dimc = dimension();
	821	};
	822
	823	void condition (const vec &val) {
	824	vec mean = elem_mult (refl, pow (val, l));
	825	migamma::condition (mean);
	826	};
	827
	828	/*! UI for migamma_ref
	829
	830	The migamma_ref is constructed from a structure with fields:
	831	\code
	832	system = {
	833	type = "migamma_ref";
	834	ref = [1e-5; 1e-5; 1e-2 1e-3]; // reference vector
	835	l = 0.999; // constant l
	836	k = 0.1; // constant k
	837
	838	// == OPTIONAL ==
	839	// description of y variables
	840	y = {type="rv"; names=["y", "u"];};
	841	// description of u variable
	842	u = {type="rv"; names=[];}
	843	};
	844	\endcode
	845
	846	Result if
	847	*/
	848	void from_setting (const Setting &set);
	849
	850	// TODO dodelat void to_setting( Setting &set ) const;
	851	};
	852
	853
	854	UIREGISTER (migamma_ref);
	855
	856	/*! Log-Normal probability density
	857	only allow diagonal covariances!
	858
	859	Density of the form \f$ \log(x)\sim \mathcal{N}(\mu,\sigma^2) \f$ , i.e.
	860	\f[
	861	x \sim \frac{1}{x\sigma\sqrt{2\pi}}\exp{-\frac{1}{2\sigma^2}(\log(x)-\mu)}
	862	\f]
	863
	864	*/
	865	class elognorm: public enorm<ldmat>
	866	{
	867	public:
	868	vec sample() const {return exp (enorm<ldmat>::sample());};
	869	vec mean() const {vec var = enorm<ldmat>::variance();return exp (mu - 0.5*var);};
	870
	871	};
	872
	873	/*!
	874	\brief Log-Normal random walk
	875
	876	Mean value, \f$\mu\f$, is...
	877
	878	==== Check == vv =
	879	Standard deviation of the random walk is proportional to one \f$k\f$-th the mean.
	880	This is achieved by setting \f$\alpha=k\f$ and \f$\beta=k/\mu\f$.
	881
	882	The standard deviation of the walk is then: \f$\mu/\sqrt(k)\f$.
	883	*/
	884	class mlognorm : public mpdf_internal<elognorm>
	885	{
	886	protected:
	887	//! parameter 1/2*sigma^2
	888	double sig2;
	889
	890	//! access
	891	vec μ
	892	public:
	893	//! Constructor
	894	mlognorm() : mpdf_internal<elognorm>(),
	895	sig2 (0),
	896	mu (iepdf._mu()) {
	897	}
	898
	899	//! Set value of \c k
	900	void set_parameters (int size, double k) {
	901	sig2 = 0.5 * log (k * k + 1);
	902	iepdf.set_parameters (zeros (size), 2sig2eye (size));
	903
	904	dimc = size;
	905	};
	906
	907	void condition (const vec &val) {
	908	mu = log (val) - sig2;//elem_mult ( refl,pow ( val,l ) );
	909	};
	910
	911	/*! UI for mlognorm
	912
	913	The mlognorm is constructed from a structure with fields:
	914	\code
	915	system = {
	916	type = "mlognorm";
	917	k = 0.1; // constant k
	918	mu0 = [1., 1.];
	919
	920	// == OPTIONAL ==
	921	// description of y variables
	922	y = {type="rv"; names=["y", "u"];};
	923	// description of u variable
	924	u = {type="rv"; names=[];}
	925	};
	926	\endcode
	927
	928	*/
	929	void from_setting (const Setting &set);
	930
	931	// TODO dodelat void to_setting( Setting &set ) const;
	932
	933	};
	934
	935	UIREGISTER (mlognorm);
	936
	937	/*! inverse Wishart density defined on Choleski decomposition
	938
	939	*/
	940	class eWishartCh : public epdf
	941	{
	942	protected:
	943	//! Upper-Triagle of Choleski decomposition of \f$ \Psi \f$
	944	chmat Y;
	945	//! dimension of matrix \f$ \Psi \f$
	946	int p;
	947	//! degrees of freedom \f$ \nu \f$
	948	double delta;
	949	public:
	950	void set_parameters (const mat &Y0, const double delta0) {Y = chmat (Y0);delta = delta0; p = Y.rows(); dim = p * p; }
	951	mat sample_mat() const {
	952	mat X = zeros (p, p);
	953
	954	//sample diagonal
	955	for (int i = 0;i < p;i++) {
	956	GamRNG.setup (0.5* (delta - i) , 0.5); // no +1 !! index if from 0
	957	#pragma omp critical
	958	X (i, i) = sqrt (GamRNG());
	959	}
	960	//do the rest
	961	for (int i = 0;i < p;i++) {
	962	for (int j = i + 1;j < p;j++) {
	963	#pragma omp critical
	964	X (i, j) = NorRNG.sample();
196	965	}
197	966	}
198		//!@}
199
200		vec sample() const;
201		vec mean() const;
202		vec variance() const;
203
204		//! LS estimate of \f$\theta\f$
205		vec est_theta() const;
206
207		//! Covariance of the LS estimate
208		ldmat est_theta_cov() const;
209
210		void mean_mat ( mat &M, mat&R ) const;
211		//! In this instance, val= [theta, r]. For multivariate instances, it is stored columnwise val = [theta_1 theta_2 ... r_1 r_2 ]
212		double evallog_nn ( const vec &val ) const;
213		double lognc () const;
214		void pow ( double p ) {V=p;nu=p;};
215
216		//! \name Access attributes
217		//!@{
218
219		ldmat& _V() {return V;}
220		const ldmat& _V() const {return V;}
221		double& _nu() {return nu;}
222		const double& _nu() const {return nu;}
223		void from_setting(const Setting &set){
224		UI::get(nu, set, "nu", UI::compulsory);
225		UI::get(dimx, set, "dimx", UI::compulsory);
226		mat V;
227		UI::get(V,set,"V", UI::compulsory);
228		set_parameters(dimx, V, nu);
229		RV* rv=UI::build<RV>(set,"rv", UI::compulsory);
230		set_rv(*rv);
231		delete rv;
	967	return X*Y._Ch();// return upper triangular part of the decomposition
	968	}
	969	vec sample () const {
	970	return vec (sample_mat()._data(), p*p);
	971	}
	972	//! fast access function y0 will be copied into Y.Ch.
	973	void setY (const mat &Ch0) {copy_vector (dim, Ch0._data(), Y._Ch()._data());}
	974	//! fast access function y0 will be copied into Y.Ch.
	975	void _setY (const vec &ch0) {copy_vector (dim, ch0._data(), Y._Ch()._data()); }
	976	//! access function
	977	const chmat& getY() const {return Y;}
	978	};
	979
	980	class eiWishartCh: public epdf
	981	{
	982	protected:
	983	eWishartCh W;
	984	int p;
	985	double delta;
	986	public:
	987	void set_parameters (const mat &Y0, const double delta0) {
	988	delta = delta0;
	989	W.set_parameters (inv (Y0), delta0);
	990	dim = W.dimension(); p = Y0.rows();
	991	}
	992	vec sample() const {mat iCh; iCh = inv (W.sample_mat()); return vec (iCh._data(), dim);}
	993	void _setY (const vec &y0) {
	994	mat Ch (p, p);
	995	mat iCh (p, p);
	996	copy_vector (dim, y0._data(), Ch._data());
	997
	998	iCh = inv (Ch);
	999	W.setY (iCh);
	1000	}
	1001	virtual double evallog (const vec &val) const {
	1002	chmat X (p);
	1003	const chmat& Y = W.getY();
	1004
	1005	copy_vector (p*p, val._data(), X._Ch()._data());
	1006	chmat iX (p);X.inv (iX);
	1007	// compute
	1008	// \frac{ \|\Psi\|^{m/2}\|X\|^{-(m+p+1)/2}e^{-tr(\Psi X^{-1})/2} }{ 2^{mp/2}\Gamma_p(m/2)},
	1009	mat M = Y.to_mat() * iX.to_mat();
	1010
	1011	double log1 = 0.5 * p * (2 * Y.logdet()) - 0.5 * (delta + p + 1) * (2 * X.logdet()) - 0.5 * trace (M);
	1012	//Fixme! Multivariate gamma omitted!! it is ok for sampling, but not otherwise!!
	1013
	1014	/* if (0) {
	1015	mat XX=X.to_mat();
	1016	mat YY=Y.to_mat();
	1017
	1018	double log2 = 0.5plog(det(YY))-0.5(delta+p+1)log(det(XX))-0.5trace(YYinv(XX));
	1019	cout << log1 << "," << log2 << endl;
	1020	}*/
	1021	return log1;
	1022	};
	1023
	1024	};
	1025
	1026	class rwiWishartCh : public mpdf
	1027	{
	1028	protected:
	1029	eiWishartCh eiW;
	1030	//!square root of \f$ \nu-p-1 \f$ - needed for computation of \f$ \Psi \f$ from conditions
	1031	double sqd;
	1032	//reference point for diagonal
	1033	vec refl;
	1034	double l;
	1035	int p;
	1036
	1037	public:
	1038	rwiWishartCh() : eiW(),
	1039	sqd (0), l (0), p (0) {
	1040	set_ep (eiW);
	1041	}
	1042
	1043	void set_parameters (int p0, double k, vec ref0, double l0) {
	1044	p = p0;
	1045	double delta = 2 / (k * k) + p + 3;
	1046	sqd = sqrt (delta - p - 1);
	1047	l = l0;
	1048	refl = pow (ref0, 1 - l);
	1049
	1050	eiW.set_parameters (eye (p), delta);
	1051	dimc = eiW.dimension();
	1052	}
	1053	void condition (const vec &c) {
	1054	vec z = c;
	1055	int ri = 0;
	1056	for (int i = 0;i < p*p;i += (p + 1)) {//trace diagonal element
	1057	z (i) = pow (z (i), l) * refl (ri);
	1058	ri++;
232	1059	}
233		//!@}
234		};
235		UIREGISTER(egiw);
236
237		/*! \brief Dirichlet posterior density
238
239		Continuous Dirichlet density of \f$n\f$-dimensional variable \f$x\f$
240		\f[
241		f(x\|\beta) = \frac{\Gamma[\gamma]}{\prod_{i=1}^{n}\Gamma(\beta_i)} \prod_{i=1}^{n}x_i^{\beta_i-1}
242		\f]
243		where \f$\gamma=\sum_i \beta_i\f$.
244		*/
245		class eDirich: public eEF
246		{
247		protected:
248		//!sufficient statistics
249		vec beta;
250		public:
251		//!\name Constructors
252		//!@{
253
254		eDirich () : eEF ( ) {};
255		eDirich ( const eDirich &D0 ) : eEF () {set_parameters ( D0.beta );};
256		eDirich ( const vec &beta0 ) {set_parameters ( beta0 );};
257		void set_parameters ( const vec &beta0 )
258		{
259		beta= beta0;
260		dim = beta.length();
261		}
262		//!@}
263
264		vec sample() const {it_error ( "Not implemented" );return vec_1 ( 0.0 );};
265		vec mean() const {return beta/sum(beta);};
266		vec variance() const {double gamma =sum(beta); return elem_mult ( beta, ( beta+1 ) ) / ( gamma* ( gamma+1 ) );}
267		//! In this instance, val is ...
268		double evallog_nn ( const vec &val ) const
269		{
270		double tmp; tmp= ( beta-1 ) *log ( val );
271		// it_assert_debug ( std::isfinite ( tmp ),"Infinite value" );
272		return tmp;
273		};
274		double lognc () const
275		{
276		double tmp;
277		double gam=sum ( beta );
278		double lgb=0.0;
279		for ( int i=0;i<beta.length();i++ ) {lgb+=lgamma ( beta ( i ) );}
280		tmp= lgb-lgamma ( gam );
281		// it_assert_debug ( std::isfinite ( tmp ),"Infinite value" );
282		return tmp;
283		};
284		//!access function
285		vec& _beta() {return beta;}
286		//!Set internal parameters
287		};
288
289		//! \brief Estimator for Multinomial density
290		class multiBM : public BMEF
291		{
292		protected:
293		//! Conjugate prior and posterior
294		eDirich est;
295		//! Pointer inside est to sufficient statistics
296		vec β
297		public:
298		//!Default constructor
299		multiBM ( ) : BMEF ( ),est ( ),beta ( est._beta() )
300		{
301		if ( beta.length() >0 ) {last_lognc=est.lognc();}
302		else{last_lognc=0.0;}
303		}
304		//!Copy constructor
305		multiBM ( const multiBM &B ) : BMEF ( B ),est ( B.est ),beta ( est._beta() ) {}
306		//! Sets sufficient statistics to match that of givefrom mB0
307		void set_statistics ( const BM* mB0 ) {const multiBM* mB=dynamic_cast<const multiBM*> ( mB0 ); beta=mB->beta;}
308		void bayes ( const vec &dt )
309		{
310		if ( frg<1.0 ) {beta*=frg;last_lognc=est.lognc();}
311		beta+=dt;
312		if ( evalll ) {ll=est.lognc()-last_lognc;}
313		}
314		double logpred ( const vec &dt ) const
315		{
316		eDirich pred ( est );
317		vec &beta = pred._beta();
318
319		double lll;
320		if ( frg<1.0 )
321		{beta*=frg;lll=pred.lognc();}
322		else
323		if ( evalll ) {lll=last_lognc;}
324		else{lll=pred.lognc();}
325
326		beta+=dt;
327		return pred.lognc()-lll;
328		}
329		void flatten ( const BMEF* B )
330		{
331		const multiBM* E=dynamic_cast<const multiBM*> ( B );
332		// sum(beta) should be equal to sum(B.beta)
333		const vec &Eb=E->beta;//const_cast<multiBM*> ( E )->_beta();
334		beta*= ( sum ( Eb ) /sum ( beta ) );
335		if ( evalll ) {last_lognc=est.lognc();}
336		}
337		const epdf& posterior() const {return est;};
338		const eDirich* _e() const {return &est;};
339		void set_parameters ( const vec &beta0 )
340		{
341		est.set_parameters ( beta0 );
342		if ( evalll ) {last_lognc=est.lognc();}
343		}
344		};
345
346		/*!
347		\brief Gamma posterior density
348
349		Multivariate Gamma density as product of independent univariate densities.
350		\f[
351		f(x\|\alpha,\beta) = \prod f(x_i\|\alpha_i,\beta_i)
352		\f]
353		*/
354
355		class egamma : public eEF
356		{
357		protected:
358		//! Vector \f$\alpha\f$
359		vec alpha;
360		//! Vector \f$\beta\f$
361		vec beta;
362		public :
363		//! \name Constructors
364		//!@{
365		egamma ( ) :eEF ( ), alpha ( 0 ), beta ( 0 ) {};
366		egamma ( const vec &a, const vec &b ) {set_parameters ( a, b );};
367		void set_parameters ( const vec &a, const vec &b ) {alpha=a,beta=b;dim = alpha.length();};
368		//!@}
369
370		vec sample() const;
371		//! TODO: is it used anywhere?
372		// mat sample ( int N ) const;
373		double evallog ( const vec &val ) const;
374		double lognc () const;
375		//! Returns poiter to alpha and beta. Potentially dengerous: use with care!
376		vec& _alpha() {return alpha;}
377		vec& _beta() {return beta;}
378		vec mean() const {return elem_div ( alpha,beta );}
379		vec variance() const {return elem_div ( alpha,elem_mult ( beta,beta ) ); }
380
381		//! Load from structure with elements:
382		//! \code
383		//! { alpha = [...]; // vector of alpha
384		//! beta = [...]; // vector of beta
385		//! rv = {class="RV",...} // description
386		//! }
387		//! \endcode
388		//!@}
389		void from_setting(const Setting &set){
390		epdf::from_setting(set); // reads rv
391		UI::get(alpha,set,"alpha", UI::compulsory);
392		UI::get(beta,set,"beta", UI::compulsory);
393		validate();
394		}
395		void validate(){
396		it_assert(alpha.length() ==beta.length(), "parameters do not match");
397		dim =alpha.length();
398		}
399		};
400		UIREGISTER(egamma);
401		/*!
402		\brief Inverse-Gamma posterior density
403
404		Multivariate inverse-Gamma density as product of independent univariate densities.
405		\f[
406		f(x\|\alpha,\beta) = \prod f(x_i\|\alpha_i,\beta_i)
407		\f]
408
409		Vector \f$\beta\f$ has different meaning (in fact it is 1/beta as used in definition of iG)
410
411		Inverse Gamma can be converted to Gamma using \f[
412		x\sim iG(a,b) => 1/x\sim G(a,1/b)
413		\f]
414		This relation is used in sampling.
415		*/
416
417		class eigamma : public egamma
418		{
419		protected:
420		public :
421		//! \name Constructors
422		//! All constructors are inherited
423		//!@{
424		//!@}
425
426		vec sample() const {return 1.0/egamma::sample();};
427		//! Returns poiter to alpha and beta. Potentially dangerous: use with care!
428		vec mean() const {return elem_div ( beta,alpha-1 );}
429		vec variance() const {vec mea=mean(); return elem_div ( elem_mult ( mea,mea ),alpha-2 );}
430		};
431		/*
432		//! Weighted mixture of epdfs with external owned components.
433		class emix : public epdf {
434		protected:
	1060
	1061	eiW._setY (sqd*z);
	1062	}
	1063	};
	1064
	1065	//! Switch between various resampling methods.
	1066	enum RESAMPLING_METHOD { MULTINOMIAL = 0, STRATIFIED = 1, SYSTEMATIC = 3 };
	1067	/*!
	1068	\brief Weighted empirical density
	1069
	1070	Used e.g. in particle filters.
	1071	*/
	1072	class eEmp: public epdf
	1073	{
	1074	protected :
	1075	//! Number of particles
435	1076	int n;
436		vec &w;
437		Array<epdf*> Coms;
438		public:
439		//! Default constructor
440		emix ( const RV &rv, vec &w0): epdf(rv), n(w0.length()), w(w0), Coms(n) {};
441		void set_parameters( int &i, double wi, epdf* ep){w(i)=wi;Coms(i)=ep;}
442		vec mean(){vec pom; for(int i=0;i<n;i++){pom+=Coms(i)->mean()*w(i);} return pom;};
443		vec sample() {it_error ( "Not implemented" );return 0;}
444		};
445		*/
446
447		//! Uniform distributed density on a rectangular support
448
449		class euni: public epdf
450		{
451		protected:
452		//! lower bound on support
453		vec low;
454		//! upper bound on support
455		vec high;
456		//! internal
457		vec distance;
458		//! normalizing coefficients
459		double nk;
460		//! cache of log( \c nk )
461		double lnk;
462		public:
463		//! \name Constructors
464		//!@{
465		euni ( ) :epdf ( ) {}
466		euni ( const vec &low0, const vec &high0 ) {set_parameters ( low0,high0 );}
467		void set_parameters ( const vec &low0, const vec &high0 )
468		{
469		distance = high0-low0;
470		it_assert_debug ( min ( distance ) >0.0,"bad support" );
471		low = low0;
472		high = high0;
473		nk = prod ( 1.0/distance );
474		lnk = log ( nk );
475		dim = low.length();
476		}
477		//!@}
478
479		double eval ( const vec &val ) const {return nk;}
480		double evallog ( const vec &val ) const {
481		if (any(val<low) && any(val>high)) {return inf;}
482		else return lnk;
483		}
484		vec sample() const
485		{
486		vec smp ( dim );
487		#pragma omp critical
488		UniRNG.sample_vector ( dim ,smp );
489		return low+elem_mult ( distance,smp );
490		}
491		//! set values of \c low and \c high
492		vec mean() const {return ( high-low ) /2.0;}
493		vec variance() const {return ( pow ( high,2 ) +pow ( low,2 ) +elem_mult ( high,low ) ) /3.0;}
494		//! Load from structure with elements:
495		//! \code
496		//! { high = [...]; // vector of upper bounds
497		//! low = [...]; // vector of lower bounds
498		//! rv = {class="RV",...} // description of RV
499		//! }
500		//! \endcode
501		//!@}
502		void from_setting(const Setting &set){
503		epdf::from_setting(set); // reads rv and rvc
504
505		UI::get(high,set,"high", UI::compulsory);
506		UI::get(low,set,"low", UI::compulsory);
507		}
508		};
509
510
511		/*!
512		\brief Normal distributed linear function with linear function of mean value;
513
514		Mean value \f$mu=A*rvc+mu_0\f$.
515		*/
516		template<class sq_T, template <typename> class TEpdf =enorm>
517		class mlnorm : public mpdf_internal< TEpdf<sq_T> >
518		{
519		protected:
520		//! Internal epdf that arise by conditioning on \c rvc
521		mat A;
522		vec mu_const;
523		// vec& _mu; //cached epdf.mu; !!!!!! WHY NOT?
524		public:
525		//! \name Constructors
526		//!@{
527		mlnorm():mpdf_internal< TEpdf<sq_T> >() {};
528		mlnorm (const mat &A, const vec &mu0, const sq_T &R ) : mpdf_internal< TEpdf<sq_T> >()
529		{
530		set_parameters ( A,mu0,R );
531		}
532
533		//! Set \c A and \c R
534		void set_parameters ( const mat &A0, const vec &mu0, const sq_T &R0 ){
535		it_assert_debug ( A0.rows() ==mu0.length(),"" );
536		it_assert_debug ( A0.rows() ==R0.rows(),"" );
537
538		this->iepdf.set_parameters(zeros(A0.rows()), R0);
539		A = A0;
540		mu_const = mu0;
541		this->dimc = A0.cols();
542		}
543		//!@}
544		//! Set value of \c rvc . Result of this operation is stored in \c epdf use function \c _ep to access it.
545		void condition ( const vec &cond ){
546		this->iepdf._mu()= A*cond + mu_const;
547		//R is already assigned;
548		}
549
550		//!access function
551		vec& _mu_const() {return mu_const;}
552		//!access function
553		mat& _A() {return A;}
554		//!access function
555		mat _R(){ return this->iepdf._R().to_mat(); }
556
557		template<typename sq_M>
558		friend std::ostream &operator<< ( std::ostream &os, mlnorm<sq_M,enorm> &ml );
559
560		void from_setting(const Setting &set){
561		mpdf::from_setting(set);
562
563		UI::get(A,set,"A", UI::compulsory);
564		UI::get(mu_const,set,"const", UI::compulsory);
565		mat R0;
566		UI::get(R0,set,"R", UI::compulsory);
567		set_parameters(A,mu_const,R0);
568		};
569		};
570		UIREGISTER(mlnorm<ldmat>);
571		UIREGISTER(mlnorm<fsqmat>);
572		UIREGISTER(mlnorm<chmat>);
573
574		//! Mpdf with general function for mean value
575		template<class sq_T>
576		class mgnorm : public mpdf_internal< enorm< sq_T > >
577		{
578		protected:
579		// vec μ WHY NOT?
580		fnc* g;
581		public:
582		//!default constructor
583		mgnorm():mpdf_internal<enorm<sq_T> >() { }
584		//!set mean function
585		inline void set_parameters ( fnc* g0, const sq_T &R0 );
586		inline void condition ( const vec &cond );
587
588
589		/*! UI for mgnorm
590
591		The mgnorm is constructed from a structure with fields:
592		\code
593		system = {
594		type = "mgnorm";
595		// function for mean value evolution
596		g = {type="fnc"; ... }
597
598		// variance
599		R = [1, 0,
600		0, 1];
601		// --OR --
602		dR = [1, 1];
603
604		// == OPTIONAL ==
605
606		// description of y variables
607		y = {type="rv"; names=["y", "u"];};
608		// description of u variable
609		u = {type="rv"; names=[];}
610		};
611		\endcode
612
613		Result if
614		*/
615
616		void from_setting( const Setting &set )
617		{
618		fnc* g = UI::build<fnc>( set, "g", UI::compulsory );
619
620		mat R;
621		vec dR;
622		if ( UI::get( dR, set, "dR" ) )
623		R=diag(dR);
624		else
625		UI::get( R, set, "R", UI::compulsory);
626
627		set_parameters(g,R);
628		}
629		};
630
631		UIREGISTER(mgnorm<chmat>);
632
633
634		/*! (Approximate) Student t density with linear function of mean value
635
636		The internal epdf of this class is of the type of a Gaussian (enorm).
637		However, each conditioning is trying to assure the best possible approximation by taking into account the zeta function. See [] for reference.
638
639		Perhaps a moment-matching technique?
640		*/
641		class mlstudent : public mlnorm<ldmat,enorm>
642		{
643		protected:
644		ldmat Lambda;
645		ldmat &_R;
646		ldmat Re;
647		public:
648		mlstudent ( ) :mlnorm<ldmat,enorm> (),
649		Lambda (), _R ( iepdf._R() ) {}
650		void set_parameters ( const mat &A0, const vec &mu0, const ldmat &R0, const ldmat& Lambda0 )
651		{
652		it_assert_debug ( A0.rows() ==mu0.length(),"" );
653		it_assert_debug ( R0.rows() ==A0.rows(),"" );
654
655		iepdf.set_parameters ( mu0,Lambda ); //
656		A = A0;
657		mu_const = mu0;
658		Re=R0;
659		Lambda = Lambda0;
660		}
661		void condition ( const vec &cond )
662		{
663		iepdf._mu() = A*cond + mu_const;
664		double zeta;
665		//ugly hack!
666		if ( ( cond.length() +1 ) ==Lambda.rows() )
667		{
668		zeta = Lambda.invqform ( concat ( cond, vec_1 ( 1.0 ) ) );
669		}
670		else
671		{
672		zeta = Lambda.invqform ( cond );
673		}
674		_R = Re;
675		_R*= ( 1+zeta );// / ( nu ); << nu is in Re!!!!!!
676		};
677
678		};
679		/*!
680		\brief Gamma random walk
681
682		Mean value, \f$\mu\f$, of this density is given by \c rvc .
683		Standard deviation of the random walk is proportional to one \f$k\f$-th the mean.
684		This is achieved by setting \f$\alpha=k\f$ and \f$\beta=k/\mu\f$.
685
686		The standard deviation of the walk is then: \f$\mu/\sqrt(k)\f$.
687		*/
688		class mgamma : public mpdf_internal<egamma>
689		{
690		protected:
691
692		//! Constant \f$k\f$
693		double k;
694
695		//! cache of iepdf.beta
696		vec &_beta;
697
698		public:
699		//! Constructor
700		mgamma():mpdf_internal<egamma>(), k(0),
701		_beta(iepdf._beta()) {
702		}
703
704		//! Set value of \c k
705		void set_parameters(double k, const vec &beta0);
706
707		void condition ( const vec &val ) {_beta=k/val;};
708		//! Load from structure with elements:
709		//! \code
710		//! { alpha = [...]; // vector of alpha
711		//! k = 1.1; // multiplicative constant k
712		//! rv = {class="RV",...} // description of RV
713		//! rvc = {class="RV",...} // description of RV in condition
714		//! }
715		//! \endcode
716		//!@}
717		void from_setting(const Setting &set){
718		mpdf::from_setting(set); // reads rv and rvc
719		vec betatmp; // ugly but necessary
720		UI::get(betatmp,set,"beta", UI::compulsory);
721		UI::get(k,set,"k", UI::compulsory);
722		set_parameters(k,betatmp);
723		}
724		};
725		UIREGISTER(mgamma);
726
727		/*!
728		\brief Inverse-Gamma random walk
729
730		Mean value, \f$ \mu \f$, of this density is given by \c rvc .
731		Standard deviation of the random walk is proportional to one \f$ k \f$-th the mean.
732		This is achieved by setting \f$ \alpha=\mu/k^2+2 \f$ and \f$ \beta=\mu(\alpha-1)\f$.
733
734		The standard deviation of the walk is then: \f$ \mu/\sqrt(k)\f$.
735		*/
736		class migamma : public mpdf_internal<eigamma>
737		{
738		protected:
739		//! Constant \f$k\f$
740		double k;
741
742		//! cache of iepdf.alpha
743		vec &_alpha;
744
745		//! cache of iepdf.beta
746		vec &_beta;
747
748		public:
749		//! \name Constructors
750		//!@{
751		migamma():mpdf_internal<eigamma>(),
752		k(0),
753		_alpha(iepdf._alpha()),
754		_beta(iepdf._beta()) {
755		}
756
757		migamma(const migamma &m):mpdf_internal<eigamma>(),
758		k(0),
759		_alpha(iepdf._alpha()),
760		_beta(iepdf._beta()) {
761		}
762		//!@}
763
764		//! Set value of \c k
765		void set_parameters ( int len, double k0 )
766		{
767		k=k0;
768		iepdf.set_parameters ( ( 1.0/ ( kk ) +2.0 ) ones ( len ) /alpha/, ones ( len ) /beta/ );
769		dimc = dimension();
770		};
771		void condition ( const vec &val )
772		{
773		_beta=elem_mult ( val, ( _alpha-1.0 ) );
774		};
775		};
776
777
778		/*!
779		\brief Gamma random walk around a fixed point
780
781		Mean value, \f$\mu\f$, of this density is given by a geometric combination of \c rvc and given fixed point, \f$p\f$. \f$l\f$ is the coefficient of the geometric combimation
782		\f[ \mu = \mu_{t-1} ^{l} p^{1-l}\f]
783
784		Standard deviation of the random walk is proportional to one \f$k\f$-th the mean.
785		This is achieved by setting \f$\alpha=k\f$ and \f$\beta=k/\mu\f$.
786
787		The standard deviation of the walk is then: \f$\mu/\sqrt(k)\f$.
788		*/
789		class mgamma_fix : public mgamma
790		{
791		protected:
792		//! parameter l
793		double l;
794		//! reference vector
795		vec refl;
796		public:
797		//! Constructor
798		mgamma_fix ( ) : mgamma ( ),refl () {};
799		//! Set value of \c k
800		void set_parameters ( double k0 , vec ref0, double l0 )
801		{
802		mgamma::set_parameters ( k0, ref0 );
803		refl=pow ( ref0,1.0-l0 );l=l0;
804		dimc=dimension();
805		};
806
807		void condition ( const vec &val ) {vec mean=elem_mult ( refl,pow ( val,l ) ); _beta=k/mean;};
808		};
809
810
811		/*!
812		\brief Inverse-Gamma random walk around a fixed point
813
814		Mean value, \f$\mu\f$, of this density is given by a geometric combination of \c rvc and given fixed point, \f$p\f$. \f$l\f$ is the coefficient of the geometric combimation
815		\f[ \mu = \mu_{t-1} ^{l} p^{1-l}\f]
816
817		==== Check == vv =
818		Standard deviation of the random walk is proportional to one \f$k\f$-th the mean.
819		This is achieved by setting \f$\alpha=k\f$ and \f$\beta=k/\mu\f$.
820
821		The standard deviation of the walk is then: \f$\mu/\sqrt(k)\f$.
822		*/
823		class migamma_ref : public migamma
824		{
825		protected:
826		//! parameter l
827		double l;
828		//! reference vector
829		vec refl;
830		public:
831		//! Constructor
832		migamma_ref ( ) : migamma (),refl ( ) {};
833		//! Set value of \c k
834		void set_parameters ( double k0 , vec ref0, double l0 )
835		{
836		migamma::set_parameters ( ref0.length(), k0 );
837		refl=pow ( ref0,1.0-l0 );
838		l=l0;
839		dimc = dimension();
840		};
841
842		void condition ( const vec &val )
843		{
844		vec mean=elem_mult ( refl,pow ( val,l ) );
845		migamma::condition ( mean );
846		};
847
848		/*! UI for migamma_ref
849
850		The migamma_ref is constructed from a structure with fields:
851		\code
852		system = {
853		type = "migamma_ref";
854		ref = [1e-5; 1e-5; 1e-2 1e-3]; // reference vector
855		l = 0.999; // constant l
856		k = 0.1; // constant k
857
858		// == OPTIONAL ==
859		// description of y variables
860		y = {type="rv"; names=["y", "u"];};
861		// description of u variable
862		u = {type="rv"; names=[];}
863		};
864		\endcode
865
866		Result if
867		*/
868		void from_setting( const Setting &set );
869
870		// TODO dodelat void to_setting( Setting &set ) const;
871		};
872
873
874		UIREGISTER(migamma_ref);
875
876		/*! Log-Normal probability density
877		only allow diagonal covariances!
878
879		Density of the form \f$ \log(x)\sim \mathcal{N}(\mu,\sigma^2) \f$ , i.e.
880		\f[
881		x \sim \frac{1}{x\sigma\sqrt{2\pi}}\exp{-\frac{1}{2\sigma^2}(\log(x)-\mu)}
882		\f]
883
884		*/
885		class elognorm: public enorm<ldmat>
886		{
887		public:
888		vec sample() const {return exp ( enorm<ldmat>::sample() );};
889		vec mean() const {vec var=enorm<ldmat>::variance();return exp ( mu - 0.5*var );};
890
891		};
892
893		/*!
894		\brief Log-Normal random walk
895
896		Mean value, \f$\mu\f$, is...
897
898		==== Check == vv =
899		Standard deviation of the random walk is proportional to one \f$k\f$-th the mean.
900		This is achieved by setting \f$\alpha=k\f$ and \f$\beta=k/\mu\f$.
901
902		The standard deviation of the walk is then: \f$\mu/\sqrt(k)\f$.
903		*/
904		class mlognorm : public mpdf_internal<elognorm>
905		{
906		protected:
907		//! parameter 1/2*sigma^2
908		double sig2;
909
910		//! access
911		vec μ
912		public:
913		//! Constructor
914		mlognorm():mpdf_internal<elognorm>(),
915		sig2(0),
916		mu(iepdf._mu()) {
917		}
918
919		//! Set value of \c k
920		void set_parameters ( int size, double k )
921		{
922		sig2 = 0.5log ( kk+1 );
923		iepdf.set_parameters ( zeros ( size ),2sig2eye ( size ) );
924
925		dimc = size;
926		};
927
928		void condition ( const vec &val )
929		{
930		mu=log ( val )-sig2;//elem_mult ( refl,pow ( val,l ) );
931		};
932
933		/*! UI for mlognorm
934
935		The mlognorm is constructed from a structure with fields:
936		\code
937		system = {
938		type = "mlognorm";
939		k = 0.1; // constant k
940		mu0 = [1., 1.];
941
942		// == OPTIONAL ==
943		// description of y variables
944		y = {type="rv"; names=["y", "u"];};
945		// description of u variable
946		u = {type="rv"; names=[];}
947		};
948		\endcode
949
950		*/
951		void from_setting( const Setting &set );
952
953		// TODO dodelat void to_setting( Setting &set ) const;
954
955		};
956
957		UIREGISTER(mlognorm);
958
959		/*! inverse Wishart density defined on Choleski decomposition
960
961		*/
962		class eWishartCh : public epdf
963		{
964		protected:
965		//! Upper-Triagle of Choleski decomposition of \f$ \Psi \f$
966		chmat Y;
967		//! dimension of matrix \f$ \Psi \f$
968		int p;
969		//! degrees of freedom \f$ \nu \f$
970		double delta;
971		public:
972		void set_parameters ( const mat &Y0, const double delta0 ) {Y=chmat ( Y0 );delta=delta0; p=Y.rows(); dim = p*p; }
973		mat sample_mat() const
974		{
975		mat X=zeros ( p,p );
976
977		//sample diagonal
978		for ( int i=0;i<p;i++ )
979		{
980		GamRNG.setup ( 0.5* ( delta-i ) , 0.5 ); // no +1 !! index if from 0
981		#pragma omp critical
982		X ( i,i ) =sqrt ( GamRNG() );
983		}
984		//do the rest
985		for ( int i=0;i<p;i++ )
986		{
987		for ( int j=i+1;j<p;j++ )
988		{
989		#pragma omp critical
990		X ( i,j ) =NorRNG.sample();
991		}
992		}
993		return X*Y._Ch();// return upper triangular part of the decomposition
994		}
995		vec sample () const
996		{
997		return vec ( sample_mat()._data(),p*p );
998		}
999		//! fast access function y0 will be copied into Y.Ch.
1000		void setY ( const mat &Ch0 ) {copy_vector ( dim,Ch0._data(), Y._Ch()._data() );}
1001		//! fast access function y0 will be copied into Y.Ch.
1002		void _setY ( const vec &ch0 ) {copy_vector ( dim, ch0._data(), Y._Ch()._data() ); }
1003		//! access function
1004		const chmat& getY()const {return Y;}
1005		};
1006
1007		class eiWishartCh: public epdf
1008		{
1009		protected:
1010		eWishartCh W;
1011		int p;
1012		double delta;
1013		public:
1014		void set_parameters ( const mat &Y0, const double delta0) {
1015		delta = delta0;
1016		W.set_parameters ( inv ( Y0 ),delta0 );
1017		dim = W.dimension(); p=Y0.rows();
1018		}
1019		vec sample() const {mat iCh; iCh=inv ( W.sample_mat() ); return vec ( iCh._data(),dim );}
1020		void _setY ( const vec &y0 )
1021		{
1022		mat Ch ( p,p );
1023		mat iCh ( p,p );
1024		copy_vector ( dim, y0._data(), Ch._data() );
1025
1026		iCh=inv ( Ch );
1027		W.setY ( iCh );
1028		}
1029		virtual double evallog ( const vec &val ) const {
1030		chmat X(p);
1031		const chmat& Y=W.getY();
1032
1033		copy_vector(p*p,val._data(),X._Ch()._data());
1034		chmat iX(p);X.inv(iX);
1035		// compute
1036		// \frac{ \|\Psi\|^{m/2}\|X\|^{-(m+p+1)/2}e^{-tr(\Psi X^{-1})/2} }{ 2^{mp/2}\Gamma_p(m/2)},
1037		mat M=Y.to_mat()*iX.to_mat();
1038
1039		double log1 = 0.5p(2Y.logdet())-0.5(delta+p+1)(2X.logdet())-0.5*trace(M);
1040		//Fixme! Multivariate gamma omitted!! it is ok for sampling, but not otherwise!!
1041
1042		/* if (0) {
1043		mat XX=X.to_mat();
1044		mat YY=Y.to_mat();
1045
1046		double log2 = 0.5plog(det(YY))-0.5(delta+p+1)log(det(XX))-0.5trace(YYinv(XX));
1047		cout << log1 << "," << log2 << endl;
1048		}*/
1049		return log1;
1050		};
1051
1052		};
1053
1054		class rwiWishartCh : public mpdf
1055		{
1056		protected:
1057		eiWishartCh eiW;
1058		//!square root of \f$ \nu-p-1 \f$ - needed for computation of \f$ \Psi \f$ from conditions
1059		double sqd;
1060		//reference point for diagonal
1061		vec refl;
1062		double l;
1063		int p;
1064
1065		public:
1066		rwiWishartCh():eiW(),
1067		sqd(0), l(0), p(0) {
1068		set_ep(eiW);
1069		}
1070
1071		void set_parameters ( int p0, double k, vec ref0, double l0 )
1072		{
1073		p=p0;
1074		double delta = 2/(k*k)+p+3;
1075		sqd=sqrt ( delta-p-1 );
1076		l=l0;
1077		refl=pow(ref0,1-l);
1078
1079		eiW.set_parameters ( eye ( p ),delta );
1080		dimc=eiW.dimension();
1081		}
1082		void condition ( const vec &c ) {
1083		vec z=c;
1084		int ri=0;
1085		for(int i=0;i<p*p;i+=(p+1)){//trace diagonal element
1086		z(i) = pow(z(i),l)*refl(ri);
1087		ri++;
1088		}
1089
1090		eiW._setY ( sqd*z );
1091		}
1092		};
1093
1094		//! Switch between various resampling methods.
1095		enum RESAMPLING_METHOD { MULTINOMIAL = 0, STRATIFIED = 1, SYSTEMATIC = 3 };
1096		/*!
1097		\brief Weighted empirical density
1098
1099		Used e.g. in particle filters.
1100		*/
1101		class eEmp: public epdf
1102		{
1103		protected :
1104		//! Number of particles
1105		int n;
1106		//! Sample weights \f$w\f$
1107		vec w;
1108		//! Samples \f$x^{(i)}, i=1..n\f$
1109		Array<vec> samples;
1110		public:
1111		//! \name Constructors
1112		//!@{
1113		eEmp ( ) :epdf ( ),w ( ),samples ( ) {};
1114		//! copy constructor
1115		eEmp ( const eEmp &e ) : epdf ( e ), w ( e.w ), samples ( e.samples ) {};
1116		//!@}
1117
1118		//! Set samples and weights
1119		void set_statistics ( const vec &w0, const epdf* pdf0 );
1120		//! Set samples and weights
1121		void set_statistics ( const epdf* pdf0 , int n ) {set_statistics ( ones ( n ) /n,pdf0 );};
1122		//! Set sample
1123		void set_samples ( const epdf* pdf0 );
1124		//! Set sample
1125		void set_parameters ( int n0, bool copy=true ) {n=n0; w.set_size ( n0,copy );samples.set_size ( n0,copy );};
1126		//! Potentially dangerous, use with care.
1127		vec& _w() {return w;};
1128		//! Potentially dangerous, use with care.
1129		const vec& _w() const {return w;};
1130		//! access function
1131		Array<vec>& _samples() {return samples;};
1132		//! access function
1133		const Array<vec>& _samples() const {return samples;};
1134		//! 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.
1135		ivec resample ( RESAMPLING_METHOD method=SYSTEMATIC );
1136		//! inherited operation : NOT implemneted
1137		vec sample() const {it_error ( "Not implemented" );return 0;}
1138		//! inherited operation : NOT implemneted
1139		double evallog ( const vec &val ) const {it_error ( "Not implemented" );return 0.0;}
1140		vec mean() const
1141		{
1142		vec pom=zeros ( dim );
1143		for ( int i=0;i<n;i++ ) {pom+=samples ( i ) *w ( i );}
1144		return pom;
1145		}
1146		vec variance() const
1147		{
1148		vec pom=zeros ( dim );
1149		for ( int i=0;i<n;i++ ) {pom+=pow ( samples ( i ),2 ) *w ( i );}
1150		return pom-pow ( mean(),2 );
1151		}
1152		//! For this class, qbounds are minimum and maximum value of the population!
1153		void qbounds ( vec &lb, vec &ub, double perc=0.95 ) const
1154		{
1155		// lb in inf so than it will be pushed below;
1156		lb.set_size ( dim );
1157		ub.set_size ( dim );
1158		lb = std::numeric_limits<double>::infinity();
1159		ub = -std::numeric_limits<double>::infinity();
1160		int j;
1161		for ( int i=0;i<n;i++ )
1162		{
1163		for ( j=0;j<dim; j++ )
1164		{
1165		if ( samples ( i ) ( j ) <lb ( j ) ) {lb ( j ) =samples ( i ) ( j );}
1166		if ( samples ( i ) ( j ) >ub ( j ) ) {ub ( j ) =samples ( i ) ( j );}
1167		}
	1077	//! Sample weights \f$w\f$
	1078	vec w;
	1079	//! Samples \f$x^{(i)}, i=1..n\f$
	1080	Array<vec> samples;
	1081	public:
	1082	//! \name Constructors
	1083	//!@{
	1084	eEmp () : epdf (), w (), samples () {};
	1085	//! copy constructor
	1086	eEmp (const eEmp &e) : epdf (e), w (e.w), samples (e.samples) {};
	1087	//!@}
	1088
	1089	//! Set samples and weights
	1090	void set_statistics (const vec &w0, const epdf &pdf0);
	1091	//! Set samples and weights
	1092	void set_statistics (const epdf &pdf0 , int n) {set_statistics (ones (n) / n, pdf0);};
	1093	//! Set sample
	1094	void set_samples (const epdf* pdf0);
	1095	//! Set sample
	1096	void set_parameters (int n0, bool copy = true) {n = n0; w.set_size (n0, copy);samples.set_size (n0, copy);};
	1097	//! Potentially dangerous, use with care.
	1098	vec& _w() {return w;};
	1099	//! Potentially dangerous, use with care.
	1100	const vec& _w() const {return w;};
	1101	//! access function
	1102	Array<vec>& _samples() {return samples;};
	1103	//! access function
	1104	const Array<vec>& _samples() const {return samples;};
	1105	//! 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.
	1106	ivec resample (RESAMPLING_METHOD method = SYSTEMATIC);
	1107	//! inherited operation : NOT implemneted
	1108	vec sample() const {it_error ("Not implemented");return 0;}
	1109	//! inherited operation : NOT implemneted
	1110	double evallog (const vec &val) const {it_error ("Not implemented");return 0.0;}
	1111	vec mean() const {
	1112	vec pom = zeros (dim);
	1113	for (int i = 0;i < n;i++) {pom += samples (i) * w (i);}
	1114	return pom;
	1115	}
	1116	vec variance() const {
	1117	vec pom = zeros (dim);
	1118	for (int i = 0;i < n;i++) {pom += pow (samples (i), 2) * w (i);}
	1119	return pom -pow (mean(), 2);
	1120	}
	1121	//! For this class, qbounds are minimum and maximum value of the population!
	1122	void qbounds (vec &lb, vec &ub, double perc = 0.95) const {
	1123	// lb in inf so than it will be pushed below;
	1124	lb.set_size (dim);
	1125	ub.set_size (dim);
	1126	lb = std::numeric_limits<double>::infinity();
	1127	ub = -std::numeric_limits<double>::infinity();
	1128	int j;
	1129	for (int i = 0;i < n;i++) {
	1130	for (j = 0;j < dim; j++) {
	1131	if (samples (i) (j) < lb (j)) {lb (j) = samples (i) (j);}
	1132	if (samples (i) (j) > ub (j)) {ub (j) = samples (i) (j);}
1168	1133	}
1169	1134	}
1170		};
	1135	}
	1136	};
1171	1137
1172	1138
1173	1139	////////////////////////
1174	1140
1175		template<class sq_T>
1176		~~void enorm<sq_T>::set_parameters ( const vec &mu0, const sq_T &R0~~ )
1177		{
	1141	template<class sq_T>
	1142	void enorm<sq_T>::set_parameters (const vec &mu0, const sq_T &R0)
	1143	{
1178	1144	//Fixme test dimensions of mu0 and R0;
1179		mu = mu0;
1180		R = R0;
1181		validate();
1182		};
1183
1184		template<class sq_T>
1185		void enorm<sq_T>::from_setting(const Setting &set){
1186		epdf::from_setting(set); //reads rv
1187
1188		UI::get(mu,set,"mu", UI::compulsory);
1189		mat Rtmp;// necessary for conversion
1190		UI::get(Rtmp,set,"R", UI::compulsory);
1191		R=Rtmp; // conversion
1192		validate();
1193		}
1194
1195		template<class sq_T>
1196		void enorm<sq_T>::dupdate ( mat &v, double nu )
1197		{
1198		//
1199		};
	1145	mu = mu0;
	1146	R = R0;
	1147	validate();
	1148	};
	1149
	1150	template<class sq_T>
	1151	void enorm<sq_T>::from_setting (const Setting &set)
	1152	{
	1153	epdf::from_setting (set); //reads rv
	1154
	1155	UI::get (mu, set, "mu", UI::compulsory);
	1156	mat Rtmp;// necessary for conversion
	1157	UI::get (Rtmp, set, "R", UI::compulsory);
	1158	R = Rtmp; // conversion
	1159	validate();
	1160	}
	1161
	1162	template<class sq_T>
	1163	void enorm<sq_T>::dupdate (mat &v, double nu)
	1164	{
	1165	//
	1166	};
1200	1167
1201	1168	// template<class sq_T>
…	…
1204	1171	// };
1205	1172
1206		template<class sq_T>
1207		vec enorm<sq_T>::sample() const
1208		{
1209		~~vec x ( dim~~ );
	1173	template<class sq_T>
	1174	vec enorm<sq_T>::sample() const
	1175	{
	1176	vec x (dim);
1210	1177	#pragma omp critical
1211		~~NorRNG.sample_vector ( dim,x~~ );
1212		~~vec smp = R.sqrt_mult ( x~~ );
1213
1214		smp += mu;
1215		return smp;
1216		};
	1178	NorRNG.sample_vector (dim, x);
	1179	vec smp = R.sqrt_mult (x);
	1180
	1181	smp += mu;
	1182	return smp;
	1183	};
1217	1184
1218	1185	// template<class sq_T>
…	…
1224	1191	// };
1225	1192
1226		template<class sq_T>
1227		~~double enorm<sq_T>::evallog_nn ( const vec &val~~ ) const
1228		{
1229		// 1.83787706640935 = log(2pi)
1230		~~double tmp=-0.5* ( R.invqform ( mu-val )~~ );// - lognc();
1231		return tmp;
1232		};
1233
1234		template<class sq_T>
1235		inline double enorm<sq_T>::lognc () const
1236		{
1237		// 1.83787706640935 = log(2pi)
1238		~~double tmp=0.5* ( R.cols() * 1.83787706640935 +R.logdet()~~ );
1239		return tmp;
1240		};
	1193	template<class sq_T>
	1194	double enorm<sq_T>::evallog_nn (const vec &val) const
	1195	{
	1196	// 1.83787706640935 = log(2pi)
	1197	double tmp = -0.5 * (R.invqform (mu - val));// - lognc();
	1198	return tmp;
	1199	};
	1200
	1201	template<class sq_T>
	1202	inline double enorm<sq_T>::lognc () const
	1203	{
	1204	// 1.83787706640935 = log(2pi)
	1205	double tmp = 0.5 * (R.cols() * 1.83787706640935 + R.logdet());
	1206	return tmp;
	1207	};
1241	1208
1242	1209
…	…
1267	1234
1268	1235
1269		template<class sq_T>
1270		~~enorm<sq_T>* enorm<sq_T>::marginal ( const RV &rvn~~ ) const
1271		{
1272		~~it_assert_debug ( isnamed(), "rv description is not assigned"~~ );
1273		~~ivec irvn = rvn.dataind ( rv~~ );
1274
1275		~~sq_T Rn ( R,irvn );~~ //select rows and columns of R
1276
1277		enorm<sq_T>* tmp = new enorm<sq_T>;
1278		~~tmp->set_rv ( rvn~~ );
1279		~~tmp->set_parameters ( mu ( irvn ), Rn~~ );
1280		return tmp;
1281		}
1282
1283		template<class sq_T>
1284		~~mpdf* enorm<sq_T>::condition ( const RV &rvn~~ ) const
1285		{
1286
1287		~~it_assert_debug ( isnamed(),"rvs are not assigned"~~ );
1288
1289		~~RV rvc = rv.subt ( rvn~~ );
1290		~~it_assert_debug ( ( rvc._dsize() +rvn._dsize() ==rv._dsize() ),"wrong rvn"~~ );
1291		//Permutation vector of the new R
1292		~~ivec irvn = rvn.dataind ( rv~~ );
1293		~~ivec irvc = rvc.dataind ( rv~~ );
1294		~~ivec perm=concat ( irvn , irvc~~ );
1295		~~sq_T Rn ( R,perm~~ );
1296
1297		//fixme - could this be done in general for all sq_T?
1298		~~mat S=~~Rn.to_mat();
1299		//fixme
1300		~~int n=rvn._dsize()-~~1;
1301		~~int end=R.rows()-~~1;
1302		~~mat S11 = S.get ( 0,n, 0, n~~ );
1303		~~mat S12 = S.get ( 0, n , rvn._dsize(), end~~ );
1304		~~mat S22 = S.get ( rvn._dsize(), end, rvn._dsize(), end~~ );
1305
1306		~~vec mu1 = mu ( irvn~~ );
1307		~~vec mu2 = mu ( irvc~~ );
1308		~~mat A=S12*inv ( S22~~ );
1309		~~sq_T R_n ( S11 - A *S12.T()~~ );
1310
1311		~~mlnorm<sq_T>* tmp=new mlnorm<sq_T> (~~ );
1312		~~tmp->set_rv ( rvn ); tmp->set_rvc ( rvc~~ );
1313		~~tmp->set_parameters ( A,mu1-A*mu2,R_n~~ );
1314		return tmp;
1315		}
1316
1317		////
1318		///////
1319		template<class sq_T>
1320		~~void mgnorm<sq_T >::set_parameters ( fnc* g0, const sq_T &R0 ) {g=g0; this->iepdf.set_parameters ( zeros ( g->dimension() ), R0~~ );}
1321		template<class sq_T>
1322		~~void mgnorm<sq_T >::condition ( const vec &cond ) {this->iepdf._mu()=g->eval ( cond~~ );};
1323
1324		template<class sq_T>
1325		~~std::ostream &operator<< ( std::ostream &os, mlnorm<sq_T> &ml~~ )
1326		{
1327		~~os << "A:"<< ml.A<<~~endl;
1328		~~os << "mu:"<< ml.mu_const<<~~endl;
1329		~~os << "R:" << ml.iepdf->_R().to_mat() <<~~endl;
1330		return os;
1331		};
	1236	template<class sq_T>
	1237	enorm<sq_T>* enorm<sq_T>::marginal (const RV &rvn) const
	1238	{
	1239	it_assert_debug (isnamed(), "rv description is not assigned");
	1240	ivec irvn = rvn.dataind (rv);
	1241
	1242	sq_T Rn (R, irvn); //select rows and columns of R
	1243
	1244	enorm<sq_T>* tmp = new enorm<sq_T>;
	1245	tmp->set_rv (rvn);
	1246	tmp->set_parameters (mu (irvn), Rn);
	1247	return tmp;
	1248	}
	1249
	1250	template<class sq_T>
	1251	mpdf* enorm<sq_T>::condition (const RV &rvn) const
	1252	{
	1253
	1254	it_assert_debug (isnamed(), "rvs are not assigned");
	1255
	1256	RV rvc = rv.subt (rvn);
	1257	it_assert_debug ( (rvc._dsize() + rvn._dsize() == rv._dsize()), "wrong rvn");
	1258	//Permutation vector of the new R
	1259	ivec irvn = rvn.dataind (rv);
	1260	ivec irvc = rvc.dataind (rv);
	1261	ivec perm = concat (irvn , irvc);
	1262	sq_T Rn (R, perm);
	1263
	1264	//fixme - could this be done in general for all sq_T?
	1265	mat S = Rn.to_mat();
	1266	//fixme
	1267	int n = rvn._dsize() - 1;
	1268	int end = R.rows() - 1;
	1269	mat S11 = S.get (0, n, 0, n);
	1270	mat S12 = S.get (0, n , rvn._dsize(), end);
	1271	mat S22 = S.get (rvn._dsize(), end, rvn._dsize(), end);
	1272
	1273	vec mu1 = mu (irvn);
	1274	vec mu2 = mu (irvc);
	1275	mat A = S12 * inv (S22);
	1276	sq_T R_n (S11 - A *S12.T());
	1277
	1278	mlnorm<sq_T>* tmp = new mlnorm<sq_T> ();
	1279	tmp->set_rv (rvn); tmp->set_rvc (rvc);
	1280	tmp->set_parameters (A, mu1 - A*mu2, R_n);
	1281	return tmp;
	1282	}
	1283
	1284	////
	1285	///////
	1286	template<class sq_T>
	1287	void mgnorm<sq_T >::set_parameters (fnc* g0, const sq_T &R0) {g = g0; this->iepdf.set_parameters (zeros (g->dimension()), R0);}
	1288	template<class sq_T>
	1289	void mgnorm<sq_T >::condition (const vec &cond) {this->iepdf._mu() = g->eval (cond);};
	1290
	1291	template<class sq_T>
	1292	std::ostream &operator<< (std::ostream &os, mlnorm<sq_T> &ml)
	1293	{
	1294	os << "A:" << ml.A << endl;
	1295	os << "mu:" << ml.mu_const << endl;
	1296	os << "R:" << ml._R() << endl;
	1297	return os;
	1298	};
1332	1299
1333	1300	}

Context Navigation

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

Legend:

library/bdm/stat/exp_family.h

Download in other formats: