Context Navigation

exp_family.h @ 487

Revision 487, 34.5 kB (checked in by smidl, 15 years ago)
1st step of mpdf redesign - BROKEN compile
Property svn:eol-style set to `native`

Line
1	/*!
2	\file
3	\brief Probability distributions for Exponential Family models.
4	\author Vaclav Smidl.
5
6	-----------------------------------
7	BDM++ - C++ library for Bayesian Decision Making under Uncertainty
8
9	Using IT++ for numerical operations
10	-----------------------------------
11	*/
12
13	#ifndef EF_H
14	#define EF_H
15
16
17	#include "../shared_ptr.h"
18	#include "../base/bdmbase.h"
19	#include "../math/chmat.h"
20
21	namespace bdm
22	{
23
24
25	//! Global Uniform_RNG
26	extern Uniform_RNG UniRNG;
27	//! Global Normal_RNG
28	extern Normal_RNG NorRNG;
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:
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();
62	}
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();
122	}
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;
196	}
197	}
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;
232	}
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:
435	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	}
1168	}
1169	}
1170	};
1171
1172
1173	////////////////////////
1174
1175	template<class sq_T>
1176	void enorm<sq_T>::set_parameters ( const vec &mu0, const sq_T &R0 )
1177	{
1178	//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	};
1200
1201	// template<class sq_T>
1202	// void enorm<sq_T>::tupdate ( double phi, mat &vbar, double nubar ) {
1203	// //
1204	// };
1205
1206	template<class sq_T>
1207	vec enorm<sq_T>::sample() const
1208	{
1209	vec x ( dim );
1210	#pragma omp critical
1211	NorRNG.sample_vector ( dim,x );
1212	vec smp = R.sqrt_mult ( x );
1213
1214	smp += mu;
1215	return smp;
1216	};
1217
1218	// template<class sq_T>
1219	// double enorm<sq_T>::eval ( const vec &val ) const {
1220	// double pdfl,e;
1221	// pdfl = evallog ( val );
1222	// e = exp ( pdfl );
1223	// return e;
1224	// };
1225
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	};
1241
1242
1243	// template<class sq_T>
1244	// vec mlnorm<sq_T>::samplecond (const vec &cond, double &lik ) {
1245	// this->condition ( cond );
1246	// vec smp = epdf.sample();
1247	// lik = epdf.eval ( smp );
1248	// return smp;
1249	// }
1250
1251	// template<class sq_T>
1252	// mat mlnorm<sq_T>::samplecond (const vec &cond, vec &lik, int n ) {
1253	// int i;
1254	// int dim = rv.count();
1255	// mat Smp ( dim,n );
1256	// vec smp ( dim );
1257	// this->condition ( cond );
1258	//
1259	// for ( i=0; i<n; i++ ) {
1260	// smp = epdf.sample();
1261	// lik ( i ) = epdf.eval ( smp );
1262	// Smp.set_col ( i ,smp );
1263	// }
1264	//
1265	// return Smp;
1266	// }
1267
1268
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	};
1332
1333	}
1334	#endif //EF_H

Note: See TracBrowser for help on using the browser.

Context Navigation

root/library/bdm/stat/exp_family.h @ 487

Download in other formats: