Context Navigation

← Previous Changeset
Next Changeset →

Changeset 213

Timestamp:

11/25/08 12:24:42 (15 years ago)

Author:

smidl

Message:

Merging - new experiment

Files:

: 6 modified

CMakeLists.txt (modified) (1 diff)
bdm/estim/merger.cpp (modified) (6 diffs)
mpdm/CMakeLists.txt (modified) (1 diff)
mpdm/merg_2a.cpp (modified) (3 diffs)
mpdm/merg_pred.cpp (modified) (6 diffs)
tests/merger_2d_test.cpp (modified) (1 diff)

Legend:

: Unmodified
: Added
: Removed

CMakeLists.txt

r211	r213
27	27	#This is for UNIX makefile which does only one release at a time.
28	28	SET(CMAKE_BUILD_TYPE Debug)
29		SET(CMAKE_CXX_FLAGS_DEBUG "-Wall -g -DNDEBUG -O3 -pipe -Wno-unknown-pragmas")
	29	# SET(CMAKE_CXX_FLAGS_DEBUG "-Wall -g -DNDEBUG -O3 -pipe -Wno-unknown-pragmas")
	30	SET(CMAKE_CXX_FLAGS_DEBUG "-Wall -g -pipe -Wno-unknown-pragmas")
30	31	INCLUDE(CMakeLocal.txt OPTIONAL)
31	32	ENDIF(WIN32)

bdm/estim/merger.cpp

r211	r213
6	6	int nu=lW.rows();
7	7	vec mu = sum ( lW ) /nu; //mean of logs
8		vec lam = sum ( pow ( lW,2 ) )-nu*pow ( mu,2 );
	8	// vec lam = sum ( pow ( lW,2 ) )-nu*pow ( mu,2 ); ======= numerically unsafe!
	9	vec lam = sum ( pow ( lW-outer_product ( ones ( lW.rows() ),mu ),2 ) );
9	10	double coef=0.0;
10	11	vec sq2bl=sqrt ( 2betalam ); //this term is everywhere
…	…
61	62	char str[100];
62	63
63		epdf* Mpred=Mix.predictor~~(rv~~);
	64	epdf* Mpred=Mix.predictor ( rv );
64	65	vec Mix_pdf ( Ns );
65	66	while ( !converged ) {
…	…
70	71	delete Mpred;
71	72	Mpred = Mix.predictor ( rv ); // Allocation => must be deleted at the end!!
72
	73
	74	// This will be active only later in iterations!!!
73	75	if ( 1./sum_sqr ( w ) <0.5*Ns ) {
74	76	// Generate new samples
…	…
76	78	for ( int i=0;i<Ns;i++ ) {
77	79	//////////// !!!!!!!!!!!!!
78		if ( Smp ( i ) ( ~~1 ) <0 ) {Smp ( i ) ( 1~~ ) = 0.01; }
	80	if ( Smp ( i ) ( 2 ) <0 ) {Smp ( i ) ( 2 ) = 0.01; }
79	81	set_col_part ( Smp_ex,i,Smp ( i ) );
80	82	}
81		~~{cout<<"Resampling =" << 1./sum_sqr ( w ) << endl;}~~
	83	if(0){cout<<"Resampling =" << 1./sum_sqr ( w ) << endl;
82	84	cout << sum ( Smp_ex,2 ) /Ns <<endl;
83		cout << Smp_ex*Smp_ex.T()~~/Ns << endl;~~
	85	cout << Smp_ex*Smp_ex.T() /Ns << endl;}
84	86	}
85	87	// sprintf ( str,"Mpred_mean%d",niter );
…	…
172	174	}
173	175	else {
174
	176	it_file itf ( "merg_err.it" );
	177	itf << Name ( "w" ) << w;
175	178	}
176	179
…	…
188	191	converged = ( niter>20 );
189	192	}
	193	delete Mpred;
190	194	cout << endl;
191	195

mpdm/CMakeLists.txt

r204	r213
7	7	## Save all needed libraries in variable BdmLibs
8	8	#SET(BdmLibs bdm itpp_debug)
9		~~SET(BdmLibs bdm itpp)~~
10		~~SET(ITppLibs itpp)~~
11	9
12	10	IF(WIN32)
13		SET(BdmLibs ${BdmLibs} libacml_dll)
14		SET(ITppLibs ${ITppLibs} libacml_dll)
	11	SET(AddLib libacml_dll)
15	12	ENDIF(WIN32)
	13
	14	# Define macro for testing a file
	15	MACRO(TEST FN)
	16	add_executable (${FN} ${FN}.cpp)
	17	target_link_libraries (${FN} debug itpp_debug)
	18	target_link_libraries (${FN} optimized itpp)
	19	target_link_libraries (${FN} bdm ${AddLib})
	20	ENDMACRO(TEST)
16	21
17	22	# Add executable called "helloDemo" that is built from the source files
18	23	# "demo.cxx" and "demo_b.cxx". The extensions are automatically found.
19		add_executable (merg_pred merg_pred.cpp)
20		target_link_libraries (merg_pred ${BdmLibs})
21
22		add_executable (merg_giw merg_giw.cpp)
23		target_link_libraries (merg_giw ${BdmLibs})
24
25		add_executable (merger_iter_cond merger_iter_cond.cpp)
26		target_link_libraries (merger_iter_cond ${BdmLibs})
	24	TEST(merg_pred)
	25	TEST(merg_giw)
	26	TEST(merger_iter_cond)
	27	TEST(merg_2a)

mpdm/merg_2a.cpp

r211	r213
10	10	using std::endl;
11	11
	12	/*! The purpose of this experiment is to test merging fragmental pdfs between two participants sharing the same parameter, a. However, this parameter has a different role in each model.
	13
	14	Lets assume that:
	15	u -> \|a,r\| -> y -> \|a,r\| -> z
	16	*/
	17
12	18	int main() {
13	19	// Setup model
14	20	RV y ( "{y }" );
15	21	RV u ( "{u }" );
16		~~RV um = u; um.t(-1);~~
17	22	RV z ( "{z }" );
18		RV a ("{a }");
19		RV b ("{b }");
20		RV c ("{c }");
	23	RV a ("{a }");
	24	RV b ("{b }"); RV ab = a; ab.add(b);
	25	RV c ("{c }"); RV ac = a; ac.add(c);
21	26	RV r ("{r }");
22	27
23	28	double at = 1.5;
24		double bt = 0.8;
25		double ct = ~~0.50~~;
26		double ~~sig = 0.1~~0;
	29	double bt = 0.5;
	30	double ct = -0.5;
	31	double rt = 0.30;
27	32	// Full system
28	33	vec thy =vec_2 ( at,bt ); //Simulated system - zero for constant term
29		vec thu =vec_2 ( at,ct ); //Simulated system - zero for constant term
30		vec Thy = concat ( thy, vec_1(~~sig~~) ); //Full parameter
31		vec Th~~u = concat ( thu, vec_1(sig~~) ); //Full parameter
	34	vec thz =vec_2 ( at,ct ); //Simulated system - zero for constant term
	35	vec Thy = concat ( thy, vec_1(rt) ); //Full parameter
	36	vec Thz = concat ( thz, vec_1(rt) ); //Full parameter
32	37
33	38	//ARX constructor
34	39	mat V0 = 0.001*eye ( 3 ); V0 ( 0,0 ) = 1; //
35	40
36		ARX P1 ( concat ( a~~, concat(b,r)~~ ), V0, -1 );
37		ARX P2 ( concat ( a~~, concat(c,r)~~ ), V0, -1 );
	41	ARX P1 ( concat ( ab,r ), V0, -1 );
	42	ARX P2 ( concat ( ac,r ), V0, -1 );
38	43
39	44	//Test estimation
…	…
42	47
43	48	// Logging
44		dirfilelog L ( "exp/merg_2a",ndat );
45		int Li_Data = L.add ( RV ( "{Y U1 U2 }" ), "" );
46		int Li_LL = L.add ( RV ( "{G M }" ), "LL" );
47		int Li_P1m = L.add ( RV ( "{a b r }" ), "P1" );
48		int Li_P2m = L.add ( RV ( "{a c r }" ), "P2" );
	49	dirfilelog L ( "exp/merg_2a",3 );
	50	int Li_Data = L.add ( RV ( "{U Y Z }" ), "" );
	51	// int Li_LL = L.add ( RV ( "{P1 P2 M1 M2 }" ), "LL" );
	52	int Li_P1m = L.add ( concat ( ab,r ), "P1m" );
	53	int Li_P2m = L.add ( concat ( ac,r ), "P2m" );
	54	int Li_Mm = L.add ( concat ( ab,concat(r,c) ), "Mm" );
49	55	L.init();
50	56
	57	vec Ut ( ndat );
51	58	vec Yt ( ndat );
	59	vec Zt ( ndat );
52	60	vec yt ( 1 );
53	61
54	62	//Proposal
55		enorm<ldmat> g0 ( a~~1 ); g0.set_parameters ( "1 ",mat ( "~~1" ) );
	63	enorm<ldmat> g0 ( ab ); g0.set_parameters ( "1 1 ",mat ( "1 0; 0 1" ) );
56	64	egamma g1 ( r ); g1.set_parameters ( "2 ", "2" );
57		enorm<ldmat> g2 ( ~~a2 ); g2.set_parameters ( "1~~ ",mat ( "1" ) );
	65	enorm<ldmat> g2 ( c ); g2.set_parameters ( "1 ",mat ( "1" ) );
58	66
59		Array<const epdf*> A ( 3 ); A ( 0 ) = &g0; A ( 1 ) =&g1; A ~~( 2 ) = &g2;~~
	67	Array<const epdf*> A ( 3 ); A ( 0 ) = &g0; A ( 1 ) =&g1; A(2) = &g2;
60	68	eprod G0 ( A );
61	69
62		for ( t=0; t<ndat; t++ ) {
	70	vec rgru(2);
	71	vec rgry(2);
	72	Yt(0) = 0.1;
	73	Ut(0) = 0.0;
	74	for ( t=1; t<ndat; t++ ) {
63	75	// True system
64		rgrg ( 0 ) = pow ( sin ( ( t/40.0 ) *pi ),3 );
65		rgrg ( 1 ) = pow ( cos ( ( t/40.0 +0.1 ) *pi ),3 );
66
67		Yt ( t ) = thgrgrg + sqr NorRNG();
	76	Ut ( t ) = pow ( sin ( ( t/40.0 ) *pi ),3 );
	77	rgru(0) = Ut(t); rgru(1) = Ut(t-1);
	78	Yt ( t ) = thyrgru + rt NorRNG();
	79	rgry(0) = Yt(t); rgry(1) = Yt(t-1);
	80	Zt ( t ) = thzrgry + rt NorRNG();
68	81
69	82	// Bayes for all
70		P1.bayes ( concat ( Yt ( t ),vec_1 ( rgrg ( 0 ) ) ) );
71		P2.bayes ( concat ( Yt ( t ),vec_1 ( rgrg ( 1 ) ) ) );
72		PG.bayes ( concat ( Yt ( t ),rgrg ) );
73
74
75		// crippling PGk by substituting zeros.
76		/* ldmat &Vk=const_cast<egiw*>(PGk._e())->_V(); //PG ldmat does not like 0!
77		mat fVk=PG._e()->_V().to_mat();
78		fVk(1,2) = 0.0;
79		fVk(2,1) = 0.0;
80		Vk = ldmat(fVk);
81		*/ //PGk is now krippled
	83	P1.bayes ( concat ( Yt ( t ),rgru ) );
	84	P2.bayes ( concat ( Zt ( t ),rgry ) );
82	85
83	86	// Merge estimates
…	…
89	92	//M2.set_parameters ( 100.0, 1000,3 ); //M2._Mix().set_method(QB);
90	93	M.merge ( &G0 );
91		~~//M2.merge ( &G0 );~~
92	94
93	95	//Likelihood
94	96	yt ( 0 ) = Yt ( t );
95	97
96		~~// LLs ( 0 ) = P1._e()->evalpdflog ( get_vec(Th, "1 2") );~~
97		~~// LLs ( 1 ) = P2._e()->evalpdflog ( get_vec(Th, "3 2") );~~
98		~~LLs ( 0 ) = PG._e()->evalpdflog ( Th );~~
99		~~LLs ( 1 ) = M._Mix().logpred ( concat ( Thj, vec_1 ( 1.0 ) ) );~~
100		~~// LLs ( 3 ) = M2._Mix().logpred ( concat(Thj, vec_1(1.0)) );~~
101		~~L.logit ( Li_LL, LLs ); //log-normal~~
102
103	98	//Logger
104		L.logit ( Li_Data, vec_3 ( Yt ( t ), rgrg ( 0 ), rgrg ( 1 ) ) );
105		L.logit ( Li_P1m, P1._e()->mean() );
106		L.logit ( Li_P2m, P2._e()->mean() );
107		L.logit ( Li_Gm, PG._e()->mean() );
108		L.logit ( Li_Mm, M.mean() );
109
	99	L.logit(Li_Data, vec_3(Ut(t), Yt(t), Zt(t)));
	100	L.logit(Li_P1m, P1._e()->mean());
	101	L.logit(Li_P2m, P2._e()->mean());
	102	L.logit(Li_Mm, M.mean());
110	103	L.step ( );
111	104
112		~~cout << "Vg: " << PG._e()->_V().to_mat() <<endl;~~
113		~~vec mea = M.mean();~~
114		~~cout << "Ve: " << M.variance() <<endl;~~
115	105	}
116	106	L.finalize( );
117		L.itsave ( "merg_~~egiw~~.it" );
	107	L.itsave ( "merg_2a.it" );
118	108	cout << endl;
119	109	}

mpdm/merg_pred.cpp

r211	r213
15	15	RV u1 ( "{u1 }" );
16	16	RV u2 ( "{u2 }" );
17		RV ym=y; ym.t~~(-1~~);
18		RV yy = y; yy.add~~(ym~~);
19		RV uu=u1; uu.add~~(u2~~);
	17	RV ym=y; ym.t ( -1 );
	18	RV yy = y; yy.add ( ym );
	19	RV uu=u1; uu.add ( u2 );
20	20
21	21	// Full system
…	…
44	44
45	45	// Logging
46		dirfilelog L ( "exp/merg",~~ndat~~ );
	46	dirfilelog L ( "exp/merg",10 );
47	47
48	48	int Li_Eth1 = L.add ( thri,"P1" );
…	…
56	56
57	57	vec Yt ( ndat );
58		vec yt(1);
	58	vec yt ( 1 );
59	59
60	60	Yt.set_subvector ( 0,randn ( ord ) ); //initial values
61		vec rgrg ( thrg.count() -1); // constant terms are in!
62		vec rgr1 ( thri.count() -1);
63		vec rgr2 ( thri.count() -1);
	61	vec rgrg ( thrg.count() -1 ); // constant terms are in!
	62	vec rgr1 ( thri.count() -1 );
	63	vec rgr2 ( thri.count() -1 );
64	64
65		vec PredLLs(5);
66		vec PostLLs(3);
67		vec PredLLs_m=zeros(5);
	65	vec PredLLs ( 5 );
	66	vec PostLLs ( 3 );
	67	vec PredLLs_m=zeros ( 5 );
68	68	ivec ind_r1 = "0 1 3";
69	69	ivec ind_r2 = "0 2 3";
…	…
73	73	rgrg ( 0 ) =Yt ( t-1 );
74	74	rgrg ( 1 ) = pow(sin ( ( t/40.0 ) *pi ),3);
75		rgrg ( 2 ) = pow(cos ( ( t/40.0 ) *pi ),3);
76		rgrg (3) = 1.0; // constant term
77
78		rgr1~~(0) = rgrg(0); rgr1(1) = rgrg(1); rgr1(2) = rgrg(3~~); // no u2
79		rgr2~~(0) = rgrg(0); rgr2(1) = rgrg(2); rgr2(2) = rgrg(3~~); // no u1
	75	rgrg ( 2 ) = pow(cos ( ( t/60.0 ) *pi ),3);
	76	rgrg ( 3 ) = 1.0; // constant term
	77
	78	rgr1 ( 0 ) = rgrg ( 0 ); rgr1 ( 1 ) = rgrg ( 1 ); rgr1 ( 2 ) = rgrg ( 3 ); // no u2
	79	rgr2 ( 0 ) = rgrg ( 0 ); rgr2 ( 1 ) = rgrg ( 2 ); rgr2 ( 2 ) = rgrg ( 3 ); // no u1
80	80
81	81	Yt ( t ) = thgrgrg + sqr NorRNG();
82	82
83	83	// Test predictors
84		if (t>2){
85		mlnorm<ldmat>* P1p = P1.predictor(y,concat(ym,u1));
86		mlnorm<ldmat>* P2p = P2.predictor(y,concat(ym,u2));
87		mlnorm<ldmat>* Pgp = PG.predictor(y,concat(ym,uu));
88
89		Array<mpdf*> A(2); A(0)=P1p;A(1)=P2p;
90		merger M(A);
91		enorm<ldmat> g0(concat(yy,uu)); g0.set_parameters("0 0 0 0 ",3*eye(4));
92		M.set_parameters(10000000.0, 100,1);
93		M.merge(&g0);
94
95		yt(0) = Yt(t);
96		double P1pl = P1p->evallogcond(yt,rgr1);
97		double P2pl = P2p->evallogcond(yt,rgr2);
98		double PGpl = Pgp->evallogcond(yt,rgrg);
	84	if ( t>2 ) {
	85	mlnorm<ldmat>* P1p = P1.predictor ( y,concat ( ym,u1 ) );
	86	mlnorm<ldmat>* P2p = P2.predictor ( y,concat ( ym,u2 ) );
	87	mlnorm<ldmat>* Pgp = PG.predictor ( y,concat ( ym,uu ) );
	88
	89	Array<mpdf*> A ( 2 ); A ( 0 ) =P1p;A ( 1 ) =P2p;
	90	merger M ( A );
	91	enorm<ldmat> g0 ( concat ( yy,uu ) ); g0.set_parameters ( "0 0 0 0 ",3*eye ( 4 ) );
	92	M.set_parameters ( 1e8, 101,1 );
	93	M.merge ( &g0 );
	94	yt ( 0 ) = Yt ( t );
	95	double P1pl = P1p->evallogcond ( yt,rgr1 );
	96	double P2pl = P2p->evallogcond ( yt,rgr2 );
	97	double PGpl = Pgp->evallogcond ( yt,rgrg );
99	98	{
100		cout << "yt: " << yt << endl;
	99	cout << "T: " << t << "yt: " << yt << endl;
101	100	cout << "yt_1: " << P1p->_epdf().mean() << endl;
	101	cout << *P1p <<endl;
102	102	cout << "yt_2: " << P2p->_epdf().mean() << endl;
103	103	cout << "yt_G: " << P2p->_epdf().mean() << endl;
…	…
106	106	double cP2pl;
107	107	{
108		ARX* Apred = (~~ARX*)M._Mix()._Coms(0~~);
109		enorm<ldmat>* MP= Apred->predictor~~(concat(yy,uu)~~);
110		enorm<ldmat>* mP1p = (~~enorm<ldmat>*)MP->marginal(concat(yy,u1)~~);
111		enorm<ldmat>* mP2p = (~~enorm<ldmat>*)MP->marginal(concat(yy,u2)~~);
112		mlnorm<ldmat>* cP1p = (~~mlnorm<ldmat>*)mP1p->condition(y~~);
113		mlnorm<ldmat>* cP2p = (~~mlnorm<ldmat>*)mP2p->condition(y~~);
	108	ARX* Apred = ( ARX* ) M._Mix()._Coms ( 0 );
	109	enorm<ldmat>* MP= Apred->predictor ( concat ( yy,uu ) );
	110	enorm<ldmat>* mP1p = ( enorm<ldmat>* ) MP->marginal ( concat ( yy,u1 ) );
	111	enorm<ldmat>* mP2p = ( enorm<ldmat>* ) MP->marginal ( concat ( yy,u2 ) );
	112	mlnorm<ldmat>* cP1p = ( mlnorm<ldmat>* ) mP1p->condition ( y );
	113	mlnorm<ldmat>* cP2p = ( mlnorm<ldmat>* ) mP2p->condition ( y );
114	114
115		cP1pl = cP1p->evallogcond~~(yt,rgr1);~~
116		cP2pl = cP2p->evallogcond~~(yt,rgr2);~~
117
	115	cP1pl = cP1p->evallogcond ( yt,rgr1 );
	116	cP2pl = cP2p->evallogcond ( yt,rgr2 );
	117
118	118	cout << "ytm1: " << cP1p->_epdf().mean() << endl;
119	119	cout << "ytm2: " << cP2p->_epdf().mean() << endl;
	120	cout << *cP2p << endl;
120	121	}
121	122
122	123	PredLLs *=frg;
123		PredLLs += ~~log(concat(vec_3(P1pl, P2pl, PGpl), vec_2(cP1pl,cP2pl))~~);
124		L.logit~~(Li_Pred, PredLLs~~); //log-normal
125
	124	PredLLs += concat ( vec_3 ( P1pl, P2pl, PGpl ), vec_2 ( cP1pl,cP2pl ) );
	125	L.logit ( Li_Pred, PredLLs ); //log-normal
	126
126	127	delete P1p;
127	128	delete P2p;
128	129	delete Pgp;
129	130	}
130
131		~~// 1st~~
132		~~P1.bayes ( concat(Yt(t),rgr1) );~~
133		~~// 2nd~~
134		~~P2.bayes ( concat(Yt(t),rgr2) );~~
135	131
136		//Global
137		PG.bayes ( concat ( Yt ( t ),rgrg ) );
138
	132	if ( t<100 ) {
	133	// 1st
	134	P1.bayes ( concat ( Yt ( t ),rgr1 ) );
	135	// 2nd
	136	P2.bayes ( concat ( Yt ( t ),rgr2 ) );
	137
	138	//Global
	139	PG.bayes ( concat ( Yt ( t ),rgrg ) );
	140	}
139	141	//Merger
140	142	}
…	…
146	148	PostLLs += vec_3 ( P1._ll(), P2._ll(), PG._ll() );
147	149	L.logit ( Li_LL, PostLLs );
148		L.step ( );
	150	L.step ( );
149	151	}
150	152	L.finalize( );

tests/merger_2d_test.cpp

r211	r213
44	44	merger M ( A );
45	45	enorm<fsqmat> g0(xy);
46		g0.set_parameters(vec("1 1"),mat("1 ~~0; 0 1~~"));
	46	g0.set_parameters(vec("1 1"),mat("10 0; 0 10"));
47	47
48	48	M.set_parameters(1e8,400,1);

Download in other formats: