Context Navigation

← Previous Change
Next Change →

Changeset 145 for bdm

Timestamp:

08/20/08 15:41:21 (16 years ago)

Author:

smidl

Message:

Oprava dokumentace

Location:

bdm

Files:

: 8 modified

estim/arx.h (modified) (1 diff)
itpp_ext.cpp (modified) (3 diffs)
itpp_ext.h (modified) (2 diffs)
stat/emix.cpp (modified) (1 diff)
stat/emix.h (modified) (4 diffs)
stat/libBM.cpp (modified) (4 diffs)
stat/libBM.h (modified) (7 diffs)
stat/loggers.h (modified) (1 diff)

Legend:

: Unmodified
: Added
: Removed

bdm/estim/arx.h

r135	r145
51	51	//! Full constructor
52	52	ARX (RV &rv, mat &V0, double &nu0, double frg0=1.0) : BM(rv),est(rv,V0,nu0), V(est._V()), nu(est._nu()), frg(frg0){last_lognc=est.lognc();tll=0.0;};
	53	//! Set sufficient statistics
53	54	void set_parameters(mat &V0, double &nu0){est._V()=V0;est._nu()=nu0;last_lognc=est.lognc();tll=last_lognc;}
	55	//! Returns sufficient statistics
54	56	void get_parameters(mat &V0, double &nu0){V0=est._V().to_mat(); nu0=est._nu();}
55	57	//! Here \f$dt = [y_t psi_t] \f$.

bdm/itpp_ext.cpp

r141	r145
2	2	// C++ Implementation: itpp_ext
3	3	//
4		// Description:
	4	// Description:
5	5	//
6	6	//
…	…
16	16
17	17	extern "C" { /* QR factorization of a general matrix A */
18		~~void dgeqrf_(~~int m, int n, double a, int lda, double tau, double work,
19		~~int lwork, int info~~);
	18	void dgeqrf_ ( int m, int n, double a, int lda, double tau, double work,
	19	int lwork, int info );
20	20	};
21	21
22	22	namespace itpp {
23		Array<int> to_Arr(const ivec &indices)
24		{
25		Array<int> a(indices.size());
26		for (int i = 0; i < a.size(); i++) {
27		a(i) = indices(i);
28		}
29		return a;
30		}
31
32		ivec linspace ( int from, int to ) {
33		int n=to-from+1;
34		int i;
35		it_assert_debug ( n>0,"wrong linspace" );
36		ivec iv ( n ); for (i=0;i<n;i++) iv(i)=from+i;
37		return iv;
38		};
	23	Array<int> to_Arr ( const ivec &indices ) {
	24	Array<int> a ( indices.size() );
	25	for ( int i = 0; i < a.size(); i++ ) {
	26	a ( i ) = indices ( i );
	27	}
	28	return a;
	29	}
	30
	31	ivec linspace ( int from, int to ) {
	32	int n=to-from+1;
	33	int i;
	34	it_assert_debug ( n>0,"wrong linspace" );
	35	ivec iv ( n ); for ( i=0;i<n;i++ ) iv ( i ) =from+i;
	36	return iv;
	37	};
	38
	39	void set_subvector ( vec &ov, const ivec &iv, const vec &v )
	40	{
	41	it_assert_debug ( ( iv.length() <=v.length() ),
	42	"Vec<>::set_subvector(ivec, vec<Num_T>): Indexing out "
	43	"of range of v" );
	44	for ( int i = 0; i < iv.length(); i++ ) {
	45	it_assert_debug ( iv ( i ) <ov.length(),
	46	"Vec<>::set_subvector(ivec, vec<Num_T>): Indexing out "
	47	"of range of v" );
	48	ov ( iv ( i ) ) = v ( i );
	49	}
	50	}
39	51
40	52	//Gamma
41	53
42		Gamma_RNG::Gamma_RNG(double a, double b)
43		{
44		setup(a,b);
45		}
46
	54	Gamma_RNG::Gamma_RNG ( double a, double b ) {
	55	setup ( a,b );
	56	}
	57
	58
	59	bvec operator& ( const bvec &a, const bvec &b ) {
	60	it_assert_debug ( b.size() ==a.size(), "operator&(): Vectors of different lengths" );
	61
	62	bvec temp ( a.size() );
	63	for ( int i = 0;i < a.size();i++ ) {
	64	temp ( i ) = a ( i ) & b ( i );
	65	}
	66	return temp;
	67	}
	68
	69	bvec operator\| ( const bvec &a, const bvec &b ) {
	70	it_assert_debug ( b.size() !=a.size(), "operator&(): Vectors of different lengths" );
	71
	72	bvec temp ( a.size() );
	73	for ( int i = 0;i < a.size();i++ ) {
	74	temp ( i ) = a ( i ) \| b ( i );
	75	}
	76	return temp;
	77	}
47	78	#define log std::log
48	79	#define exp std::exp
…	…
50	81	#define R_FINITE std::isfinite
51	82
52		double Gamma_RNG::sample()
53		{
54		//A copy of rgamma code from the R package!!
55		//
56
57		/* Constants : */
58		const static double sqrt32 = 5.656854;
59		const static double exp_m1 = 0.36787944117144232159;/* exp(-1) = 1/e */
60
61		/* Coefficients q[k] - for q0 = sum(q[k]*a^(-k))
62		* Coefficients a[k] - for q = q0+(tt/2)sum(a[k]*v^k)
63		* Coefficients e[k] - for exp(q)-1 = sum(e[k]*q^k)
64		*/
65		const static double q1 = 0.04166669;
66		const static double q2 = 0.02083148;
67		const static double q3 = 0.00801191;
68		const static double q4 = 0.00144121;
69		const static double q5 = -7.388e-5;
70		const static double q6 = 2.4511e-4;
71		const static double q7 = 2.424e-4;
72
73		const static double a1 = 0.3333333;
74		const static double a2 = -0.250003;
75		const static double a3 = 0.2000062;
76		const static double a4 = -0.1662921;
77		const static double a5 = 0.1423657;
78		const static double a6 = -0.1367177;
79		const static double a7 = 0.1233795;
80
81		/* State variables [FIXME for threading!] :*/
82		static double aa = 0.;
83		static double aaa = 0.;
84		static double s, s2, d; /* no. 1 (step 1) */
85		static double q0, b, si, c;/* no. 2 (step 4) */
86
87		double e, p, q, r, t, u, v, w, x, ret_val;
88		double a=alpha;
89		double scale=1.0/beta;
90
91		if (!R_FINITE(a) \|\| !R_FINITE(scale) \|\| a < 0.0 \|\| scale <= 0.0)
92		{it_error("Gamma_RNG wrong parameters");}
93
94		if (a < 1.) { /* GS algorithm for parameters a < 1 */
95		if(a == 0)
96		return 0.;
97		e = 1.0 + exp_m1 * a;
98		for(;;) { //VS repeat
99		p = e * unif_rand();
100		if (p >= 1.0) {
101		x = -log((e - p) / a);
102		if (exp_rand() >= (1.0 - a) * log(x))
103		break;
104		} else {
105		x = exp(log(p) / a);
106		if (exp_rand() >= x)
107		break;
108		}
109		}
110		return scale * x;
111		}
112
113		/* --- a >= 1 : GD algorithm --- */
114
115		/* Step 1: Recalculations of s2, s, d if a has changed */
116		if (a != aa) {
117		aa = a;
118		s2 = a - 0.5;
119		s = sqrt(s2);
120		d = sqrt32 - s * 12.0;
121		}
122		/* Step 2: t = standard normal deviate,
123		x = (s,1/2) -normal deviate. */
124
125		/* immediate acceptance (i) */
126		t = norm_rand();
127		x = s + 0.5 * t;
128		ret_val = x * x;
129		if (t >= 0.0)
130		return scale * ret_val;
131
132		/* Step 3: u = 0,1 - uniform sample. squeeze acceptance (s) */
133		u = unif_rand();
134		if ((d * u) <= (t * t * t))
135		return scale * ret_val;
136
137		/* Step 4: recalculations of q0, b, si, c if necessary */
138
139		if (a != aaa) {
140		aaa = a;
141		r = 1.0 / a;
142		q0 = ((((((q7 * r + q6) * r + q5) * r + q4) * r + q3) * r
143		+ q2) * r + q1) * r;
144
145		/* Approximation depending on size of parameter a */
146		/* The constants in the expressions for b, si and c */
147		/* were established by numerical experiments */
148
149		if (a <= 3.686) {
150		b = 0.463 + s + 0.178 * s2;
151		si = 1.235;
152		c = 0.195 / s - 0.079 + 0.16 * s;
153		} else if (a <= 13.022) {
154		b = 1.654 + 0.0076 * s2;
155		si = 1.68 / s + 0.275;
156		c = 0.062 / s + 0.024;
157		} else {
158		b = 1.77;
159		si = 0.75;
160		c = 0.1515 / s;
161		}
162		}
163		/* Step 5: no quotient test if x not positive */
164
165		if (x > 0.0) {
166		/* Step 6: calculation of v and quotient q */
167		v = t / (s + s);
168		if (fabs(v) <= 0.25)
169		q = q0 + 0.5 * t * t * ((((((a7 * v + a6) * v + a5) * v + a4) * v
170		+ a3) * v + a2) * v + a1) * v;
171		else
172		q = q0 - s * t + 0.25 * t * t + (s2 + s2) * log(1.0 + v);
173
174
175		/* Step 7: quotient acceptance (q) */
176		if (log(1.0 - u) <= q)
177		return scale * ret_val;
178		}
179
180		for(;;) { //VS repeat
181		/* Step 8: e = standard exponential deviate
182		* u = 0,1 -uniform deviate
183		* t = (b,si)-double exponential (laplace) sample */
184		e = exp_rand();
185		u = unif_rand();
186		u = u + u - 1.0;
187		if (u < 0.0)
188		t = b - si * e;
189		else
190		t = b + si * e;
191		/* Step 9: rejection if t < tau(1) = -0.71874483771719 */
192		if (t >= -0.71874483771719) {
193		/* Step 10: calculation of v and quotient q */
194		v = t / (s + s);
195		if (fabs(v) <= 0.25)
196		q = q0 + 0.5 * t * t *
197		((((((a7 * v + a6) * v + a5) * v + a4) * v + a3) * v
198		+ a2) * v + a1) * v;
199		else
200		q = q0 - s * t + 0.25 * t * t + (s2 + s2) * log(1.0 + v);
201		/* Step 11: hat acceptance (h) */
202		/* (if q not positive go to step 8) */
203		if (q > 0.0) {
204		// TODO: w = expm1(q);
205		w = exp(q)-1;
206		/* ^^^^^ original code had approximation with rel.err < 2e-7 */
207		/* if t is rejected sample again at step 8 */
208		if ((c * fabs(u)) <= (w * exp(e - 0.5 * t * t)))
209		break;
210		}
211		}
212		} /* repeat .. until `t' is accepted */
213		x = s + 0.5 * t;
214		return scale * x * x;
	83	double Gamma_RNG::sample() {
	84	//A copy of rgamma code from the R package!!
	85	//
	86
	87	/* Constants : */
	88	const static double sqrt32 = 5.656854;
	89	const static double exp_m1 = 0.36787944117144232159;/* exp(-1) = 1/e */
	90
	91	/* Coefficients q[k] - for q0 = sum(q[k]*a^(-k))
	92	* Coefficients a[k] - for q = q0+(tt/2)sum(a[k]*v^k)
	93	* Coefficients e[k] - for exp(q)-1 = sum(e[k]*q^k)
	94	*/
	95	const static double q1 = 0.04166669;
	96	const static double q2 = 0.02083148;
	97	const static double q3 = 0.00801191;
	98	const static double q4 = 0.00144121;
	99	const static double q5 = -7.388e-5;
	100	const static double q6 = 2.4511e-4;
	101	const static double q7 = 2.424e-4;
	102
	103	const static double a1 = 0.3333333;
	104	const static double a2 = -0.250003;
	105	const static double a3 = 0.2000062;
	106	const static double a4 = -0.1662921;
	107	const static double a5 = 0.1423657;
	108	const static double a6 = -0.1367177;
	109	const static double a7 = 0.1233795;
	110
	111	/* State variables [FIXME for threading!] :*/
	112	static double aa = 0.;
	113	static double aaa = 0.;
	114	static double s, s2, d; /* no. 1 (step 1) */
	115	static double q0, b, si, c;/* no. 2 (step 4) */
	116
	117	double e, p, q, r, t, u, v, w, x, ret_val;
	118	double a=alpha;
	119	double scale=1.0/beta;
	120
	121	if ( !R_FINITE ( a ) \|\| !R_FINITE ( scale ) \|\| a < 0.0 \|\| scale <= 0.0 )
	122	{it_error ( "Gamma_RNG wrong parameters" );}
	123
	124	if ( a < 1. ) { /* GS algorithm for parameters a < 1 */
	125	if ( a == 0 )
	126	return 0.;
	127	e = 1.0 + exp_m1 * a;
	128	for ( ;; ) { //VS repeat
	129	p = e * unif_rand();
	130	if ( p >= 1.0 ) {
	131	x = -log ( ( e - p ) / a );
	132	if ( exp_rand() >= ( 1.0 - a ) * log ( x ) )
	133	break;
	134	}
	135	else {
	136	x = exp ( log ( p ) / a );
	137	if ( exp_rand() >= x )
	138	break;
	139	}
	140	}
	141	return scale * x;
	142	}
	143
	144	/* --- a >= 1 : GD algorithm --- */
	145
	146	/* Step 1: Recalculations of s2, s, d if a has changed */
	147	if ( a != aa ) {
	148	aa = a;
	149	s2 = a - 0.5;
	150	s = sqrt ( s2 );
	151	d = sqrt32 - s * 12.0;
	152	}
	153	/* Step 2: t = standard normal deviate,
	154	x = (s,1/2) -normal deviate. */
	155
	156	/* immediate acceptance (i) */
	157	t = norm_rand();
	158	x = s + 0.5 * t;
	159	ret_val = x * x;
	160	if ( t >= 0.0 )
	161	return scale * ret_val;
	162
	163	/* Step 3: u = 0,1 - uniform sample. squeeze acceptance (s) */
	164	u = unif_rand();
	165	if ( ( d * u ) <= ( t * t * t ) )
	166	return scale * ret_val;
	167
	168	/* Step 4: recalculations of q0, b, si, c if necessary */
	169
	170	if ( a != aaa ) {
	171	aaa = a;
	172	r = 1.0 / a;
	173	q0 = ( ( ( ( ( ( q7 * r + q6 ) * r + q5 ) * r + q4 ) * r + q3 ) * r
	174	+ q2 ) * r + q1 ) * r;
	175
	176	/* Approximation depending on size of parameter a */
	177	/* The constants in the expressions for b, si and c */
	178	/* were established by numerical experiments */
	179
	180	if ( a <= 3.686 ) {
	181	b = 0.463 + s + 0.178 * s2;
	182	si = 1.235;
	183	c = 0.195 / s - 0.079 + 0.16 * s;
	184	}
	185	else if ( a <= 13.022 ) {
	186	b = 1.654 + 0.0076 * s2;
	187	si = 1.68 / s + 0.275;
	188	c = 0.062 / s + 0.024;
	189	}
	190	else {
	191	b = 1.77;
	192	si = 0.75;
	193	c = 0.1515 / s;
	194	}
	195	}
	196	/* Step 5: no quotient test if x not positive */
	197
	198	if ( x > 0.0 ) {
	199	/* Step 6: calculation of v and quotient q */
	200	v = t / ( s + s );
	201	if ( fabs ( v ) <= 0.25 )
	202	q = q0 + 0.5 * t * t * ( ( ( ( ( ( a7 * v + a6 ) * v + a5 ) * v + a4 ) * v
	203	+ a3 ) * v + a2 ) * v + a1 ) * v;
	204	else
	205	q = q0 - s * t + 0.25 * t * t + ( s2 + s2 ) * log ( 1.0 + v );
	206
	207
	208	/* Step 7: quotient acceptance (q) */
	209	if ( log ( 1.0 - u ) <= q )
	210	return scale * ret_val;
	211	}
	212
	213	for ( ;; ) { //VS repeat
	214	/* Step 8: e = standard exponential deviate
	215	* u = 0,1 -uniform deviate
	216	* t = (b,si)-double exponential (laplace) sample */
	217	e = exp_rand();
	218	u = unif_rand();
	219	u = u + u - 1.0;
	220	if ( u < 0.0 )
	221	t = b - si * e;
	222	else
	223	t = b + si * e;
	224	/* Step 9: rejection if t < tau(1) = -0.71874483771719 */
	225	if ( t >= -0.71874483771719 ) {
	226	/* Step 10: calculation of v and quotient q */
	227	v = t / ( s + s );
	228	if ( fabs ( v ) <= 0.25 )
	229	q = q0 + 0.5 * t * t *
	230	( ( ( ( ( ( a7 * v + a6 ) * v + a5 ) * v + a4 ) * v + a3 ) * v
	231	+ a2 ) * v + a1 ) * v;
	232	else
	233	q = q0 - s * t + 0.25 * t * t + ( s2 + s2 ) * log ( 1.0 + v );
	234	/* Step 11: hat acceptance (h) */
	235	/* (if q not positive go to step 8) */
	236	if ( q > 0.0 ) {
	237	// TODO: w = expm1(q);
	238	w = exp ( q )-1;
	239	/* ^^^^^ original code had approximation with rel.err < 2e-7 */
	240	/* if t is rejected sample again at step 8 */
	241	if ( ( c * fabs ( u ) ) <= ( w * exp ( e - 0.5 * t * t ) ) )
	242	break;
	243	}
	244	}
	245	} /* repeat .. until `t' is accepted */
	246	x = s + 0.5 * t;
	247	return scale * x * x;
	248	}
	249
	250
	251	bool qr ( const mat &A, mat &R ) {
	252	int info;
	253	int m = A.rows();
	254	int n = A.cols();
	255	int lwork = n;
	256	int k = std::min ( m, n );
	257	vec tau ( k );
	258	vec work ( lwork );
	259
	260	R = A;
	261
	262	// perform workspace query for optimum lwork value
	263	int lwork_tmp = -1;
	264	dgeqrf_ ( &m, &n, R._data(), &m, tau._data(), work._data(), &lwork_tmp,
	265	&info );
	266	if ( info == 0 ) {
	267	lwork = static_cast<int> ( work ( 0 ) );
	268	work.set_size ( lwork, false );
	269	}
	270	dgeqrf_ ( &m, &n, R._data(), &m, tau._data(), work._data(), &lwork, &info );
	271
	272	// construct R
	273	for ( int i=0; i<m; i++ )
	274	for ( int j=0; j<std::min ( i,n ); j++ )
	275	R ( i,j ) = 0;
	276
	277	return ( info==0 );
	278	}
	279
215	280	}
216
217
218		~~bool qr(const mat &A, mat &R)~~
219		{
220		~~int info;~~
221		~~int m = A.rows();~~
222		~~int n = A.cols();~~
223		~~int lwork = n;~~
224		~~int k = std::min(m, n);~~
225		~~vec tau(k);~~
226		~~vec work(lwork);~~
227
228		~~R = A;~~
229
230		~~// perform workspace query for optimum lwork value~~
231		~~int lwork_tmp = -1;~~
232		~~dgeqrf_(&m, &n, R._data(), &m, tau._data(), work._data(), &lwork_tmp,~~
233		~~&info);~~
234		~~if (info == 0) {~~
235		~~lwork = static_cast<int>(work(0));~~
236		~~work.set_size(lwork, false);~~
237		}
238		~~dgeqrf_(&m, &n, R._data(), &m, tau._data(), work._data(), &lwork, &info);~~
239
240		~~// construct R~~
241		~~for (int i=0; i<m; i++)~~
242		~~for (int j=0; j<std::min(i,n); j++)~~
243		~~R(i,j) = 0;~~
244
245		~~return (info==0);~~
246		}
247
248		}

bdm/itpp_ext.h

r91	r145
13	13
14	14	#ifndef ITEX_H
15		#define ITEX_H
	15	#define ITEX_H
16	16
17	17	using std::cout;
…	…
19	19
20	20	namespace itpp {
21		Array<int> to_Arr ( const ivec &indices );
22		ivec linspace ( int from, int to );
	21	Array<int> to_Arr ( const ivec &indices );
	22	ivec linspace ( int from, int to );
23	23
24		/*!
25		\brief Gamma distribution
26		\ingroup randgen
27		*/
28		class Gamma_RNG {
29		public:
	24	bvec operator& ( const bvec &a, const bvec &b );
	25	bvec operator\| ( const bvec &a, const bvec &b );
	26
	27	// template<class Num_T>
	28	// void set_subvector(vec &ov, ivec &iv, const Vec<Num_T> &v);
	29
	30	void set_subvector ( vec &ov, const ivec &iv, const vec &v );
	31
	32	/*!
	33	\brief Gamma distribution
	34	\ingroup randgen
	35	*/
	36	class Gamma_RNG {
	37	public:
30	38	//! constructor. Set lambda.
31		Gamma_RNG ( double a=1.0, double b=1.0 );
32		//! Set lambda
33		void setup ( double a0, double b0 ) { alpha=a0; beta=b0;}
34		//! get lambda
35		double get_setup() const;
36		//! Get one sample.
37		double operator() () { return sample(); }
38		//! Get a sample vector.
39		vec operator() ( int n );
40		//! Get a sample matrix.
41		mat operator() ( int h, int w );
42		protected:
43		private:
44		//!
45		double sample();
46		//!
47		double alpha;
48		//!
49		double beta;
50		//!
51		Random_Generator RNG;
52		Normal_RNG NRNG;
53		//! compatibility with R
54		inline double exp_rand() {return -std::log ( RNG.random_01() );}
55		inline double unif_rand() {return RNG.random_01();}
56		inline double norm_rand() {return NRNG.sample();}
	39	Gamma_RNG ( double a=1.0, double b=1.0 );
	40	//! Set lambda
	41	void setup ( double a0, double b0 ) { alpha=a0; beta=b0;}
	42	//! get lambda
	43	double get_setup() const;
	44	//! Get one sample.
	45	double operator() () { return sample(); }
	46	//! Get a sample vector.
	47	vec operator() ( int n );
	48	//! Get a sample matrix.
	49	mat operator() ( int h, int w );
	50	protected:
	51	private:
	52	//!
	53	double sample();
	54	//!
	55	double alpha;
	56	//!
	57	double beta;
	58	//!
	59	Random_Generator RNG;
	60	Normal_RNG NRNG;
	61	//! compatibility with R
	62	inline double exp_rand() {return -std::log ( RNG.random_01() );}
	63	inline double unif_rand() {return RNG.random_01();}
	64	inline double norm_rand() {return NRNG.sample();}
57	65
58		};
	66	};
59	67
60		bool qr ( const mat &A, mat &R );
	68	bool qr ( const mat &A, mat &R );
61	69
62	70	}
	71
	72
63	73	#endif //ITEX_H

bdm/stat/emix.cpp

r141	r145
18	18
19	19	int i=0;
20		while ( ( w ( i ) <u0 ) && ( i< ( w.length()-1 ) ) ) {i++;}
	20	while ( ( cumDist ( i ) <u0 ) && ( i< ( w.length()-1 ) ) ) {i++;}
21	21
22	22	return Coms ( i )->sample();

bdm/stat/emix.h

r141	r145
30	30
31	31	*/
32		class emix : public epdf
33		{
	32	class emix : public epdf {
34	33	protected:
35	34	//! weights of the components
…	…
39	38	public:
40	39	//!Default constructor
41		emix ~~( RV &rv ) : epdf ( rv~~ ) {};
	40	emix(RV &rv) : epdf(rv) {};
42	41	//! Set weights \c w and components \c R
43		void set_parameters ~~( const vec &w, const Array<epdf*> &Coms~~ );
	42	void set_parameters(const vec &w, const Array<epdf*> &Coms);
44	43
45	44	vec sample() const;
46		vec mean() const
47		{
48		int i; vec mu=zeros ( rv.count() );
49		for ( i=0;i<w.length();i++ ) {mu+=w ( i ) *Coms ( i )->mean(); }
	45	vec mean() const {
	46	int i; vec mu = zeros(rv.count());
	47	for (i = 0;i < w.length();i++) {mu += w(i) * Coms(i)->mean(); }
50	48	return mu;
51	49	}
52		double evalpdflog ( const vec &val ) const
53		{
	50	double evalpdflog(const vec &val) const {
54	51	int i;
55		double sum=0.0;
56		for ( ~~i=0;i<w.length();i++ ) {sum+=w ( i ) *Coms ( i )->evalpdflog ( val~~ );}
57		return log ~~( sum~~ );
	52	double sum = 0.0;
	53	for (i = 0;i < w.length();i++) {sum += w(i) * Coms(i)->evalpdflog(val);}
	54	return log(sum);
58	55	};
59	56
…	…
65	62	/*! \brief Chain rule decomposition of epdf
66	63
	64	Probability density in the form of Chain-rule decomposition:
	65	\[
	66	f(x_1,x_2,x_3) = f(x_1\|x_2,x_3)f(x_2,x_3)f(x_3)
	67	\]
	68	Note that
	69	*/
	70	class eprod: public epdf {
	71	protected:
	72	//
	73	int n;
	74	// pointers to epdfs
	75	Array<epdf*> epdfs;
	76	Array<mpdf*> mpdfs;
	77	//
	78	Array<ivec> rvinds;
	79	Array<ivec> rvcinds;
	80	//! Indicate independence of its factors
	81	bool independent;
	82	//! Indicate internal creation of mpdfs which must be destroyed
	83	bool intermpdfs;
	84	public:
	85	//!Constructor from list of eFacs or list of mFacs
	86	eprod(Array<mpdf*> mFacs): epdf(RV()), n(mFacs.length()), epdfs(n), mpdfs(mFacs), rvinds(n), rvcinds(n) {
	87	int i;
	88	intermpdfs = false;
	89	for (i = 0;i < n;i++) {
	90	rv.add(mpdfs(i)->_rv()); //add rv to common rvs.
	91	epdfs(i) = &(mpdfs(i)->_epdf()); // add pointer to epdf
	92	};
	93	independent = true;
	94	//test rvc of mpdfs and fill rvinds
	95	for (i = 0;i < n;i++) {
	96	// find ith rv in common rv
	97	rvinds(i) = mpdfs(i)->_rv().dataind(rv);
	98	// find ith rvc in common rv
	99	rvcinds(i) = mpdfs(i)->_rvc().dataind(rv);
	100	if (rvcinds(i).length()>0) {independent = false;}
	101	}
	102
	103	};
	104	eprod(Array<epdf*> eFacs): epdf(RV()), n(eFacs.length()), epdfs(eFacs), mpdfs(n), rvinds(n), rvcinds(n) {
	105	int i;
	106	intermpdfs = true;
	107	for (i = 0;i < n;i++) {
	108	if (rv.add(eFacs(i)->_rv())) {it_error("Incompatible eFacs.rv!");} //add rv to common rvs.
	109	mpdfs(i) = new mepdf(*(epdfs(i)));
	110	epdfs(i) = &(mpdfs(i)->_epdf()); // add pointer to epdf
	111	};
	112	independent = true;
	113	//test rvc of mpdfs and fill rvinds
	114	for (i = 0;i < n;i++) {
	115	// find ith rv in common rv
	116	rvinds(i) = mpdfs(i)->_rv().dataind(rv);
	117	}
	118	};
	119
	120	double evalpdflog(const vec &val) const {
	121	int i;
	122	double res = 0.0;
	123	for (i = n - 1;i > 0;i++) {
	124	if (rvcinds(i).length() > 0)
	125	{mpdfs(i)->condition(val(rvcinds(i)));}
	126	// add logarithms
	127	res += epdfs(i)->evalpdflog(val(rvinds(i)));
	128	}
	129	}
	130	vec sample () const{
	131	vec smp=zeros(rv.count());
	132	for (int i = (n - 1);i >= 0;i--) {
	133	if (rvcinds(i).length() > 0)
	134	{mpdfs(i)->condition(smp(rvcinds(i)));}
	135	set_subvector(smp,rvinds(i), epdfs(i)->sample());
	136	}
	137	return smp;
	138	}
	139	vec mean() const{
	140	vec tmp(rv.count());
	141	if (independent) {
	142	for (int i=0;i<n;i++) {
	143	vec pom = epdfs(i)->mean();
	144	set_subvector(tmp,rvinds(i), pom);
	145	}
	146	}
	147	else {
	148	int N=50*rv.count();
	149	it_warning("eprod.mean() computed by sampling");
	150	tmp = zeros(rv.count());
	151	for (int i=0;i<N;i++){ tmp += sample();}
	152	tmp /=N;
	153	}
	154	return tmp;
	155	}
	156	~eprod(){if (intermpdfs) {for (int i=0;i<n;i++){delete mpdfs(i);}}};
	157	};
	158
	159	/*! \brief Mixture of mpdfs with constant weights, all mpdfs are of equal type
67	160
68	161	*/
69		class eprod: public epdf
70		{
71		protected:
72		Array<epdf*> epdfs;
73		Array<mpdf*> mpdfs;
74		public:
75
76
77		};
78
79		/*! \brief Mixture of mpdfs with constant weights
80
81		*/
82		class mmix : public mpdf
83		{
	162	class mmix : public mpdf {
84	163	protected:
85	164	//! Component (epdfs)
…	…
89	168	public:
90	169	//!Default constructor
91		mmix ~~( RV &rv, RV &rvc ) : mpdf ( rv, rvc ), Epdf ( rv ) {ep=~~&Epdf;};
	170	mmix(RV &rv, RV &rvc) : mpdf(rv, rvc), Epdf(rv) {ep = &Epdf;};
92	171	//! Set weights \c w and components \c R
93		void set_parameters ( const vec &w, const Array<mpdf*> &Coms )
94		{
95		Array<epdf*> Eps ( Coms.length() );
	172	void set_parameters(const vec &w, const Array<mpdf*> &Coms) {
	173	Array<epdf*> Eps(Coms.length());
96	174
97		for ( int i=0;i<Coms.length();i++ )
98		{
99		Eps ( i ) =& ( Coms ( i )->_epdf() );
	175	for (int i = 0;i < Coms.length();i++) {
	176	Eps(i) = & (Coms(i)->_epdf());
100	177	}
101		Epdf.set_parameters ~~( w,Eps~~ );
	178	Epdf.set_parameters(w, Eps);
102	179	};
103	180
104		void condition ( const vec &cond )
105		{
106		for ( int i=0;i<Coms.length();i++ ) {Coms ( i )->condition ( cond );}
	181	void condition(const vec &cond) {
	182	for (int i = 0;i < Coms.length();i++) {Coms(i)->condition(cond);}
107	183	};
108	184	};

bdm/stat/libBM.cpp

r102	r145
25	25	times = in_times;
26	26	tsize = 0;
27		for ( i~~=0;i<len;i++ ) {tsize+=~~sizes ( i );}
	27	for ( i = 0;i < len;i++ ) {tsize += sizes ( i );}
28	28	};
29	29
…	…
32	32	}
33	33
34		RV::RV ~~() : tsize ( 0 ),~~len ( 0 ) {};
	34	RV::RV() : tsize ( 0 ), len ( 0 ) {};
35	35
36		~~void~~ RV::add ( const RV &rv2 ) {
	36	bool RV::add ( const RV &rv2 ) {
37	37	// TODO
38		tsize+=rv2.tsize;
39		len +=rv2.len;
40		ids=concat ( ids,rv2.ids );
41		sizes=concat ( sizes,rv2.sizes );
42		times=concat ( times,rv2.times );
43		names=concat ( names,rv2.names );
	38	ivec ind = rv2.findself ( *this ); //should be -1 all the time
	39	ivec index = itpp::find(ind==-1);
	40
	41	if ( index.length() < rv2.len ) { //conflict
	42	ids = concat ( ids, rv2.ids(index) );
	43	sizes = concat ( sizes, rv2.sizes(index) );
	44	times = concat ( times, rv2.times(index) );
	45	names = concat ( names, rv2.names(to_Arr(index)) );
	46	}
	47	else {
	48	ids = concat ( ids, rv2.ids );
	49	sizes = concat ( sizes, rv2.sizes );
	50	times = concat ( times, rv2.times );
	51	names = concat ( names, rv2.names );
	52	}
	53	tsize = sum(sizes);
	54	len = ids.length();
	55	return (index.length()<rv2.len);
	56
44	57	// return *this;
45	58	};
…	…
58	71	}
59	72
60		RV RV::subselect ( ivec ind ) {
	73	RV RV::subselect ( ivec ind ) const {
61	74	return RV ( ids ( ind ), names ( to_Arr ( ind ) ), sizes ( ind ), times ( ind ) );
62	75	}
63	76
64		void RV::t ( int delta ) { times +=delta;}
	77	void RV::t ( int delta ) { times += delta;}
65	78
66		RV RV::operator() ( ivec ind ) {
	79	RV RV::operator() ( ivec ind ) const {
67	80	return RV ( ids ( ind ), names ( to_Arr ( ind ) ), sizes ( ind ), times ( ind ) );
68	81	}
69	82
70		bool RV::equal ( RV rv2 ) const {
71		return ( ids~~==rv2.ids ) && ( times == rv2.times ) && ( sizes==~~rv2.sizes );
	83	bool RV::equal ( const RV &rv2 ) const {
	84	return ( ids == rv2.ids ) && ( times == rv2.times ) && ( sizes == rv2.sizes );
72	85	}
73	86
74	87	mat epdf::sampleN ( int N ) const {
75		mat X =~~zeros ( rv.count(),~~N );
76		for ( int i~~=0;i<N;i++ ) X.set_col ( i,this->sample() );~~
	88	mat X = zeros ( rv.count(), N );
	89	for ( int i = 0;i < N;i++ ) X.set_col ( i, this->sample() );
77	90	return X;
78	91	};
…	…
88	101	}
89	102
90		ivec RV::indexlist() {
91		ivec indlist ( tsize );
	103	str RV::tostr() const {
	104	ivec idlist ( tsize );
	105	ivec tmlist ( tsize );
92	106	int i;
93	107	int pos = 0;
94		for ( i=0;i<len;i++ ) {
95		indlist.set_subvector ( pos,pos+sizes ( i )-1, ids ( i ) );
	108	for ( i = 0;i < len;i++ ) {
	109	idlist.set_subvector ( pos, pos + sizes ( i ) - 1, ids ( i ) );
	110	tmlist.set_subvector ( pos, pos + sizes ( i ) - 1, times ( i ) );
	111	pos += sizes ( i );
96	112	}
97		return indlist;
	113	return str ( idlist, tmlist );
	114	}
	115
	116	ivec RV::dataind ( RV rv2 ) const {
	117	str str2 = rv2.tostr();
	118	ivec res ( 0 );
	119	ivec part;
	120	int i;
	121	for ( i = 0;i < len;i++ ) {
	122	part = itpp::find ( ( str2.ids == ids ( i ) ) & ( str2.times == times ( i ) ) );
	123	res = concat ( res, part );
	124	}
	125	return res;
	126	}
	127
	128	RV RV::subt ( const RV rv2 ) const {
	129	ivec res = this->findself ( rv2 ); // nonzeros
	130	ivec valid = itpp::find ( res == -1 ); //-1 => value not found => it remains
	131	return ( *this ) ( valid ); //keep those that were not found in rv2
	132	}
	133
	134	ivec RV::findself ( const RV &rv2 ) const {
	135	int i, j;
	136	ivec tmp = -ones_i ( len );
	137	for ( i = 0;i < len;i++ ) {
	138	for ( j = 0;j < rv2.length();j++ ) {
	139	if ( ( ids ( i ) == rv2.ids ( j ) ) & ( times ( i ) == rv2.times ( j ) ) ) {
	140	tmp ( i ) = j;
	141	break;
	142	}
	143	}
	144	}
	145	return tmp;
98	146	}
99	147

bdm/stat/libBM.h

r118	r145
18	18
19	19	using namespace itpp;
	20
	21	//! Structure of RV (used internally)
	22	class str{
	23	public:
	24	ivec ids;
	25	ivec times;
	26	str(ivec ids0, ivec times0):ids(ids0),times(times0){
	27	it_assert_debug(times0.length()==ids0.length(),"Incompatible input");
	28	};
	29	};
20	30
21	31	/*!
…	…
61	71	//TODO why not inline and later??
62	72
63		//! Find indexes of ~~another rv in self~~
64		ivec find ~~( RV rv2 )~~;
	73	//! Find indexes of self in another rv, \return ivec of the same size as self.
	74	ivec findself (const RV &rv2 ) const;
65	75	//! Compare if \c rv2 is identical to this \c RV
66		bool equal ( RV rv2 ) const;
67		//! Add (concat) another variable to the current one
68		void add ( const RV &rv2 );
69		//! Add (concat) another variable to the current one
70		friend RV concat ( const RV &rv1, const RV &rv2 );
	76	bool equal (const RV &rv2 ) const;
	77	//! Add (concat) another variable to the current one, \return 0 if all rv2 were added, 1 if rv2 is in conflict
	78	bool add ( const RV &rv2 );
71	79	//! Subtract another variable from the current one
72		RV subt ( ~~RV rv2 )~~;
	80	RV subt ( const RV rv2 ) const;
73	81	//! Select only variables at indeces ind
74		RV subselect ( ivec ind );
	82	RV subselect ( ivec ind ) const;
75	83	//! Select only variables at indeces ind
76		RV operator() ( ivec ind );
77		//! ~~Generate new \c RV with~~ \c time shifted by delta.
	84	RV operator() ( ivec ind ) const;
	85	//! Shift \c time shifted by delta.
78	86	void t ( int delta );
79		//! generate a list of indeces, i.e. which
80		ivec indexlist();
	87	//! generate \c str from rv, by expanding sizes
	88	str tostr() const;
	89	//! generate indeces into \param crv data vector that form data vector of self.
	90	ivec dataind(RV crv) const;
81	91
82	92	//!access function
…	…
92	102	std::string name ( int at ) {return names ( at );};
93	103	};
	104
	105	//! Concat two random variables
	106	RV concat ( const RV &rv1, const RV &rv2 );
94	107
95	108
…	…
146	159	//! Destructor for future use;
147	160	virtual ~epdf() {};
148		//! access function
149		RV ~~_rv() const~~ {return rv;}
	161	//! access function, possibly dangerous!
	162	RV& _rv() {return rv;}
150	163	};
151	164
…	…
170	183	vec temp= ep->sample();
171	184	ll=ep->evalpdflog ( temp );return temp;};
172		//! Returns N samples from the density conditioned on \c cond, \f$x \sim epdf(rv\|cond)\f$. \param cond is numeric value of \c rv \param ll is a return value of log-likelihood of the sample.
	185	//! Returns \param N samples from the density conditioned on \c cond, \f$x \sim epdf(rv\|cond)\f$. \param cond is numeric value of \c rv \param ll is a return value of log-likelihood of the sample.
173	186	virtual mat samplecond ( vec &cond, vec &ll, int N ) {
174	187	this->condition ( cond );
…	…
190	203	//! access function
191	204	RV _rvc() {return rvc;}
	205	//! access function
	206	RV _rv() {return rv;}
192	207	//!access function
193	208	epdf& _epdf() {return *ep;}
…	…
201	216	public:
202	217	//!Default constructor
203		mepdf ( ~~const RV &rv, const RV &rvc, epdf* em ) :mpdf ( rv,rvc ) {ep=~~em;};
	218	mepdf (epdf &em ) :mpdf ( em._rv(),RV() ) {ep=&em;};
204	219	};
205	220

bdm/stat/loggers.h

r102	r145
84	84	void step(bool final=false) {if ( ind<maxlen ) ind++; else it_error ( "memlog::ind is too high;" );}
85	85	void logit ( int id, vec v ) {vectors ( id ).set_row ( ind,v );}
	86	//! Save values into an itfile named after \c fname.
86	87	void itsave(const char* fname);
87	88	};

Context Navigation

Changeset 145 for bdm

Legend:

Download in other formats: